! $Id: fisrtilp_tr.f90 5285 2024-10-28 13:33:29Z abarral $ SUBROUTINE fisrtilp_tr(dtime, paprs, pplay, t, q, ratqs, d_t, d_q, d_ql, & rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, & frac_nucl, prfl, psfl, rhcl) ! relative humidity in clear sky (needed for aer optical ! properties; aeropt.F) USE dimphy USE print_control_mod, ONLY: lunout USE yomcst_mod_h USE yoethf_mod_h IMPLICIT NONE ! ====================================================================== ! Auteur(s): Z.X. Li (LMD/CNRS) ! Date: le 20 mars 1995 ! Objet: condensation et precipitation stratiforme. ! schema de nuage ! ====================================================================== ! ====================================================================== ! Arguments: REAL dtime ! intervalle du temps (s) REAL paprs(klon, klev+1) ! pression a inter-couche REAL pplay(klon, klev) ! pression au milieu de couche REAL t(klon, klev) ! temperature (K) REAL q(klon, klev) ! humidite specifique (kg/kg) REAL d_t(klon, klev) ! incrementation de la temperature (K) REAL d_q(klon, klev) ! incrementation de la vapeur d'eau REAL d_ql(klon, klev) ! incrementation de l'eau liquide REAL rneb(klon, klev) ! fraction nuageuse REAL radliq(klon, klev) ! eau liquide utilisee dans rayonnements REAL rain(klon) ! pluies (mm/s) REAL snow(klon) ! neige (mm/s) REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s) ! jq For aerosol opt properties needed (see aeropt.F) REAL rhcl(klon, klev) ! AA ! Coeffients de fraction lessivee : pour OFF-LINE REAL pfrac_nucl(klon, klev) REAL pfrac_1nucl(klon, klev) REAL pfrac_impa(klon, klev) ! Fraction d'aerosols lessivee par impaction et par nucleation ! POur ON-LINE REAL frac_impa(klon, klev) REAL frac_nucl(klon, klev) ! AA ! Options du programme: REAL seuil_neb ! un nuage existe vraiment au-dela PARAMETER (seuil_neb=0.001) REAL ct ! inverse du temps pour qu'un nuage precipite PARAMETER (ct=1./1800.) REAL cl ! seuil de precipitation PARAMETER (cl=2.6E-4) ! cc PARAMETER (cl=2.3e-4) ! cc PARAMETER (cl=2.0e-4) INTEGER ninter ! sous-intervals pour la precipitation PARAMETER (ninter=5) LOGICAL evap_prec ! evaporation de la pluie PARAMETER (evap_prec=.TRUE.) REAL coef_eva PARAMETER (coef_eva=2.0E-05) LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur PARAMETER (calcrat=.TRUE.) REAL zx_min, rat_max PARAMETER (zx_min=1.0, rat_max=0.01) REAL zx_max, rat_min PARAMETER (zx_max=0.1, rat_min=0.3) REAL zx LOGICAL cpartiel ! condensation partielle PARAMETER (cpartiel=.TRUE.) REAL t_coup PARAMETER (t_coup=234.0) ! Variables locales: INTEGER i, k, n, kk REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 REAL zrfl(klon), zrfln(klon), zqev, zqevt REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq REAL ztglace, zt(klon) INTEGER nexpo ! exponentiel pour glace/eau REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon) REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon) LOGICAL appel1er SAVE appel1er !$OMP THREADPRIVATE(appel1er) ! --------------------------------------------------------------- ! AA Variables traceurs: ! AA Provisoire !!! Parametres alpha du lessivage ! AA A priori on a 4 scavenging # possibles REAL a_tr_sca(4) SAVE a_tr_sca !$OMP THREADPRIVATE(a_tr_sca) ! Variables intermediaires REAL zalpha_tr REAL zfrac_lessi REAL zprec_cond(klon) ! AA ! --------------------------------------------------------------- ! Fonctions en ligne: REAL fallv ! vitesse de chute pour crystaux de glace REAL zzz include "FCTTRE.h" fallv(zzz) = 3.29/2.0*((zzz)**0.16) ! cc fallv (zzz) = 3.29/3.0 * ((zzz)**0.16) ! cc fallv (zzz) = 3.29 * ((zzz)**0.16) DATA appel1er/.TRUE./ IF (appel1er) THEN WRITE (lunout, *) 'fisrtilp, calcrat:', calcrat WRITE (lunout, *) 'fisrtilp, ninter:', ninter WRITE (lunout, *) 'fisrtilp, evap_prec:', evap_prec WRITE (lunout, *) 'fisrtilp, cpartiel:', cpartiel IF (abs(dtime/real(ninter)-360.0)>0.001) THEN WRITE (lunout, *) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime WRITE (lunout, *) 'Je prefere un sous-intervalle de 6 minutes' CALL abort END IF appel1er = .FALSE. ! AA initialiation provisoire a_tr_sca(1) = -0.5 a_tr_sca(2) = -0.5 a_tr_sca(3) = -0.5 a_tr_sca(4) = -0.5 ! AA Initialisation a 1 des coefs des fractions lessivees DO k = 1, klev DO i = 1, klon pfrac_nucl(i, k) = 1. pfrac_1nucl(i, k) = 1. pfrac_impa(i, k) = 1. END DO END DO END IF ! test sur appel1er ! MAf Initialisation a 0 de zoliq DO i = 1, klon zoliq(i) = 0. END DO ! Determiner les nuages froids par leur temperature ztglace = rtt - 15.0 nexpo = 6 ! cc nexpo = 1 ! Initialiser les sorties: DO k = 1, klev + 1 DO i = 1, klon prfl(i, k) = 0.0 psfl(i, k) = 0.0 END DO END DO DO k = 1, klev DO i = 1, klon d_t(i, k) = 0.0 d_q(i, k) = 0.0 d_ql(i, k) = 0.0 rneb(i, k) = 0.0 radliq(i, k) = 0.0 frac_nucl(i, k) = 1. frac_impa(i, k) = 1. END DO END DO DO i = 1, klon rain(i) = 0.0 snow(i) = 0.0 END DO ! Initialiser le flux de precipitation a zero DO i = 1, klon zrfl(i) = 0.0 zneb(i) = seuil_neb END DO ! AA Pour plus de securite zalpha_tr = 0. zfrac_lessi = 0. ! AA---------------------------------------------------------- ! Boucle verticale (du haut vers le bas) DO k = klev, 1, -1 ! AA---------------------------------------------------------- DO i = 1, klon zt(i) = t(i, k) zq(i) = q(i, k) END DO ! Calculer l'evaporation de la precipitation IF (evap_prec) THEN DO i = 1, klon IF (zrfl(i)>0.) THEN IF (thermcep) THEN zdelta = max(0., sign(1.,rtt-zt(i))) zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k) zqs(i) = min(0.5, zqs(i)) zcor = 1./(1.-retv*zqs(i)) zqs(i) = zqs(i)*zcor ELSE IF (zt(i)<t_coup) THEN zqs(i) = qsats(zt(i))/pplay(i, k) ELSE zqs(i) = qsatl(zt(i))/pplay(i, k) END IF END IF zqev = max(0.0, (zqs(i)-zq(i))*zneb(i)) zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* & (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*zt(i)*rd/rg zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ & (paprs(i,k)-paprs(i,k+1)) zqev = min(zqev, zqevt) zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, & k+1)))*dtime zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, & k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i)) zrfl(i) = zrfln(i) END IF END DO END IF ! Calculer Qs et L/Cp*dQs/dT: IF (thermcep) THEN DO i = 1, klon zdelta = max(0., sign(1.,rtt-zt(i))) zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i)) zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k) zqs(i) = min(0.5, zqs(i)) zcor = 1./(1.-retv*zqs(i)) zqs(i) = zqs(i)*zcor zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor) END DO ELSE DO i = 1, klon IF (zt(i)<t_coup) THEN zqs(i) = qsats(zt(i))/pplay(i, k) zdqs(i) = dqsats(zt(i), zqs(i)) ELSE zqs(i) = qsatl(zt(i))/pplay(i, k) zdqs(i) = dqsatl(zt(i), zqs(i)) END IF END DO END IF ! Determiner la condensation partielle et calculer la quantite ! de l'eau condensee: IF (cpartiel) THEN DO i = 1, klon zdelq = ratqs(i, k)*zq(i) rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq) zqn(i) = (zq(i)+zdelq+zqs(i))/2.0 IF (rneb(i,k)<=0.0) zqn(i) = 0.0 IF (rneb(i,k)>=1.0) zqn(i) = zq(i) rneb(i, k) = max(0.0, min(1.0,rneb(i,k))) zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)/(1.+zdqs(i)) ! --Olivier rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i) IF (rneb(i,k)<=0.0) rhcl(i, k) = zq(i)/zqs(i) IF (rneb(i,k)>=1.0) rhcl(i, k) = 1.0 ! --fin END DO ELSE DO i = 1, klon IF (zq(i)>zqs(i)) THEN rneb(i, k) = 1.0 ELSE rneb(i, k) = 0.0 END IF zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i)) END DO END IF DO i = 1, klon zq(i) = zq(i) - zcond(i) zt(i) = zt(i) + zcond(i)*rlvtt/rcpd END DO ! Partager l'eau condensee en precipitation et eau liquide nuageuse DO i = 1, klon IF (rneb(i,k)>0.0) THEN zoliq(i) = zcond(i) zrho(i) = pplay(i, k)/zt(i)/rd zdz(i) = (paprs(i,k)-paprs(i,k+1))/(zrho(i)*rg) zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace) zfice(i) = min(max(zfice(i),0.0), 1.0) zfice(i) = zfice(i)**nexpo zneb(i) = max(rneb(i,k), seuil_neb) radliq(i, k) = zoliq(i)/real(ninter+1) END IF END DO DO n = 1, ninter DO i = 1, klon IF (rneb(i,k)>0.0) THEN zchau(i) = ct*dtime/real(ninter)*zoliq(i)* & (1.0-exp(-(zoliq(i)/zneb(i)/cl)**2))*(1.-zfice(i)) zrhol(i) = zrho(i)*zoliq(i)/zneb(i) zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)*fallv(zrhol(i))* & zfice(i) ztot(i) = zchau(i) + zfroi(i) IF (zneb(i)==seuil_neb) ztot(i) = 0.0 ztot(i) = min(max(ztot(i),0.0), zoliq(i)) zoliq(i) = max(zoliq(i)-ztot(i), 0.0) radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1) END IF END DO END DO DO i = 1, klon IF (rneb(i,k)>0.0) THEN d_ql(i, k) = zoliq(i) zrfl(i) = zrfl(i) + max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k & +1))/(rg*dtime) END IF IF (zt(i)<rtt) THEN psfl(i, k) = zrfl(i) ELSE prfl(i, k) = zrfl(i) END IF END DO ! Calculer les tendances de q et de t: DO i = 1, klon d_q(i, k) = zq(i) - q(i, k) d_t(i, k) = zt(i) - t(i, k) END DO ! AA--------------- Calcul du lessivage stratiforme ------------- DO i = 1, klon zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k+1))/ & rg IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN ! AA lessivage nucleation LMD5 dans la couche elle-meme IF (t(i,k)>=ztglace) THEN zalpha_tr = a_tr_sca(3) ELSE zalpha_tr = a_tr_sca(4) END IF zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i)) pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi) frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi ! nucleation avec un facteur -1 au lieu de -0.5 zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i)) pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi) END IF END DO ! boucle sur i ! AA Lessivage par impaction dans les couches en-dessous DO kk = k - 1, 1, -1 DO i = 1, klon IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN IF (t(i,kk)>=ztglace) THEN zalpha_tr = a_tr_sca(1) ELSE zalpha_tr = a_tr_sca(2) END IF zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i)) pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi) frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi END IF END DO END DO ! AA---------------------------------------------------------- ! FIN DE BOUCLE SUR K END DO ! AA----------------------------------------------------------- ! Pluie ou neige au sol selon la temperature de la 1ere couche DO i = 1, klon IF ((t(i,1)+d_t(i,1))<rtt) THEN snow(i) = zrfl(i) ELSE rain(i) = zrfl(i) END IF END DO RETURN END SUBROUTINE fisrtilp_tr