SUBROUTINE conlmd (dtime, paprs, pplay, t, q, conv_q, s d_t, d_q, rain, snow, ibas, itop) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 c Objet: Schema de convection utilis'e dans le modele du LMD c Ajustement humide (Manabe) + Ajustement convectif (Kuo) c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOETHF.h" c c Arguments: c REAL dtime ! pas d'integration (s) REAL paprs(klon,klev+1) ! pression inter-couche (Pa) REAL pplay(klon,klev) ! pression au milieu de couche (Pa) REAL t(klon,klev) ! temperature (K) REAL q(klon,klev) ! humidite specifique (kg/kg) REAL conv_q(klon,klev) ! taux de convergence humidite (g/g/s) c REAL d_t(klon,klev) ! incrementation temperature REAL d_q(klon,klev) ! incrementation humidite REAL rain(klon) ! pluies (mm/s) REAL snow(klon) ! neige (mm/s) INTEGER ibas(klon) ! niveau du bas INTEGER itop(klon) ! niveau du haut c LOGICAL usekuo ! utiliser convection profonde (schema Kuo) PARAMETER (usekuo=.TRUE.) c REAL d_t_bis(klon,klev) REAL d_q_bis(klon,klev) REAL rain_bis(klon) REAL snow_bis(klon) INTEGER ibas_bis(klon) INTEGER itop_bis(klon) REAL d_ql(klon,klev), d_ql_bis(klon,klev) REAL rneb(klon,klev), rneb_bis(klon,klev) c INTEGER i, k REAL zlvdcp, zlsdcp, zdelta, zz, za, zb c ccc CALL fiajh ! ancienne version de Convection Manabe CALL conman ! nouvelle version de Convection Manabe e (dtime, paprs, pplay, t, q, s d_t, d_q, d_ql, rneb, s rain, snow, ibas, itop) c IF (usekuo) THEN ccc CALL fiajc ! ancienne version de Convection Kuo CALL conkuo ! nouvelle version de Convection Kuo e (dtime, paprs, pplay, t, q, conv_q, s d_t_bis, d_q_bis, d_ql_bis, rneb_bis, s rain_bis, snow_bis, ibas_bis, itop_bis) DO k = 1, klev DO i = 1, klon d_t(i,k) = d_t(i,k) + d_t_bis(i,k) d_q(i,k) = d_q(i,k) + d_q_bis(i,k) d_ql(i,k) = d_ql(i,k) + d_ql_bis(i,k) ENDDO ENDDO DO i = 1, klon rain(i) = rain(i) + rain_bis(i) snow(i) = snow(i) + snow_bis(i) ibas(i) = MIN(ibas(i),ibas_bis(i)) itop(i) = MAX(itop(i),itop_bis(i)) ENDDO ENDIF c c L'eau liquide convective est dispersee dans l'air: c DO k = 1, klev DO i = 1, klon zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q(i,k)) zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q(i,k)) zdelta = MAX(0.,SIGN(1.,RTT-t(i,k))) zz = d_ql(i,k) ! re-evap. de l'eau liquide zb = MAX(0.0,zz) za = - MAX(0.0,zz) * (zlvdcp*(1.-zdelta)+zlsdcp*zdelta) d_t(i,k) = d_t(i,k) + za d_q(i,k) = d_q(i,k) + zb ENDDO ENDDO c RETURN END SUBROUTINE conman (dtime, paprs, pplay, t, q, s d_t, d_q, d_ql, rneb, s rain, snow, ibas, itop) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19970324 c Objet: ajustement humide convectif avec la possibilite de faire c l'ajustement sur une fraction de la maille. c Methode: On impose une distribution uniforme pour la vapeur d'eau c au sein d'une maille. On applique la procedure d'ajustement c successivement a la totalite, 75%, 50%, 25% et 5% de la maille c jusqu'a ce que l'ajustement a lieu. J'espere que ceci augmente c les activites convectives et corrige le biais "trop froid et sec" c du modele. c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c REAL dtime ! pas d'integration (s) REAL t(klon,klev) ! temperature (K) REAL q(klon,klev) ! humidite specifique (kg/kg) REAL paprs(klon,klev+1) ! pression inter-couche (Pa) REAL pplay(klon,klev) ! pression au milieu de couche (Pa) c REAL d_t(klon,klev) ! incrementation temperature REAL d_q(klon,klev) ! incrementation humidite REAL d_ql(klon,klev) ! incrementation eau liquide REAL rneb(klon,klev) ! nebulosite REAL rain(klon) ! pluies (mm/s) REAL snow(klon) ! neige (mm/s) INTEGER ibas(klon) ! niveau du bas INTEGER itop(klon) ! niveau du haut c LOGICAL afaire(klon) ! .TRUE. implique l'ajustement LOGICAL accompli(klon) ! .TRUE. si l'ajustement est effectif c INTEGER nb ! nombre de sous-fractions a considere PARAMETER (nb=1) ccc PARAMETER (nb=3) c REAL ratqs ! largeur de la distribution pour vapeur d'eau PARAMETER (ratqs=0.05) c REAL w_q(klon,klev) REAL w_d_t(klon,klev), w_d_q(klon,klev), w_d_ql(klon,klev) REAL w_rneb(klon,klev) REAL w_rain(klon), w_snow(klon) INTEGER w_ibas(klon), w_itop(klon) REAL zq1, zq2 INTEGER i, k, n c REAL t_coup PARAMETER (t_coup=234.0) REAL zdp1, zdp2 REAL zqs1, zqs2, zdqs1, zdqs2 REAL zgamdz REAL zflo ! flotabilite REAL zsat ! sur-saturation REAL zdelta, zcor, zcvm5 LOGICAL imprim c INTEGER ncpt SAVE ncpt REAL frac(nb) ! valeur de la maille fractionnelle SAVE frac INTEGER opt_cld(nb) ! option pour le modele nuageux SAVE opt_cld LOGICAL appel1er SAVE appel1er c c Fonctions thermodynamiques: c #include "YOETHF.h" #include "FCTTRE.h" c DATA frac / 1.0 / DATA opt_cld / 4 / ccc DATA frac / 1.0, 0.50, 0.25/ ccc DATA opt_cld / 4, 4, 4/ c DATA appel1er /.TRUE./ DATA ncpt /0/ c IF (appel1er) THEN PRINT*, 'conman, nb:', nb PRINT*, 'conman, frac:', frac PRINT*, 'conman, opt_cld:', opt_cld appel1er = .FALSE. ENDIF c c Initialiser les sorties a zero: c 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 ENDDO ENDDO DO i = 1, klon ibas(i) = klev itop(i) = 1 rain(i) = 0.0 snow(i) = 0.0 ENDDO c c S'il n'y a pas d'instabilite conditionnelle, c pas la penne de se fatiguer: c DO i = 1, klon afaire(i) = .FALSE. ENDDO DO k = 1, klev-1 DO i = 1, klon IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-t(i,k))) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k)) zqs1= R2ES*FOEEW(t(i,k),zdelta)/pplay(i,k) zqs1=MIN(0.5,zqs1) zcor=1./(1.-RETV*zqs1) zqs1=zqs1*zcor zdqs1 =FOEDE(t(i,k),zdelta,zcvm5,zqs1,zcor) c zdelta=MAX(0.,SIGN(1.,RTT-t(i,k+1))) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k+1)) zqs2= R2ES*FOEEW(t(i,k+1),zdelta)/pplay(i,k+1) zqs2=MIN(0.5,zqs2) zcor=1./(1.-RETV*zqs2) zqs2=zqs2*zcor zdqs2 =FOEDE(t(i,k+1),zdelta,zcvm5,zqs2,zcor) ELSE IF (t(i,k) .LT. t_coup) THEN zqs1= qsats(t(i,k)) / pplay(i,k) zdqs1= dqsats(t(i,k),zqs1) c zqs2= qsats(t(i,k+1)) / pplay(i,k+1) zdqs2= dqsats(t(i,k+1),zqs2) ELSE zqs1= qsatl(t(i,k)) / pplay(i,k) zdqs1= dqsatl(t(i,k),zqs1) c zqs2= qsatl(t(i,k+1)) / pplay(i,k+1) zdqs2= dqsatl(t(i,k+1),zqs2) ENDIF ENDIF zdp1 = paprs(i,k) - paprs(i,k+1) zdp2 = paprs(i,k+1) - paprs(i,k+2) zgamdz = - (pplay(i,k)-pplay(i,k+1))/paprs(i,k+1)/RCPD . *( RD*(t(i,k)*zdp1+t(i,k+1)*zdp2)/(zdp1+zdp2) . +RLVTT*(zqs1*zdp1+zqs2*zdp2)/(zdp1+zdp2) . ) / (1.0+(zdqs1*zdp1+zdqs2*zdp2)/(zdp1+zdp2) ) zflo = t(i,k) + zgamdz - t(i,k+1) zsat = (q(i,k)-zqs1)*zdp1 + (q(i,k+1)-zqs2)*zdp2 IF (zflo.GT.0.0) afaire(i) = .TRUE. c erreur IF (zflo.GT.0.0 .AND. zsat.GT.0.0) afaire(i) = .TRUE. ENDDO ENDDO c imprim = MOD(ncpt,48).EQ.0 DO 99999 n = 1, nb c DO k = 1, klev DO i = 1, klon IF (afaire(i)) THEN zq1 = q(i,k) * (1.0-ratqs) zq2 = q(i,k) * (1.0+ratqs) w_q(i,k) = zq2 - frac(n)/2.0 * (zq2-zq1) ENDIF ENDDO ENDDO c CALL conmanv (dtime, paprs, pplay, t, w_q, e afaire, opt_cld(n), s w_d_t, w_d_q, w_d_ql, w_rneb, s w_rain, w_snow, w_ibas, w_itop,accompli,imprim) DO k = 1, klev DO i = 1, klon IF (afaire(i) .AND. accompli(i)) THEN d_t(i,k) = w_d_t(i,k) * frac(n) d_q(i,k) = w_d_q(i,k) * frac(n) d_ql(i,k) = w_d_ql(i,k) * frac(n) IF (NINT(w_rneb(i,k)).EQ.1) rneb(i,k) = frac(n) ENDIF ENDDO ENDDO DO i = 1, klon IF (afaire(i) .AND. accompli(i)) THEN rain(i) = w_rain(i) * frac(n) snow(i) = w_snow(i) * frac(n) ibas(i) = MIN(ibas(i),w_ibas(i)) itop(i) = MAX(itop(i),w_itop(i)) ENDIF ENDDO DO i = 1, klon IF(afaire(i) .AND. accompli(i)) afaire(i) = .FALSE. ENDDO c 99999 CONTINUE c ncpt = ncpt + 1 c RETURN END SUBROUTINE conmanv (dtime, paprs, pplay, t, q, e afaire, opt_cld, s d_t, d_q, d_ql, rneb, s rain, snow, ibas, itop,accompli,imprim) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 c Objet: ajustement humide (convection proposee par Manabe). c Pour une colonne verticale, il peut avoir plusieurs blocs c necessitant l'ajustement. ibas est le bas du plus bas bloc c et itop est le haut du plus haut bloc c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c c Arguments: c REAL dtime ! pas d'integration (s) REAL t(klon,klev) ! temperature (K) REAL q(klon,klev) ! humidite specifique (kg/kg) REAL paprs(klon,klev+1) ! pression inter-couche (Pa) REAL pplay(klon,klev) ! pression au milieu de couche (Pa) INTEGER opt_cld ! comment traiter l'eau liquide LOGICAL afaire(klon) ! .TRUE. si le point est a faire (Input) LOGICAL imprim ! .T. pour imprimer quelques diagnostiques c REAL d_t(klon,klev) ! incrementation temperature REAL d_q(klon,klev) ! incrementation humidite REAL d_ql(klon,klev) ! incrementation eau liquide REAL rneb(klon,klev) ! nebulosite REAL rain(klon) ! pluies (mm/s) REAL snow(klon) ! neige (mm/s) INTEGER ibas(klon) ! niveau du bas INTEGER itop(klon) ! niveau du haut LOGICAL accompli(klon) ! .TRUE. si l'ajustement a eu lieu (Output) c c Quelques options: c LOGICAL new_top ! re-calculer sommet quand re-ajustement est fait PARAMETER (new_top=.FALSE.) LOGICAL evap_prec ! evaporation de pluie au-dessous de convection PARAMETER (evap_prec=.TRUE.) REAL coef_eva PARAMETER (coef_eva=1.0E-05) REAL t_coup PARAMETER (t_coup=234.0) REAL seuil_vap PARAMETER (seuil_vap=1.0E-10) LOGICAL old_tau ! implique precip nulle, si vrai. PARAMETER (old_tau=.FALSE.) REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande PARAMETER (dpmin=0.15, tomax=0.97) REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible PARAMETER (dpmax=0.30, tomin=0.05) REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to PARAMETER (deep_sig=0.50, deep_to=0.05) LOGICAL exigent ! implique un calcul supplementaire pour Qs PARAMETER (exigent=.FALSE.) c INTEGER kbase PARAMETER (kbase=0) c c Variables locales: c INTEGER nexpo INTEGER i, k, k1min, k1max, k2min, k2max, is REAL zgamdz(klon,klev-1) REAL zt(klon,klev), zq(klon,klev) REAL zqs(klon,klev), zdqs(klon,klev) REAL zqmqsdp(klon,klev) REAL ztnew(klon,klev), zqnew(klon,klev) REAL zcond(klon), zvapo(klon), zrapp(klon) REAL zrfl(klon), zrfln, zqev, zqevt REAL zsat(klon) ! sur-saturation REAL zflo(klon) ! flotabilite REAL za(klon), zb(klon), zc(klon) INTEGER k1(klon), k2(klon) REAL zdelta, zcor, zcvm5 REAL delp(klon,klev) LOGICAL possible(klon), todo(klon), etendre(klon) LOGICAL aller(klon), todobis(klon) REAL zalfa INTEGER nbtodo, nbdone c c Fonctions thermodynamiques: c #include "YOETHF.h" #include "FCTTRE.h" c DO k = 1, klev DO i = 1, klon delp(i,k) = paprs(i,k) - paprs(i,k+1) ENDDO ENDDO c c Initialiser les sorties a zero c 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 ENDDO ENDDO DO i = 1, klon ibas(i) = klev itop(i) = 1 rain(i) = 0.0 snow(i) = 0.0 accompli(i) = .FALSE. ENDDO c c Preparations c DO k = 1, klev DO i = 1, klon IF (afaire(i)) THEN zt(i,k) = t(i,k) zq(i,k) = q(i,k) c c Calculer Qs et L/Cp*dQs/dT c IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-zt(i,k))) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i,k)) zqs(i,k)= R2ES*FOEEW(zt(i,k),zdelta)/pplay(i,k) zqs(i,k)=MIN(0.5,zqs(i,k)) zcor=1./(1.-RETV*zqs(i,k)) zqs(i,k)=zqs(i,k)*zcor zdqs(i,k) =FOEDE(zt(i,k),zdelta,zcvm5,zqs(i,k),zcor) ELSE IF (zt(i,k) .LT. t_coup) THEN zqs(i,k)= qsats(zt(i,k)) / pplay(i,k) zdqs(i,k)= dqsats(zt(i,k),zqs(i,k)) ELSE zqs(i,k)= qsatl(zt(i,k)) / pplay(i,k) zdqs(i,k)= dqsatl(zt(i,k),zqs(i,k)) ENDIF ENDIF c c Calculer (q-qs)*dp zqmqsdp(i,k) = (zq(i,k)-zqs(i,k)) * delp(i,k) ENDIF ENDDO ENDDO c c-----zgama is the moist convective lapse rate (-dT/dz). c-----zgamdz(*,k) est la difference minimale autorisee des temperatures c-----entre deux couches (k et k+1), c.a.d. si T(k+1)-T(k) est inferieur c-----a zgamdz(*,k), alors ces 2 couches sont instables conditionnellement c DO k = 1, klev-1 DO i = 1, klon IF (afaire(i)) THEN zgamdz(i,k) = - (pplay(i,k)-pplay(i,k+1))/paprs(i,k+1)/RCPD . *( RD*(zt(i,k)*delp(i,k)+zt(i,k+1)*delp(i,k+1)) . /(delp(i,k)+delp(i,k+1)) . +RLVTT*(zqs(i,k)*delp(i,k)+zqs(i,k+1)*delp(i,k+1)) . /(delp(i,k)+delp(i,k+1)) . ) / (1.0+(zdqs(i,k)*delp(i,k)+zdqs(i,k+1)*delp(i,k+1)) . /(delp(i,k)+delp(i,k+1)) ) ENDIF ENDDO ENDDO c c On cherche la presence simultanee d'instabilite conditionnelle c et de sur-saturation. Sinon, pas la penne de se fatiguer: c DO i = 1, klon possible(i) = .FALSE. ENDDO DO k = 2, klev DO i = 1, klon IF (afaire(i)) THEN zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k) zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1) IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) possible(i) = .TRUE. ENDIF ENDDO ENDDO c DO i = 1, klon IF (possible(i)) THEN k1(i) = kbase k2(i) = k1(i) + 1 ENDIF ENDDO c 810 CONTINUE ! chercher le bas de la colonne a ajuster c k2min = klev DO i = 1, klon todo(i) = .FALSE. aller(i) = .TRUE. IF (possible(i)) k2min = MIN(k2min,k2(i)) ENDDO IF (k2min.EQ.klev) GOTO 860 DO k = k2min, klev-1 DO i = 1, klon IF (possible(i) .AND. k.GE.k2(i) .AND. aller(i)) THEN zflo(i) = zt(i,k) + zgamdz(i,k) - zt(i,k+1) zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k+1) IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN k1(i) = k k2(i) = k+1 todo(i) = .TRUE. aller(i) = .FALSE. ENDIF ENDIF ENDDO ENDDO DO i = 1, klon IF (possible(i).AND.aller(i)) THEN todo(i) = .FALSE. k1(i) = klev k2(i) = klev ENDIF ENDDO c CCC DO i = 1, klon CCC IF (possible(i)) THEN CCC 811 k2(i) = k2(i) + 1 CCC IF (k2(i) .GT. klev) THEN CCC todo(i) = .FALSE. CCC GOTO 812 CCC ENDIF CCC k = k2(i) CCC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k) CCC zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1) CCC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 811 CCC k1(i) = k2(i) - 1 CCC todo(i) = .TRUE. CCC ENDIF CCC 812 CONTINUE CCC ENDDO c 820 CONTINUE ! chercher le haut de la colonne c k2min = klev DO i = 1, klon aller(i) = .TRUE. IF (todo(i)) k2min = MIN(k2min,k2(i)) ENDDO IF (k2min.LT.klev) THEN DO k = k2min, klev DO i = 1, klon IF (todo(i) .AND. k.GT.k2(i) .AND. aller(i)) THEN zsat(i) = zsat(i) + zqmqsdp(i,k) zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k) IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN aller(i) = .FALSE. ELSE k2(i) = k ENDIF ENDIF ENDDO ENDDO c error is = 0 c error DO i = 1, klon c error IF(todo(i).AND.aller(i)) THEN c error is = is + 1 c error todo(i) = .FALSE. c error k2(i) = klev c error ENDIF c error ENDDO c error IF (is.GT.0) THEN c error PRINT*, "Bizard. je pourrais continuer mais j arrete" c error CALL abort c error ENDIF ENDIF c CCC DO i = 1, klon CCC IF (todo(i)) THEN CCC 821 CONTINUE CCC IF (k2(i) .EQ. klev) GOTO 822 CCC k = k2(i) + 1 CCC zsat(i) = zsat(i) + zqmqsdp(i,k) CCC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k) CCC IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) GOTO 822 CCC k2(i) = k CCC GOTO 821 CCC ENDIF CCC 822 CONTINUE CCC ENDDO c 830 CONTINUE ! faire l'ajustement en sachant k1 et k2 c is = 0 DO i = 1, klon IF (todo(i)) THEN IF (k2(i).LE.k1(i)) is = is + 1 ENDIF ENDDO IF (is.GT.0) THEN PRINT*, "Impossible: k1 trop grand ou k2 trop petit" PRINT*, "is=", is CALL abort ENDIF c k1min = klev k1max = 1 k2max = 1 DO i = 1, klon IF (todo(i)) THEN k1min = MIN(k1min,k1(i)) k1max = MAX(k1max,k1(i)) k2max = MAX(k2max,k2(i)) ENDIF ENDDO c DO i = 1, klon IF (todo(i)) THEN k = k1(i) za(i) = 0. zb(i) = ( RCPD*(1.+zdqs(i,k))*(zt(i,k)-za(i)) . -RLVTT*(zqs(i,k)-zq(i,k)) )*delp(i,k) zc(i) = delp(i,k) * RCPD*(1.+zdqs(i,k)) ENDIF ENDDO c DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i)) THEN za(i) = za(i) + zgamdz(i,k-1) zb(i) = zb(i)+(RCPD*(1.+zdqs(i,k))*(zt(i,k)-za(i)) . -RLVTT*(zqs(i,k)-zq(i,k)) ) * delp(i,k) zc(i) = zc(i) + delp(i,k)*RCPD*(1.+zdqs(i,k)) ENDIF ENDDO ENDDO c DO i = 1, klon IF (todo(i)) THEN k = k1(i) ztnew(i,k) = zb(i)/zc(i) zqnew(i,k) = zqs(i,k) + (ztnew(i,k)-zt(i,k)) . *RCPD/RLVTT*zdqs(i,k) ENDIF ENDDO c DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i)) THEN ztnew(i,k) = ztnew(i,k-1) + zgamdz(i,k-1) zqnew(i,k) = zqs(i,k) + (ztnew(i,k)-zt(i,k)) . *RCPD/RLVTT*zdqs(i,k) ENDIF ENDDO ENDDO c c Quantite de condensation produite pendant l'ajustement: c DO i = 1, klon zcond(i) = 0.0 ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN rneb(i,k) = 1.0 zcond(i) = zcond(i) + (zq(i,k)-zqnew(i,k)) *delp(i,k)/RG ENDIF ENDDO ENDDO c c Si condensation negative, effort completement perdu: c DO i = 1, klon IF (todo(i).AND.zcond(i).LE.0.) todo(i) = .FALSE. ENDDO c c L'ajustement a ete accompli, meme les calculs accessoires c ne sont pas encore faits: c DO i = 1, klon IF (todo(i)) accompli(i) = .TRUE. ENDDO c c===== c Une fois que la condensation a lieu, on doit construire un c "modele nuageux" pour partager la condensation entre l'eau c liquide nuageuse et la precipitation (leur rapport toliq c est calcule selon l'epaisseur nuageuse). Je suppose que c toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin, c et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation c lineaire entre dpmin et dpmax). c===== DO i = 1, klon IF (todo(i)) THEN toliq(i) = tomax-((paprs(i,k1(i))-paprs(i,k2(i)+1)) . /paprs(i,1)-dpmin) . *(tomax-tomin)/(dpmax-dpmin) toliq(i) = MAX(tomin,MIN(tomax,toliq(i))) IF (pplay(i,k2(i))/paprs(i,1) .LE. deep_sig) toliq(i) = deep_to IF (old_tau) toliq(i) = 1.0 ENDIF ENDDO c===== c On doit aussi determiner la distribution verticale de c l'eau nuageuse. Plusieurs options sont proposees: c c (0) La condensation precipite integralement (toliq ne sera c pas utilise). c (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle c a la vapeur d'eau locale. c (2) Elle est distribuee entre k1 et k2 avec une valeur constante. c (3) Elle est seulement distribuee aux couches ou la vapeur d'eau c est effectivement diminuee pendant le processus d'ajustement. c (4) Elle est en fonction (lineaire ou exponentielle) de la c distance (epaisseur en pression) avec le niveau k1 (la couche c k1 n'aura donc pas d'eau liquide). c===== c IF (opt_cld.EQ.0) THEN c DO i = 1, klon IF (todo(i)) zrfl(i) = zcond(i) / dtime ENDDO c ELSE IF (opt_cld.EQ.1) THEN c DO i = 1, klon IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) . zvapo(i) = zvapo(i) + zqnew(i,k)*delp(i,k)/RG ENDDO ENDDO DO i = 1, klon IF (todo(i)) THEN zrapp(i) = toliq(i) * zcond(i) / zvapo(i) zrapp(i) = MAX(0.,MIN(1.,zrapp(i))) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDIF ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN d_ql(i,k) = d_ql(i,k) + zrapp(i) * zqnew(i,k) ENDIF ENDDO ENDDO c ELSE IF (opt_cld.EQ.2) THEN c DO i = 1, klon IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) . zvapo(i) = zvapo(i) + delp(i,k)/RG ENDDO ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i) ENDIF ENDDO ENDDO DO i = 1, klon IF (todo(i)) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDDO c ELSE IF (opt_cld.EQ.3) THEN c DO i = 1, klon IF (todo(i)) zvapo(i) = 0.0 ! quantite de l'eau strictement condensee ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) . zvapo(i) = zvapo(i) + MAX(0.0,zq(i,k)-zqnew(i,k)) . * delp(i,k)/RG ENDDO ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i) .AND. . zvapo(i).GT.0.0) . d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i) . * MAX(0.0,zq(i,k)-zqnew(i,k)) ENDDO ENDDO DO i = 1, klon IF (todo(i)) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDDO c ELSE IF (opt_cld.EQ.4) THEN c nexpo = 3 ccc nexpo = 1 ! distribution lineaire c DO i = 1, klon IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse ENDDO ! (avec ponderation) DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i)) . zvapo(i) = zvapo(i) + delp(i,k) / RG . * (pplay(i,k1(i))-pplay(i,k))**nexpo ENDDO ENDDO DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.(k1(i)+1) .AND. k.LE.k2(i)) . d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i) . * (pplay(i,k1(i))-pplay(i,k))**nexpo ENDDO ENDDO DO i = 1, klon IF (todo(i)) zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDDO c ELSE ! valeur non-prevue pour opt_cld c PRINT*, "opt_cld est faux:", opt_cld CALL abort c ENDIF ! fin de opt_cld c c L'eau precipitante peut etre evaporee: c zalfa = 0.05 IF (evap_prec .AND. (k1max.GE.2)) THEN DO k = k1max-1, 1, -1 DO i = 1, klon IF (todo(i) .AND. k.LT.k1(i) .AND. zrfl(i).GT.0.0) THEN zqev = MAX (0.0, (zqs(i,k)-zq(i,k))*zalfa ) zqevt = coef_eva * (1.0-zq(i,k)/zqs(i,k))*SQRT(zrfl(i)) . * delp(i,k)/pplay(i,k)*zt(i,k)*RD/RG zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/delp(i,k) zqev = MIN (zqev, zqevt) zrfln = zrfl(i) - zqev*(delp(i,k))/RG/dtime zq(i,k) = zq(i,k) - (zrfln-zrfl(i)) . * (RG/(delp(i,k)))*dtime zt(i,k) = zt(i,k) + (zrfln-zrfl(i)) . * (RG/(delp(i,k)))*dtime . * RLVTT/RCPD/(1.0+RVTMP2*zq(i,k)) zrfl(i) = zrfln ENDIF ENDDO ENDDO ENDIF c c La temperature de la premiere couche determine la pluie ou la neige: c DO i = 1, klon IF (todo(i)) THEN IF (zt(i,1) .GT. RTT) THEN rain(i) = rain(i) + zrfl(i) ELSE snow(i) = snow(i) + zrfl(i) ENDIF ENDIF ENDDO c c Mise a jour de la temperature et de l'humidite c DO k = k1min, k2max DO i = 1, klon IF (todo(i) .AND. k.GE.k1(i) .AND. k.LE.k2(i)) THEN zt(i,k) = ztnew(i,k) zq(i,k) = zqnew(i,k) ENDIF ENDDO ENDDO c c Re-calculer certaines variables pour etendre et re-ajuster la colonne c IF (exigent) THEN DO k = 1, klev DO i = 1, klon IF (todo(i)) THEN IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-zt(i,k))) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i,k)) zqs(i,k)= R2ES*FOEEW(zt(i,k),zdelta)/pplay(i,k) zqs(i,k)=MIN(0.5,zqs(i,k)) zcor=1./(1.-RETV*zqs(i,k)) zqs(i,k)=zqs(i,k)*zcor zdqs(i,k) =FOEDE(zt(i,k),zdelta,zcvm5,zqs(i,k),zcor) ELSE IF (zt(i,k) .LT. t_coup) THEN zqs(i,k)= qsats(zt(i,k)) / pplay(i,k) zdqs(i,k)= dqsats(zt(i,k),zqs(i,k)) ELSE zqs(i,k)= qsatl(zt(i,k)) / pplay(i,k) zdqs(i,k)= dqsatl(zt(i,k),zqs(i,k)) ENDIF ENDIF ENDIF ENDDO ENDDO ENDIF c IF (exigent) THEN DO k = 1, klev-1 DO i = 1, klon IF (todo(i)) THEN zgamdz(i,k) = - (pplay(i,k)-pplay(i,k+1))/paprs(i,k+1)/RCPD . *( RD*(zt(i,k)*delp(i,k)+zt(i,k+1)*delp(i,k+1)) . /(delp(i,k)+delp(i,k+1)) . +RLVTT*(zqs(i,k)*delp(i,k)+zqs(i,k+1)*delp(i,k+1)) . /(delp(i,k)+delp(i,k+1)) . ) / (1.0+(zdqs(i,k)*delp(i,k)+zdqs(i,k+1)*delp(i,k+1)) . /(delp(i,k)+delp(i,k+1)) ) ENDIF ENDDO ENDDO ENDIF c c Puisque l'humidite a ete modifiee, on re-fait (q-qs)*dp c DO k = 1, klev DO i = 1, klon IF (todo(i)) THEN zqmqsdp(i,k) = (zq(i,k)-zqs(i,k))*delp(i,k) ENDIF ENDDO ENDDO c c Verifier si l'on peut etendre le bas de la colonne c DO i = 1, klon etendre(i) = .FALSE. ENDDO c k1max = 1 DO i = 1, klon IF (todo(i) .AND. k1(i).GT.(kbase+1)) THEN k = k1(i) zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k) zsat(i) = zqmqsdp(i,k) + zqmqsdp(i,k-1) csc voici l'ancienne ligne: csc IF (zflo(i).LE.0.0 .OR. zsat(i).LE.0.0) THEN csc sylvain: il faut RESPECTER les 2 criteres: IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN etendre(i) = .TRUE. k1(i) = k1(i) - 1 k1max = MAX(k1max,k1(i)) aller(i) = .TRUE. ENDIF ENDIF ENDDO c IF (k1max.GT.(kbase+1)) THEN DO k = k1max, kbase+1, -1 DO i = 1, klon IF (etendre(i) .AND. k.LT.k1(i) .AND. aller(i)) THEN zsat(i) = zsat(i) + zqmqsdp(i,k) zflo(i) = zt(i,k) + zgamdz(i,k) - zt(i,k+1) IF (zsat(i).LE.0.0 .OR. zflo(i).LE.0.0) THEN aller(i) = .FALSE. ELSE k1(i) = k ENDIF ENDIF ENDDO ENDDO DO i = 1, klon IF (etendre(i).AND.aller(i)) THEN k1(i) = 1 ENDIF ENDDO ENDIF c CCC DO i = 1, klon CCC IF (etendre(i)) THEN CCC 840 k = k1(i) CCC IF (k.GT.1) THEN CCC zsat(i) = zsat(i) + zqmqsdp(i,k-1) CCC zflo(i) = zt(i,k-1) + zgamdz(i,k-1) - zt(i,k) CCC IF (zflo(i).GT.0.0 .AND. zsat(i).GT.0.0) THEN CCC k1(i) = k - 1 CCC GOTO 840 CCC ENDIF CCC ENDIF CCC ENDIF CCC ENDDO c DO i = 1, klon todobis(i) = todo(i) todo(i) = .FALSE. ENDDO is = 0 DO i = 1, klon IF (etendre(i)) THEN todo(i) = .TRUE. is = is + 1 ENDIF ENDDO IF (is.GT.0) THEN IF (new_top) THEN GOTO 820 ! chercher de nouveau le sommet k2 ELSE GOTO 830 ! supposer que le sommet est celui deja trouve ENDIF ENDIF c DO i = 1, klon possible(i) = .FALSE. ENDDO is = 0 DO i = 1, klon IF (todobis(i) .AND. k2(i).LT.klev) THEN is = is + 1 possible(i) = .TRUE. ENDIF ENDDO IF (is.GT.0) GOTO 810 !on cherche en haut d'autres blocks c a ajuster a partir du sommet de la colonne precedente c 860 CONTINUE ! Calculer les tendances et diagnostiques ccc print*, "Apres 860" c DO k = 1, klev DO i = 1, klon IF (accompli(i)) THEN d_t(i,k) = zt(i,k) - t(i,k) zq(i,k) = MAX(zq(i,k),seuil_vap) d_q(i,k) = zq(i,k) - q(i,k) ENDIF ENDDO ENDDO c DO 888 i = 1, klon IF (accompli(i)) THEN DO k = 1, klev IF (rneb(i,k).GT.0.0) THEN ibas(i) = k GOTO 807 ENDIF ENDDO 807 CONTINUE DO k = klev, 1, -1 IF (rneb(i,k).GT.0.0) THEN itop(i) = k GOTO 808 ENDIF ENDDO 808 CONTINUE ENDIF 888 CONTINUE c IF (imprim) THEN nbtodo = 0 nbdone = 0 DO i = 1, klon IF (afaire(i)) nbtodo = nbtodo + 1 IF (accompli(i)) nbdone = nbdone + 1 ENDDO PRINT*, "nbTodo, nbDone=", nbtodo, nbdone ENDIF c RETURN END SUBROUTINE conkuo(dtime, paprs, pplay, t, q, conv_q, s d_t, d_q, d_ql, rneb, s rain, snow, ibas, itop) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818 c Objet: Schema de convection de type Kuo (1965). c Cette version du code peut calculer le niveau de depart c N.B. version vectorielle (le 6 oct. 1997) c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c c Arguments: c REAL dtime ! intervalle du temps (s) REAL paprs(klon,klev+1) ! pression a inter-couche (Pa) REAL pplay(klon,klev) ! pression au milieu de couche (Pa) REAL t(klon,klev) ! temperature (K) REAL q(klon,klev) ! humidite specifique REAL conv_q(klon,klev) ! taux de convergence humidite (g/g/s) c REAL d_t(klon,klev) ! incrementation temperature REAL d_q(klon,klev) ! incrementation humidite REAL d_ql(klon,klev) ! incrementation eau liquide REAL rneb(klon,klev) ! nebulosite REAL rain(klon) ! pluies (mm/s) REAL snow(klon) ! neige (mm/s) INTEGER itop(klon) ! niveau du sommet INTEGER ibas(klon) ! niveau du bas c LOGICAL ldcum(klon) ! convection existe LOGICAL todo(klon) c c Quelsques options: c LOGICAL calcfcl ! calculer le niveau de convection libre PARAMETER (calcfcl=.TRUE.) INTEGER ldepar ! niveau fixe de convection libre PARAMETER (ldepar=4) INTEGER opt_cld ! comment traiter l'eau liquide PARAMETER (opt_cld=4) ! valeur possible: 0, 1, 2, 3 ou 4 LOGICAL evap_prec ! evaporation de pluie au-dessous de convection PARAMETER (evap_prec=.TRUE.) REAL coef_eva PARAMETER (coef_eva=1.0E-05) LOGICAL new_deh ! nouvelle facon de calculer dH PARAMETER (new_deh=.FALSE.) REAL t_coup PARAMETER (t_coup=234.0) LOGICAL old_tau ! implique precipitation nulle PARAMETER (old_tau=.FALSE.) REAL toliq(klon) ! rapport entre l'eau nuageuse et l'eau precipitante REAL dpmin, tomax !Epaisseur faible, rapport eau liquide plus grande PARAMETER (dpmin=0.15, tomax=0.97) REAL dpmax, tomin !Epaisseur grande, rapport eau liquide plus faible PARAMETER (dpmax=0.30, tomin=0.05) REAL deep_sig, deep_to ! au dela de deep_sig, utiliser deep_to PARAMETER (deep_sig=0.50, deep_to=0.05) c c Variables locales: c INTEGER nexpo LOGICAL nuage(klon) INTEGER i, k, kbmin, kbmax, khmax REAL ztotal(klon,klev), zdeh(klon,klev) REAL zgz(klon,klev) REAL zqs(klon,klev) REAL zdqs(klon,klev) REAL ztemp(klon,klev) REAL zpres(klon,klev) REAL zconv(klon) ! convergence d'humidite REAL zvirt(klon) ! convergence virtuelle d'humidite REAL zfrac(klon) ! fraction convective INTEGER kb(klon), kh(klon) c REAL zcond(klon), zvapo(klon), zrapp(klon) REAL zrfl(klon), zrfln, zqev, zqevt REAL zdelta, zcvm5, zcor REAL zvar c LOGICAL appel1er SAVE appel1er c c Fonctions thermodynamiques c #include "YOETHF.h" #include "FCTTRE.h" c DATA appel1er /.TRUE./ c IF (appel1er) THEN PRINT*, 'conkuo, calcfcl:', calcfcl IF (.NOT.calcfcl) PRINT*, 'conkuo, ldepar:', ldepar PRINT*, 'conkuo, opt_cld:', opt_cld PRINT*, 'conkuo, evap_prec:', evap_prec PRINT*, 'conkuo, new_deh:', new_deh appel1er = .FALSE. ENDIF c c Initialiser les sorties a zero c DO k = 1, klev DO i = 1, klon d_q(i,k) = 0.0 d_t(i,k) = 0.0 d_ql(i,k) = 0.0 rneb(i,k) = 0.0 ENDDO ENDDO DO i = 1, klon rain(i) = 0.0 snow(i) = 0.0 ibas(i) = 0 itop(i) = 0 ENDDO c c Calculer la vapeur d'eau saturante Qs et sa derive L/Cp * dQs/dT c DO k = 1, klev DO i = 1, klon IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-t(i,k))) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k)) zqs(i,k)=R2ES*FOEEW(t(i,k),zdelta)/pplay(i,k) zqs(i,k)=MIN(0.5,zqs(i,k)) zcor=1./(1.-RETV*zqs(i,k)) zqs(i,k)=zqs(i,k)*zcor zdqs(i,k) =FOEDE(t(i,k),zdelta,zcvm5,zqs(i,k),zcor) ELSE IF (t(i,k).LT.t_coup) THEN zqs(i,k) = qsats(t(i,k))/pplay(i,k) zdqs(i,k) = dqsats(t(i,k),zqs(i,k)) ELSE zqs(i,k) = qsatl(t(i,k))/pplay(i,k) zdqs(i,k) = dqsatl(t(i,k),zqs(i,k)) ENDIF ENDIF ENDDO ENDDO c c Calculer gz (energie potentielle) c DO i = 1, klon zgz(i,1) = RD * t(i,1) / (0.5*(paprs(i,1)+pplay(i,1))) . * (paprs(i,1)-pplay(i,1)) ENDDO DO k = 2, klev DO i = 1, klon zgz(i,k) = zgz(i,k-1) . + RD * 0.5*(t(i,k-1)+t(i,k)) / paprs(i,k) . * (pplay(i,k-1)-pplay(i,k)) ENDDO ENDDO c c Calculer l'energie statique humide saturee (Cp*T + gz + L*Qs) c DO k = 1, klev DO i = 1, klon ztotal(i,k) = RCPD*t(i,k) + RLVTT*zqs(i,k) + zgz(i,k) ENDDO ENDDO c c Determiner le niveau de depart et calculer la difference de c l'energie statique humide saturee (ztotal) entre la couche c de depart et chaque couche au-dessus. c IF (calcfcl) THEN DO k = 1, klev DO i = 1, klon zpres(i,k) = pplay(i,k) ztemp(i,k) = t(i,k) ENDDO ENDDO CALL kuofcl(ztemp, q, zgz, zpres, ldcum, kb) DO i = 1, klon IF (ldcum(i)) THEN k = kb(i) IF (new_deh) THEN zdeh(i,k) = ztotal(i,k-1) - ztotal(i,k) ELSE zdeh(i,k) = RCPD * (t(i,k-1)-t(i,k)) . - RD *0.5*(t(i,k-1)+t(i,k))/paprs(i,k) . *(pplay(i,k-1)-pplay(i,k)) . + RLVTT*(zqs(i,k-1)-zqs(i,k)) ENDIF zdeh(i,k) = zdeh(i,k) * 0.5 ENDIF ENDDO DO k = 1, klev DO i = 1, klon IF (ldcum(i) .AND. k.GE.(kb(i)+1)) THEN IF (new_deh) THEN zdeh(i,k) = zdeh(i,k-1) + (ztotal(i,k-1)-ztotal(i,k)) ELSE zdeh(i,k) = zdeh(i,k-1) . + RCPD * (t(i,k-1)-t(i,k)) . - RD *0.5*(t(i,k-1)+t(i,k))/paprs(i,k) . *(pplay(i,k-1)-pplay(i,k)) . + RLVTT*(zqs(i,k-1)-zqs(i,k)) ENDIF ENDIF ENDDO ENDDO ELSE DO i = 1, klon k = ldepar kb(i) = ldepar ldcum(i) = .TRUE. IF (new_deh) THEN zdeh(i,k) = ztotal(i,k-1) - ztotal(i,k) ELSE zdeh(i,k) = RCPD * (t(i,k-1)-t(i,k)) . - RD *0.5*(t(i,k-1)+t(i,k))/paprs(i,k) . *(pplay(i,k-1)-pplay(i,k)) . + RLVTT*(zqs(i,k-1)-zqs(i,k)) ENDIF zdeh(i,k) = zdeh(i,k) * 0.5 ENDDO DO k = ldepar+1, klev DO i = 1, klon IF (new_deh) THEN zdeh(i,k) = zdeh(i,k-1) + (ztotal(i,k-1)-ztotal(i,k)) ELSE zdeh(i,k) = zdeh(i,k-1) . + RCPD * (t(i,k-1)-t(i,k)) . - RD *0.5*(t(i,k-1)+t(i,k))/paprs(i,k) . *(pplay(i,k-1)-pplay(i,k)) . + RLVTT*(zqs(i,k-1)-zqs(i,k)) ENDIF ENDDO ENDDO ENDIF c c-----Chercher le sommet du nuage c-----Calculer la convergence de l'humidite (en kg/m**2 a un facteur c-----psolpa/RG pres) du bas jusqu'au sommet du nuage. c-----Calculer la convergence virtuelle pour que toute la maille c-----deviennt nuageuse (du bas jusqu'au sommet du nuage) c DO i = 1, klon nuage(i) = .TRUE. zconv(i) = 0.0 zvirt(i) = 0.0 kh(i) = -999 ENDDO DO k = 1, klev DO i = 1, klon IF (k.GE.kb(i) .AND. ldcum(i)) THEN nuage(i) = nuage(i) .AND. zdeh(i,k).GT.0.0 IF (nuage(i)) THEN kh(i) = k zconv(i)=zconv(i)+conv_q(i,k)*dtime . *(paprs(i,k)-paprs(i,k+1)) zvirt(i)=zvirt(i)+(zdeh(i,k)/RLVTT+zqs(i,k)-q(i,k)) . *(paprs(i,k)-paprs(i,k+1)) ENDIF ENDIF ENDDO ENDDO c DO i = 1, klon todo(i) = ldcum(i) .AND. kh(i).GT.kb(i) .AND. zconv(i).GT.0.0 ENDDO c kbmin = klev kbmax = 0 khmax = 0 DO i = 1, klon IF (todo(i)) THEN kbmin = MIN(kbmin,kb(i)) kbmax = MAX(kbmax,kb(i)) khmax = MAX(khmax,kh(i)) ENDIF ENDDO c c-----Calculer la surface couverte par le nuage c DO i = 1, klon IF (todo(i)) THEN zfrac(i) = MAX(0.0,MIN(zconv(i)/zvirt(i), 1.0)) ENDIF ENDDO c c-----Calculs essentiels: c DO i = 1, klon IF (todo(i)) THEN zcond(i) = 0.0 ENDIF ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN zvar = zdeh(i,k)/(1.+zdqs(i,k)) d_t(i,k) = zvar * zfrac(i) / RCPD d_q(i,k) = (zvar*zdqs(i,k)/RLVTT+zqs(i,k)-q(i,k))*zfrac(i) . - conv_q(i,k)*dtime zcond(i) = zcond(i) - d_q(i,k) *(paprs(i,k)-paprs(i,k+1))/RG rneb(i,k) = zfrac(i) ENDIF ENDDO ENDDO c DO i = 1, klon IF (todo(i) .AND. zcond(i).LT.0.0) THEN PRINT*, 'WARNING: cond. negative (Kuo) ', . i,kb(i),kh(i), zcond(i) zcond(i) = 0.0 DO k = kb(i), kh(i) d_t(i,k) = 0.0 d_q(i,k) = 0.0 ENDDO todo(i) = .FALSE. ! effort totalement perdu ENDIF ENDDO c c===== c Une fois que la condensation a lieu, on doit construire un c "modele nuageux" pour partager la condensation entre l'eau c liquide nuageuse et la precipitation (leur rapport toliq c est calcule selon l'epaisseur nuageuse). Je suppose que c toliq=tomax quand l'epaisseur nuageuse est inferieure a dpmin, c et que toliq=tomin quand l'epaisseur depasse dpmax (interpolation c lineaire entre dpmin et dpmax). c===== DO i = 1, klon IF (todo(i)) THEN toliq(i) = tomax-((paprs(i,kb(i))-paprs(i,kh(i)+1)) . /paprs(i,1)-dpmin) . *(tomax-tomin)/(dpmax-dpmin) toliq(i) = MAX(tomin,MIN(tomax,toliq(i))) IF (pplay(i,kh(i))/paprs(i,1) .LE. deep_sig) toliq(i) = deep_to IF (old_tau) toliq(i) = 1.0 ENDIF ENDDO c===== c On doit aussi determiner la distribution verticale de c l'eau nuageuse. Plusieurs options sont proposees: c c (0) La condensation precipite integralement (toliq ne sera c pas utilise). c (1) L'eau liquide est distribuee entre k1 et k2 et proportionnelle c a la vapeur d'eau locale. c (2) Elle est distribuee entre k1 et k2 avec une valeur constante. c (3) Elle est seulement distribuee aux couches ou la vapeur d'eau c est effectivement diminuee pendant le processus d'ajustement. c (4) Elle est en fonction (lineaire ou exponentielle) de la c distance (epaisseur en pression) avec le niveau k1 (la couche c k1 n'aura donc pas d'eau liquide). c===== c IF (opt_cld.EQ.0) THEN c DO i = 1, klon IF (todo(i)) zrfl(i) = zcond(i) / dtime ENDDO c ELSE IF (opt_cld.EQ.1) THEN c DO i = 1, klon IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de vapeur d'eau ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN zvapo(i) = zvapo(i) + (q(i,k)+d_q(i,k)) . *(paprs(i,k)-paprs(i,k+1))/RG ENDIF ENDDO ENDDO DO i = 1, klon IF (todo(i)) THEN zrapp(i) = toliq(i) * zcond(i) / zvapo(i) zrapp(i) = MAX(0.,MIN(1.,zrapp(i))) ENDIF ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN d_ql(i,k) = zrapp(i) * (q(i,k)+d_q(i,k)) ENDIF ENDDO ENDDO DO i = 1, klon IF (todo(i)) THEN zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDIF ENDDO c ELSE IF (opt_cld.EQ.2) THEN c DO i = 1, klon IF (todo(i)) zvapo(i) = 0.0 ! quantite integrale de masse ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1))/RG ENDIF ENDDO ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN d_ql(i,k) = toliq(i) * zcond(i) / zvapo(i) ENDIF ENDDO ENDDO DO i = 1, klon IF (todo(i)) THEN zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDIF ENDDO c ELSE IF (opt_cld.EQ.3) THEN c DO i = 1, klon IF (todo(i)) THEN zvapo(i) = 0.0 ! quantite de l'eau strictement condensee ENDIF ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i)) THEN zvapo(i) = zvapo(i) + MAX(0.0,-d_q(i,k)) . * (paprs(i,k)-paprs(i,k+1))/RG ENDIF ENDDO ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.kb(i) .AND. k.LE.kh(i) .AND. . zvapo(i).GT.0.0) THEN d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i) . * MAX(0.0,-d_q(i,k)) ENDIF ENDDO ENDDO DO i = 1, klon IF (todo(i)) THEN zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDIF ENDDO c ELSE IF (opt_cld.EQ.4) THEN c nexpo = 3 ccc nexpo = 1 ! distribution lineaire c DO i = 1, klon IF (todo(i)) THEN zvapo(i) = 0.0 ! quantite integrale de masse (avec ponderation) ENDIF ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.(kb(i)+1) .AND. k.LE.kh(i)) THEN zvapo(i) = zvapo(i) + (paprs(i,k)-paprs(i,k+1)) / RG . * (pplay(i,kb(i))-pplay(i,k))**nexpo ENDIF ENDDO ENDDO DO k = kbmin, khmax DO i = 1, klon IF (todo(i) .AND. k.GE.(kb(i)+1) .AND. k.LE.kh(i)) THEN d_ql(i,k) = d_ql(i,k) + toliq(i) * zcond(i) / zvapo(i) . * (pplay(i,kb(i))-pplay(i,k))**nexpo ENDIF ENDDO ENDDO DO i = 1, klon IF (todo(i)) THEN zrfl(i) = (1.0-toliq(i)) * zcond(i) / dtime ENDIF ENDDO c ELSE ! valeur non-prevue pour opt_cld c PRINT*, "opt_cld est faux:", opt_cld CALL abort c ENDIF ! fin de opt_cld c c L'eau precipitante peut etre re-evaporee: c IF (evap_prec .AND. kbmax.GE.2) THEN DO k = kbmax, 1, -1 DO i = 1, klon IF (todo(i) .AND. k.LE.(kb(i)-1) .AND. zrfl(i).GT.0.0) THEN zqev = MAX (0.0, (zqs(i,k)-q(i,k))*zfrac(i) ) zqevt = coef_eva * (1.0-q(i,k)/zqs(i,k))*SQRT(zrfl(i)) . * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*t(i,k)*RD/RG zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) . * RG*dtime/(paprs(i,k)-paprs(i,k+1)) zqev = MIN (zqev, zqevt) zrfln = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) . /RG/dtime d_q(i,k) = - (zrfln-zrfl(i)) . * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime d_t(i,k) = (zrfln-zrfl(i)) . * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime . * RLVTT/RCPD zrfl(i) = zrfln ENDIF ENDDO ENDDO ENDIF c c La temperature de la premiere couche determine la pluie ou la neige: c DO i = 1, klon IF (todo(i)) THEN IF (t(i,1) .GT. RTT) THEN rain(i) = rain(i) + zrfl(i) ELSE snow(i) = snow(i) + zrfl(i) ENDIF ENDIF ENDDO c RETURN END SUBROUTINE kuofcl(pt, pq, pg, pp, LDCUM, kcbot) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927 c adaptation du code de Tiedtke du ECMWF c Objet: calculer le niveau de convection libre c (FCL: Free Convection Level) c====================================================================== c Arguments: c pt---input-R- temperature (K) c pq---input-R- vapeur d'eau (kg/kg) c pg---input-R- geopotentiel (g*z ou z est en metre) c pp---input-R- pression (Pa) c c LDCUM---output-L- Y-t-il la convection c kcbot---output-I- Niveau du bas de la convection c====================================================================== #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" #include "YOETHF.h" C REAL pt(klon,klev), pq(klon,klev), pg(klon,klev), pp(klon,klev) INTEGER kcbot(klon) LOGICAL LDCUM(klon) C REAL ztu(klon,klev), zqu(klon,klev), zlu(klon,klev) REAL zqold(klon), zbuo INTEGER is, i, k c c klab=1: on est sous le nuage convectif c klab=2: le bas du nuage convectif c klab=0: autres couches INTEGER klab(klon,klev) c c quand lflag=.true., on est sous le nuage, il faut donc appliquer c le processus d'elevation. LOGICAL lflag(klon) C DO k = 1, klev DO i = 1, klon ztu(i,k) = pt(i,k) zqu(i,k) = pq(i,k) zlu(i,k) = 0.0 klab(i,k) = 0 ENDDO ENDDO C---------------------------------------------------------------------- DO i = 1, klon klab(i,1)=1 kcbot(i)=2 LDCUM(i)=.FALSE. ENDDO C DO 290 k = 2, klev-1 c is=0 DO i = 1, klon if (klab(i,k-1).EQ.1) is = is + 1 lflag(i) = .FALSE. if (klab(i,k-1).EQ.1) lflag(i) = .TRUE. ENDDO IF (is.EQ.0) GOTO 290 c c on eleve le parcel d'air selon l'adiabatique sec c DO i = 1, klon IF (lflag(i)) THEN zqu(i,k) = zqu(i,k-1) ztu(i,k) = ztu(i,k-1) + (pg(i,k-1)-pg(i,k))/RCPD zbuo = ztu(i,k)*(1.+RETV*zqu(i,k))- . pt(i,k)*(1.+RETV*pq(i,k))+0.5 IF (zbuo.GT.0.) klab(i,k)=1 zqold(i) = zqu(i,k) ENDIF ENDDO c c on calcule la condensation eventuelle c CALL adjtq(pp(1,k), ztu(1,k), zqu(1,k), lflag, 1) c c s'il y a la condensation et la "buoyancy" force est positive c c'est bien le bas de la tour de convection c DO i=1, klon IF(lflag(i).AND.zqu(i,k).NE.zqold(i)) THEN klab(i,k) = 2 zlu(i,k) = zlu(i,k)+zqold(i)-zqu(i,k) zbuo = ztu(i,k)*(1.+RETV*zqu(i,k))- . pt(i,k)*(1.+RETV*pq(i,k))+0.5 IF (zbuo.GT.0.) THEN kcbot(i) = k LDCUM(i) = .TRUE. ENDIF ENDIF ENDDO C 290 CONTINUE C RETURN END SUBROUTINE adjtq(pp, pt, pq, LDFLAG, KCALL) IMPLICIT none c====================================================================== c Auteur(s): Z.X. Li (LMD/CNRS) date: 19940927 c adaptation du code de Tiedtke du ECMWF c Objet: ajustement entre T et Q c====================================================================== c Arguments: c pp---input-R- pression (Pa) c pt---input/output-R- temperature (K) c pq---input/output-R- vapeur d'eau (kg/kg) c====================================================================== C TO PRODUCE T,Q AND L VALUES FOR CLOUD ASCENT C C NOTE: INPUT PARAMETER KCALL DEFINES CALCULATION AS C KCALL=0 ENV. T AND QS IN*CUINI* C KCALL=1 CONDENSATION IN UPDRAFTS (E.G. CUBASE, CUASC) C KCALL=2 EVAPORATION IN DOWNDRAFTS (E.G. CUDLFS,CUDDRAF) C #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" C REAL pt(klon), pq(klon), pp(klon) LOGICAL ldflag(klon) INTEGER KCALL c REAL t_coup PARAMETER (t_coup=234.0) c REAL zcond(klon), zcond1 REAL zdelta, zcvm5, zldcp, zqsat, zcor, zdqsat INTEGER is, i #include "YOETHF.h" #include "FCTTRE.h" c DO i = 1, klon zcond(i) = 0.0 ENDDO C DO 210 i=1, klon IF (LDFLAG(i)) THEN zdelta=MAX(0.,SIGN(1.,RTT-pt(i))) zldcp = RLVTT*(1.-zdelta) + zdelta*RLSTT zldcp = zldcp / RCPD/(1.0+RVTMP2*pq(i)) IF (thermcep) THEN zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 / RCPD/(1.0+RVTMP2*pq(i)) zqsat=R2ES*FOEEW (pt(i), zdelta) / pp(i) zqsat=MIN(0.5,zqsat) zcor=1./(1.-RETV *zqsat) zqsat=zqsat*zcor zdqsat = FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor) ELSE IF (pt(i).LT.t_coup) THEN zqsat = qsats(pt(i))/pp(i) zdqsat = dqsats(pt(i),zqsat) ELSE zqsat = qsatl(pt(i))/pp(i) zdqsat = dqsatl(pt(i),zqsat) ENDIF ENDIF zcond(i)=(pq(i)-zqsat) / (1. + zdqsat) IF(KCALL.EQ.1) zcond(i)=MAX(zcond(i),0.) IF(KCALL.EQ.2) zcond(i)=MIN(zcond(i),0.) pt(i)=pt(i)+zldcp*zcond(i) pq(i)=pq(i)-zcond(i) ENDIF 210 CONTINUE C is = 0 DO i=1, klon if (zcond(i).NE.0.) is = is + 1 ENDDO IF(is.EQ.0) GOTO 230 C DO 220 i = 1, klon IF(LDFLAG(i).AND.zcond(i).NE.0.) THEN zdelta=MAX(0.,SIGN(1.,RTT-pt(i))) zldcp = RLVTT*(1.-zdelta) + zdelta*RLSTT zldcp = zldcp / RCPD/(1.0+RVTMP2*pq(i)) IF (thermcep) THEN zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 / RCPD/(1.0+RVTMP2*pq(i)) zqsat=R2ES*FOEEW (pt(i), zdelta) / pp(i) zqsat=MIN(0.5,zqsat) zcor=1./(1.-RETV *zqsat) zqsat=zqsat*zcor zdqsat = FOEDE(pt(i), zdelta, zcvm5, zqsat, zcor) ELSE IF (pt(i).LT.t_coup) THEN zqsat = qsats(pt(i))/pp(i) zdqsat = dqsats(pt(i),zqsat) ELSE zqsat = qsatl(pt(i))/pp(i) zdqsat = dqsatl(pt(i),zqsat) ENDIF ENDIF zcond1=(pq(i)-zqsat) / (1.+zdqsat) pt(i)=pt(i)+zldcp*zcond1 pq(i)=pq(i)-zcond1 END IF 220 CONTINUE C 230 CONTINUE RETURN END SUBROUTINE fiajh(dtime, paprs, pplay, t, q, . d_t, d_q, d_ql, rneb, . rain, snow, ibas, itop) IMPLICIT NONE c c Ajustement humide (Schema de convection de Manabe) C. #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c c Arguments: c REAL dtime ! intervalle du temps (s) REAL t(klon,klev) ! temperature (K) REAL q(klon,klev) ! humidite specifique (kg/kg) REAL paprs(klon,klev+1) ! pression a inter-couche (Pa) REAL pplay(klon,klev) ! pression au milieu de couche (Pa) c REAL d_t(klon,klev) ! incrementation pour la temperature REAL d_q(klon,klev) ! incrementation pour vapeur d'eau REAL d_ql(klon,klev) ! incrementation pour l'eau liquide REAL rneb(klon,klev) ! fraction nuageuse c REAL rain(klon) ! variable non utilisee REAL snow(klon) ! variable non utilisee INTEGER ibas(klon) ! variable non utilisee INTEGER itop(klon) ! variable non utilisee REAL t_coup PARAMETER (t_coup=234.0) REAL seuil_vap PARAMETER (seuil_vap=1.0E-10) c c Variables locales: c INTEGER i, k INTEGER k1, k1p, k2, k2p LOGICAL itest(klon) REAL delta_q(klon, klev) REAL cp_new_t(klev) REAL cp_delta_t(klev) REAL new_qb(klev) REAL v_cptj(klev), v_cptjk1, v_ssig REAL v_cptt(klon,klev), v_p, v_t REAL v_qs(klon,klev), v_qsd(klon,klev) REAL zq1(klon), zq2(klon) REAL gamcpdz(klon,2:klev) REAL zdp, zdpm c REAL zsat ! sur-saturation REAL zflo ! flotabilite c REAL local_q(klon,klev),local_t(klon,klev) c REAL zdelta, zcor, zcvm5 C #include "YOETHF.h" #include "FCTTRE.h" C DO k = 1, klev DO i = 1, klon local_q(i,k) = q(i,k) local_t(i,k) = t(i,k) rneb(i,k) = 0.0 d_ql(i,k) = 0.0 d_t(i,k) = 0.0 d_q(i,k) = 0.0 ENDDO ENDDO DO i = 1, klon rain(i) = 0.0 snow(i) = 0.0 ibas(i) = 0 itop(i) = 0 ENDDO c c Calculer v_qs et v_qsd: c DO k = 1, klev DO i = 1, klon v_cptt(i,k) = RCPD * local_t(i,k) v_t = local_t(i,k) v_p = pplay(i,k) c IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-v_t)) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*local_q(i,k)) v_qs(i,k)= R2ES * FOEEW(v_t,zdelta)/v_p v_qs(i,k)=MIN(0.5,v_qs(i,k)) zcor=1./(1.-RETV*v_qs(i,k)) v_qs(i,k)=v_qs(i,k)*zcor v_qsd(i,k) =FOEDE(v_t,zdelta,zcvm5,v_qs(i,k),zcor) ELSE IF (v_t.LT.t_coup) THEN v_qs(i,k) = qsats(v_t) / v_p v_qsd(i,k) = dqsats(v_t,v_qs(i,k)) ELSE v_qs(i,k) = qsatl(v_t) / v_p v_qsd(i,k) = dqsatl(v_t,v_qs(i,k)) ENDIF ENDIF ENDDO ENDDO c c Calculer Gamma * Cp * dz: (gamm est le gradient critique) c DO k = 2, klev DO i = 1, klon zdp = paprs(i,k)-paprs(i,k+1) zdpm = paprs(i,k-1)-paprs(i,k) gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) * . (v_cptt(i,k-1)*zdpm + v_cptt(i,k)*zdp) . +RLVTT /(zdpm+zdp) * . (v_qs(i,k-1)*zdpm + v_qs(i,k)*zdp) . )* (pplay(i,k-1)-pplay(i,k)) / paprs(i,k) ) . / (1.0+(v_qsd(i,k-1)*zdpm+ . v_qsd(i,k)*zdp)/(zdpm+zdp) ) ENDDO ENDDO C C------------------------------------ modification des profils instables DO 9999 i = 1, klon itest(i) = .FALSE. C k1 = 0 k2 = 1 C 810 CONTINUE ! chercher k1, le bas de la colonne k2 = k2 + 1 IF (k2 .GT. klev) GOTO 9999 zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2) zsat=(local_q(i,k2-1)-v_qs(i,k2-1))*(paprs(i,k2-1)-paprs(i,k2)) . +(local_q(i,k2)-v_qs(i,k2))*(paprs(i,k2)-paprs(i,k2+1)) IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810 k1 = k2 - 1 itest(i) = .TRUE. C 820 CONTINUE ! chercher k2, le haut de la colonne IF (k2 .EQ. klev) GOTO 821 k2p = k2 + 1 zsat=zsat +(paprs(i,k2p)-paprs(i,k2p+1)) . *(local_q(i,k2p)-v_qs(i,k2p)) zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p) IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821 k2 = k2p GOTO 820 821 CONTINUE C C------------------------------------------------------ ajustement local 830 CONTINUE ! ajustement proprement dit v_cptj(k1) = 0.0 zdp = paprs(i,k1)-paprs(i,k1+1) v_cptjk1 = ( (1.0+v_qsd(i,k1))*(v_cptt(i,k1)+v_cptj(k1)) . + RLVTT*(local_q(i,k1)-v_qs(i,k1)) ) * zdp v_ssig = zdp * (1.0+v_qsd(i,k1)) C k1p = k1 + 1 DO k = k1p, k2 zdp = paprs(i,k)-paprs(i,k+1) v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k) v_cptjk1 = v_cptjk1 + zdp . * ( (1.0+v_qsd(i, k))*(v_cptt(i,k)+v_cptj(k)) . + RLVTT*(local_q(i,k)-v_qs(i,k)) ) v_ssig = v_ssig + zdp *(1.0+v_qsd(i,k)) ENDDO C DO k = k1, k2 cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k) cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k) new_qb(k) = v_qs(i,k) + v_qsd(i,k)*cp_delta_t(k)/RLVTT local_q(i,k) = new_qb(k) local_t(i,k) = cp_new_t(k) / RCPD ENDDO C C--------------------------------------------------- sondage vers le bas C -- on redefinit les variables prognostiques dans C -- la colonne qui vient d'etre ajustee C DO k = k1, k2 v_cptt(i,k) = RCPD * local_t(i,k) v_t = local_t(i,k) v_p = pplay(i,k) C IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-v_t)) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*local_q(i,k)) v_qs(i,k)= R2ES * FOEEW(v_t,zdelta)/v_p v_qs(i,k)=MIN(0.5,v_qs(i,k)) zcor=1./(1.-RETV*v_qs(i,k)) v_qs(i,k)=v_qs(i,k)*zcor v_qsd(i,k) =FOEDE(v_t,zdelta,zcvm5,v_qs(i,k),zcor) ELSE IF (v_t.LT.t_coup) THEN v_qs(i,k) = qsats(v_t) / v_p v_qsd(i,k) = dqsats(v_t,v_qs(i,k)) ELSE v_qs(i,k) = qsatl(v_t) / v_p v_qsd(i,k) = dqsatl(v_t,v_qs(i,k)) ENDIF ENDIF ENDDO DO k = 2, klev zdpm = paprs(i,k-1) - paprs(i,k) zdp = paprs(i,k) - paprs(i,k+1) gamcpdz(i,k) = ( ( RD/RCPD /(zdpm+zdp) * . (v_cptt(i,k-1)*zdpm+v_cptt(i,k)*zdp) . +RLVTT /(zdpm+zdp) * . (v_qs(i,k-1)*zdpm+v_qs(i,k)*zdp) . )* (pplay(i,k-1)-pplay(i,k)) / paprs(i,k) ) . / (1.0+(v_qsd(i,k-1)*zdpm+v_qsd(i,k)*zdp) . /(zdpm+zdp) ) ENDDO C C Verifier si l'on peut etendre la colonne vers le bas C IF (k1 .EQ. 1) GOTO 841 ! extension echouee zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) zsat=(local_q(i,k1-1)-v_qs(i,k1-1))*(paprs(i,k1-1)-paprs(i,k1)) . + (local_q(i,k1)-v_qs(i,k1))*(paprs(i,k1)-paprs(i,k1+1)) IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! extension echouee C 840 CONTINUE k1 = k1 - 1 IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) zsat = zsat + (local_q(i,k1-1)-v_qs(i,k1-1)) . *(paprs(i,k1-1)-paprs(i,k1)) zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1) IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN GOTO 840 ELSE GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995) ENDIF 841 CONTINUE C GOTO 810 ! chercher d'autres blocks en haut C 9999 CONTINUE ! boucle sur tous les points C----------------------------------------------------------------------- c c Determiner la fraction nuageuse (hypothese: la nebulosite a lieu c a l'endroit ou la vapeur d'eau est diminuee par l'ajustement): c DO k = 1, klev DO i = 1, klon IF (itest(i)) THEN delta_q(i,k) = local_q(i,k) - q(i,k) IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0 ENDIF ENDDO ENDDO c c Distribuer l'eau condensee en eau liquide nuageuse (hypothese: c l'eau liquide est distribuee aux endroits ou la vapeur d'eau c diminue et d'une maniere proportionnelle a cet diminution): c DO i = 1, klon IF (itest(i)) THEN zq1(i) = 0.0 zq2(i) = 0.0 ENDIF ENDDO DO k = 1, klev DO i = 1, klon IF (itest(i)) THEN zdp = paprs(i,k)-paprs(i,k+1) zq1(i) = zq1(i) - delta_q(i,k) * zdp zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp ENDIF ENDDO ENDDO DO k = 1, klev DO i = 1, klon IF (itest(i)) THEN IF (zq2(i).NE.0.0) . d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i) ENDIF ENDDO ENDDO C DO k = 1, klev DO i = 1, klon local_q(i, k) = MAX(local_q(i, k), seuil_vap) ENDDO ENDDO C DO k = 1, klev DO i = 1, klon d_t(i,k) = local_t(i,k) - t(i,k) d_q(i,k) = local_q(i,k) - q(i,k) ENDDO ENDDO c RETURN END SUBROUTINE fiajc(dtime,paprs,pplay, . t, q,conv_q, . d_t, d_q, d_ql,rneb, . rain, snow, ibas, itop) IMPLICIT NONE c #include "dimensions.h" #include "dimphy.h" #include "YOMCST.h" c c Options: c INTEGER plb ! niveau de depart pour la convection PARAMETER (plb=4) c c Mystere: cette option n'est pas innocente pour les resultats ! c Qui peut resoudre ce mystere ? (Z.X.Li mars 1995) LOGICAL vector ! calcul vectorise PARAMETER (vector=.FALSE.) c REAL t_coup PARAMETER (t_coup=234.0) c c Arguments: c REAL q(klon,klev) ! humidite specifique (kg/kg) REAL t(klon,klev) ! temperature (K) REAL paprs(klon,klev+1) ! pression a inter-couche (Pa) REAL pplay(klon,klev) ! pression au milieu de couche (Pa) REAL dtime ! intervalle du temps (s) REAL conv_q(klon,klev) ! taux de convergence de l'humidite REAL rneb(klon,klev) ! fraction nuageuse REAL d_q(klon,klev) ! incrementaion pour la vapeur d'eau REAL d_ql(klon,klev) ! incrementation pour l'eau liquide REAL d_t(klon,klev) ! incrementation pour la temperature REAL rain(klon) ! variable non-utilisee REAL snow(klon) ! variable non-utilisee INTEGER itop(klon) ! variable non-utilisee INTEGER ibas(klon) ! variable non-utilisee c INTEGER kh(klon), i, k LOGICAL nuage(klon), test(klon,klev) REAL zconv(klon), zdeh(klon,klev), zvirt(klon) REAL zdqs(klon,klev), zqs(klon,klev) REAL ztt, zvar, zfrac(klon) REAL zq1(klon), zq2(klon) REAL zdelta, zcor, zcvm5 C #include "YOETHF.h" #include "FCTTRE.h" c c Initialiser les sorties: c DO k = 1, klev DO i = 1, klon rneb(i,k) = 0.0 d_ql(i,k) = 0.0 d_t(i,k) = 0.0 d_q(i,k) = 0.0 ENDDO ENDDO DO i = 1, klon itop(i) = 0 ibas(i) = 0 rain(i) = 0.0 snow(i) = 0.0 ENDDO c c Calculer Qs et L/Cp * dQs/dT: c DO k = 1, klev DO i = 1, klon ztt = t(i,k) IF (thermcep) THEN zdelta=MAX(0.,SIGN(1.,RTT-ztt)) zcvm5=R5LES*RLVTT*(1.-zdelta) + zdelta*R5IES*RLSTT zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*q(i,k)) zqs(i,k)= R2ES*FOEEW(ztt,zdelta)/pplay(i,k) zqs(i,k)=MIN(0.5,zqs(i,k)) zcor=1./(1.-RETV*zqs(i,k)) zqs(i,k)=zqs(i,k)*zcor zdqs(i,k) =FOEDE(ztt,zdelta,zcvm5,zqs(i,k),zcor) ELSE IF (ztt .LT. t_coup) THEN zqs(i,k) = qsats(ztt) / pplay(i,k) zdqs(i,k) = dqsats(ztt,zqs(i,k)) ELSE zqs(i,k) = qsatl(ztt) / pplay(i,k) zdqs(i,k) = dqsatl(ztt,zqs(i,k)) ENDIF ENDIF ENDDO ENDDO c c Determiner la difference de l'energie totale saturee: c DO i = 1, klon k = plb zdeh(i,k) = RCPD * (t(i,k-1)-t(i,k)) . - RD *0.5*(t(i,k-1)+t(i,k))/paprs(i,k) . *(pplay(i,k-1)-pplay(i,k)) . + RLVTT*(zqs(i,k-1)-zqs(i,k)) zdeh(i,k) = zdeh(i,k) * 0.5 ! on prend la moitie ENDDO DO k = plb+1, klev DO i = 1, klon zdeh(i,k) = zdeh(i,k-1) . + RCPD * (t(i,k-1)-t(i,k)) . - RD *0.5*(t(i,k-1)+t(i,k))/paprs(i,k) . *(pplay(i,k-1)-pplay(i,k)) . + RLVTT*(zqs(i,k-1)-zqs(i,k)) ENDDO ENDDO c c Determiner le sommet du nuage selon l'instabilite c Calculer les convergences d'humidite (reelle et virtuelle) c DO i = 1, klon nuage(i) = .TRUE. zconv(i) = 0.0 zvirt(i) = 0.0 kh(i) = -999 ENDDO DO k = plb, klev DO i = 1, klon nuage(i) = nuage(i) .AND. zdeh(i,k).GT.0.0 IF (nuage(i)) THEN kh(i) = k zconv(i) = zconv(i)+conv_q(i,k)*dtime . *(paprs(i,k)-paprs(i,k+1)) zvirt(i)=zvirt(i)+(zdeh(i,k)/RLVTT+zqs(i,k)-q(i,k)) . *(paprs(i,k)-paprs(i,k+1)) ENDIF ENDDO ENDDO c IF (vector) THEN c c DO k = plb, klev DO i = 1, klon IF (k.LE.kh(i) .AND. kh(i).GT.plb .AND. zconv(i).GT.0.0) THEN test(i,k) = .TRUE. zfrac(i) = MAX(0.0,MIN(zconv(i)/zvirt(i),1.0)) ELSE test(i,k) = .FALSE. ENDIF ENDDO ENDDO c DO k = plb, klev DO i = 1, klon IF (test(i,k)) THEN zvar = zdeh(i,k)/(1.0+zdqs(i,k)) d_q(i,k) = (zvar*zdqs(i,k)/RLVTT+zqs(i,k)-q(i,k))*zfrac(i) . - conv_q(i,k)*dtime d_t(i,k) = zvar * zfrac(i) / RCPD ENDIF ENDDO ENDDO c DO i = 1, klon zq1(i) = 0.0 zq2(i) = 0.0 ENDDO DO k = plb, klev DO i = 1, klon IF (test(i,k)) THEN IF (d_q(i,k).LT.0.0) rneb(i,k) = zfrac(i) zq1(i) = zq1(i) - d_q(i,k) * (paprs(i,k)-paprs(i,k+1)) zq2(i) = zq2(i) - MIN(0.0, d_q(i,k)) . * (paprs(i,k)-paprs(i,k+1)) ENDIF ENDDO ENDDO c DO k = plb, klev DO i = 1, klon IF (test(i,k)) THEN IF(zq2(i).NE.0.)d_ql(i,k)=-MIN(0.0,d_q(i,k))*zq1(i)/zq2(i) ENDIF ENDDO ENDDO c ELSE ! (.NOT. vector) c DO 999 i = 1, klon IF (kh(i).GT.plb .AND. zconv(i).GT.0.0) THEN ccc IF (kh(i).LE.plb) GOTO 999 ! il n'y a pas d'instabilite ccc IF (zconv(i).LE.0.0) GOTO 999 ! convergence insuffisante zfrac(i) = MAX(0.0,MIN(zconv(i)/zvirt(i),1.0)) DO k = plb, kh(i) zvar = zdeh(i,k)/(1.0+zdqs(i,k)) d_q(i,k) = (zvar*zdqs(i,k)/RLVTT+zqs(i,k)-q(i,k))*zfrac(i) . - conv_q(i,k)*dtime d_t(i,k) = zvar * zfrac(i) / RCPD ENDDO c zq1(i) = 0.0 zq2(i) = 0.0 DO k = plb, kh(i) IF (d_q(i,k).LT.0.0) rneb(i,k) = zfrac(i) zq1(i) = zq1(i) - d_q(i,k) * (paprs(i,k)-paprs(i,k+1)) zq2(i) = zq2(i) - MIN(0.0, d_q(i,k)) . * (paprs(i,k)-paprs(i,k+1)) ENDDO DO k = plb, kh(i) IF(zq2(i).NE.0.)d_ql(i,k)=-MIN(0.0,d_q(i,k))*zq1(i)/zq2(i) ENDDO ENDIF 999 CONTINUE c ENDIF ! fin de teste sur vector c RETURN END