/[lmdze]/trunk/phylmd/fisrtilp.f
ViewVC logotype

Diff of /trunk/phylmd/fisrtilp.f

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/libf/phylmd/fisrtilp.f90 revision 52 by guez, Fri Sep 23 12:28:01 2011 UTC trunk/phylmd/fisrtilp.f90 revision 76 by guez, Fri Nov 15 18:45:49 2013 UTC
# Line 1  Line 1 
1  SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs,d_t,d_q,d_ql,rneb, &  module fisrtilp_m
      radliq,rain,snow,pfrac_impa,pfrac_nucl,pfrac_1nucl,frac_impa, &  
      frac_nucl,prfl,psfl,rhcl)  
   
   ! From phylmd/fisrtilp.F,v 1.2 2004/11/09 16:55:40  
   ! Auteur(s): Z.X. Li (LMD/CNRS)  
   ! Date: le 20 mars 1995  
   ! Objet: condensation et precipitation stratiforme.  
   !        schema de nuage  
   
   USE dimens_m  
   USE dimphy  
   USE tracstoke  
   USE suphec_m  
   USE yoethf_m  
   USE fcttre  
   USE comfisrtilp  
   use numer_rec, only: nr_erf  
2    
3    IMPLICIT NONE    IMPLICIT NONE
4    
5    ! Arguments:  contains
6    
7    REAL, INTENT (IN) :: & ! intervalle du temps (s)                    SUBROUTINE fisrtilp(dtime, paprs, pplay, t, q, ptconv, ratqs, d_t, d_q, &
8         dtime         d_ql, rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, &
9    REAL, INTENT (IN) :: paprs(klon,klev+1) ! pression a inter-couche           frac_impa, frac_nucl, prfl, psfl, rhcl)
10    REAL, INTENT (IN):: pplay(klon,klev) ! pression au milieu de couche  
11    REAL, INTENT (IN):: t(klon,klev) ! temperature (K)      ! From phylmd/fisrtilp.F, version 1.2 2004/11/09 16:55:40
12    REAL q(klon,klev) ! humidite specifique (kg/kg)                        ! Author: Z. X. Li (LMD/CNRS), 20 mars 1995
13    REAL d_t(klon,klev) ! incrementation de la temperature (K)          
14    REAL d_q(klon,klev) ! incrementation de la vapeur d'eau                ! Objet : condensation et précipitation stratiforme, schéma de
15    REAL d_ql(klon,klev) ! incrementation de l'eau liquide                  ! nuage, schéma de condensation à grande échelle (pluie).
16    REAL rneb(klon,klev) ! fraction nuageuse                            
17    REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements        USE dimphy, ONLY: klev, klon
18    REAL rhcl(klon,klev) ! humidite relative en ciel clair                  USE suphec_m, ONLY: rcpd, rd, retv, rg, rlstt, rlvtt, rtt
19    REAL rain(klon) ! pluies (mm/s)                                        USE yoethf_m, ONLY: r2es, r5ies, r5les, rvtmp2
20    REAL snow(klon) ! neige (mm/s)                                          USE fcttre, ONLY: dqsatl, dqsats, foede, foeew, qsatl, qsats, thermcep
21    REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)      USE comfisrtilp, ONLY: cld_lc_con, cld_lc_lsc, cld_tau_con, &
22    REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)           cld_tau_lsc, coef_eva, ffallv_con, ffallv_lsc, iflag_pdf, reevap_ice
23    ! Coeffients de fraction lessivee : pour OFF-LINE      USE numer_rec_95, ONLY: nr_erf
24    
25    REAL pfrac_nucl(klon,klev)      ! Arguments:
26    REAL pfrac_1nucl(klon,klev)  
27    REAL pfrac_impa(klon,klev)      REAL, INTENT (IN):: dtime ! intervalle du temps (s)                
28        REAL, INTENT (IN):: paprs(klon, klev+1) ! pression a inter-couche  
29    ! Fraction d'aerosols lessivee par impaction et par nucleation      REAL, INTENT (IN):: pplay(klon, klev) ! pression au milieu de couche
30    ! POur ON-LINE      REAL, INTENT (IN):: t(klon, klev) ! temperature (K)
31        REAL, INTENT (IN):: q(klon, klev) ! humidite specifique (kg/kg)
32    REAL frac_impa(klon,klev)      LOGICAL ptconv(klon, klev) ! determine la largeur de distribution de vapeur
33    REAL frac_nucl(klon,klev)      REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur
34    REAL zct(klon), zcl(klon)      REAL d_t(klon, klev) ! incrementation de la temperature (K)        
35    !AA      REAL d_q(klon, klev) ! incrementation de la vapeur d'eau          
36        REAL d_ql(klon, klev) ! incrementation de l'eau liquide            
37    ! Options du programme:      REAL rneb(klon, klev) ! fraction nuageuse                          
38        REAL radliq(klon, klev) ! eau liquide utilisee dans rayonnements  
39    REAL seuil_neb ! un nuage existe vraiment au-dela      REAL rain(klon) ! pluies (mm/s)                                  
40    PARAMETER (seuil_neb=0.001)      REAL snow(klon) ! neige (mm/s)                                    
41    
42    INTEGER ninter ! sous-intervals pour la precipitation      ! Coeffients de fraction lessivee : pour OFF-LINE
43    PARAMETER (ninter=5)      REAL pfrac_impa(klon, klev)
44    LOGICAL evap_prec ! evaporation de la pluie                            REAL pfrac_nucl(klon, klev)
45    PARAMETER (evap_prec=.TRUE.)      REAL pfrac_1nucl(klon, klev)
46    REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur  
47    LOGICAL ptconv(klon,klev) ! determine la largeur de distribution de vapeur      ! Fraction d'aerosols lessivee par impaction et par nucleation
48    REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)      ! POur ON-LINE
49    REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)      REAL frac_nucl(klon, klev)
50    
51    LOGICAL cpartiel ! condensation partielle                              REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
52    PARAMETER (cpartiel=.TRUE.)      REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
53    REAL t_coup      REAL rhcl(klon, klev) ! humidite relative en ciel clair            
54    PARAMETER (t_coup=234.0)  
55        ! Local:
56    ! Variables locales:  
57        ! Fraction d'aerosols lessivee par impaction et par nucleation
58    INTEGER i, k, n, kk      ! POur ON-LINE
59    REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5      REAL frac_impa(klon, klev)
60    REAL zrfl(klon), zrfln(klon), zqev, zqevt      REAL zct(klon), zcl(klon)
61    REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq      !AA
62    REAL ztglace, zt(klon)  
63    INTEGER nexpo ! exponentiel pour glace/eau                              ! Options du programme:
64    REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)  
65    REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)      REAL seuil_neb ! un nuage existe vraiment au-dela
66        PARAMETER (seuil_neb=0.001)
67    LOGICAL appel1er  
68    SAVE appel1er      INTEGER ninter ! sous-intervals pour la precipitation
69        PARAMETER (ninter=5)
70    !---------------------------------------------------------------      LOGICAL evap_prec ! evaporation de la pluie                      
71        PARAMETER (evap_prec=.TRUE.)
72    !AA Variables traceurs:      REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
73    !AA  Provisoire !!! Parametres alpha du lessivage      REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)
74    !AA  A priori on a 4 scavenging numbers possibles  
75        LOGICAL cpartiel ! condensation partielle                        
76    REAL a_tr_sca(4)      PARAMETER (cpartiel=.TRUE.)
77    SAVE a_tr_sca      REAL t_coup
78        PARAMETER (t_coup=234.0)
79    ! Variables intermediaires  
80        ! Variables locales:
81    REAL zalpha_tr  
82    REAL zfrac_lessi      INTEGER i, k, n, kk
83    REAL zprec_cond(klon)      REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
84    !AA      REAL zrfl(klon), zrfln(klon), zqev, zqevt
85    REAL zmair, zcpair, zcpeau      REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
86    !     Pour la conversion eau-neige      REAL ztglace, zt(klon)
87    REAL zlh_solid(klon), zm_solid      INTEGER nexpo ! exponentiel pour glace/eau                        
88    !IM      REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
89    INTEGER klevm1      REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
90    !---------------------------------------------------------------  
91        LOGICAL appel1er
92    ! Fonctions en ligne:      SAVE appel1er
93    
94    REAL fallvs, fallvc ! vitesse de chute pour crystaux de glace            !---------------------------------------------------------------
95    REAL zzz  
96        !AA Variables traceurs:
97    fallvc(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_con      !AA  Provisoire !!! Parametres alpha du lessivage
98    fallvs(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_lsc      !AA  A priori on a 4 scavenging numbers possibles
99    
100    DATA appel1er/ .TRUE./      REAL a_tr_sca(4)
101    !ym      SAVE a_tr_sca
102    zdelq = 0.0  
103        ! Variables intermediaires
104    IF (appel1er) THEN  
105        REAL zalpha_tr
106       PRINT *, 'fisrtilp, ninter:', ninter      REAL zfrac_lessi
107       PRINT *, 'fisrtilp, evap_prec:', evap_prec      REAL zprec_cond(klon)
108       PRINT *, 'fisrtilp, cpartiel:', cpartiel      !AA
109       IF (abs(dtime/float(ninter)-360.0)>0.001) THEN      REAL zmair, zcpair, zcpeau
110          PRINT *, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime      !     Pour la conversion eau-neige
111          PRINT *, 'Je prefere un sous-intervalle de 6 minutes'      REAL zlh_solid(klon), zm_solid
112          !         stop 1      !IM
113       END IF      INTEGER klevm1
114       appel1er = .FALSE.      !---------------------------------------------------------------
115    
116       !AA initialiation provisoire      ! Fonctions en ligne:
117       a_tr_sca(1) = -0.5  
118       a_tr_sca(2) = -0.5      REAL fallvs, fallvc ! vitesse de chute pour crystaux de glace      
119       a_tr_sca(3) = -0.5      REAL zzz
120       a_tr_sca(4) = -0.5  
121        fallvc(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_con
122       !AA Initialisation a 1 des coefs des fractions lessivees      fallvs(zzz) = 3.29/2.0*((zzz)**0.16)*ffallv_lsc
123    
124       DO k = 1, klev      DATA appel1er/ .TRUE./
125          DO i = 1, klon      !ym
126             pfrac_nucl(i,k) = 1.      zdelq = 0.0
127             pfrac_1nucl(i,k) = 1.  
128             pfrac_impa(i,k) = 1.      IF (appel1er) THEN
129          END DO  
130       END DO         PRINT *, 'fisrtilp, ninter:', ninter
131           PRINT *, 'fisrtilp, evap_prec:', evap_prec
132           PRINT *, 'fisrtilp, cpartiel:', cpartiel
133    END IF !  test sur appel1er         IF (abs(dtime / real(ninter) - 360.) > 0.001) THEN
134    !MAf Initialisation a 0 de zoliq            PRINT *, "fisrtilp : ce n'est pas prévu, voir Z. X. Li", dtime
135    DO i = 1, klon            PRINT *, 'Je préfère un sous-intervalle de 6 minutes.'
136       zoliq(i) = 0.         END IF
137    END DO         appel1er = .FALSE.
138    ! Determiner les nuages froids par leur temperature  
139    !  nexpo regle la raideur de la transition eau liquide / eau glace.         !AA initialiation provisoire
140           a_tr_sca(1) = -0.5
141    ztglace = rtt - 15.0         a_tr_sca(2) = -0.5
142    nexpo = 6         a_tr_sca(3) = -0.5
143    !cc      nexpo = 1         a_tr_sca(4) = -0.5
144    
145    ! Initialiser les sorties:         !AA Initialisation a 1 des coefs des fractions lessivees
146    
147    DO k = 1, klev + 1         DO k = 1, klev
148       DO i = 1, klon            DO i = 1, klon
149          prfl(i,k) = 0.0               pfrac_nucl(i, k) = 1.
150          psfl(i,k) = 0.0               pfrac_1nucl(i, k) = 1.
151       END DO               pfrac_impa(i, k) = 1.
152    END DO            END DO
153           END DO
154    DO k = 1, klev  
155       DO i = 1, klon  
156          d_t(i,k) = 0.0      END IF !  test sur appel1er
157          d_q(i,k) = 0.0      !MAf Initialisation a 0 de zoliq
158          d_ql(i,k) = 0.0      DO i = 1, klon
159          rneb(i,k) = 0.0         zoliq(i) = 0.
160          radliq(i,k) = 0.0      END DO
161          frac_nucl(i,k) = 1.      ! Determiner les nuages froids par leur temperature
162          frac_impa(i,k) = 1.      !  nexpo regle la raideur de la transition eau liquide / eau glace.
163       END DO  
164    END DO      ztglace = rtt - 15.0
165    DO i = 1, klon      nexpo = 6
166       rain(i) = 0.0      !cc      nexpo = 1
167       snow(i) = 0.0  
168    END DO      ! Initialiser les sorties:
169    
170    ! Initialiser le flux de precipitation a zero      DO k = 1, klev + 1
171           DO i = 1, klon
172    DO i = 1, klon            prfl(i, k) = 0.0
173       zrfl(i) = 0.0            psfl(i, k) = 0.0
174       zneb(i) = seuil_neb         END DO
175    END DO      END DO
176    
177        DO k = 1, klev
178    !AA Pour plus de securite         DO i = 1, klon
179              d_t(i, k) = 0.0
180    zalpha_tr = 0.            d_q(i, k) = 0.0
181    zfrac_lessi = 0.            d_ql(i, k) = 0.0
182              rneb(i, k) = 0.0
183    !AA----------------------------------------------------------            radliq(i, k) = 0.0
184              frac_nucl(i, k) = 1.
185    ! Boucle verticale (du haut vers le bas)            frac_impa(i, k) = 1.
186           END DO
187    !IM : klevm1      END DO
188    klevm1 = klev - 1      DO i = 1, klon
189    DO  k = klev, 1, -1         rain(i) = 0.0
190           snow(i) = 0.0
191       !AA----------------------------------------------------------      END DO
192    
193       DO i = 1, klon      ! Initialiser le flux de precipitation a zero
194          zt(i) = t(i,k)  
195          zq(i) = q(i,k)      DO i = 1, klon
196       END DO         zrfl(i) = 0.0
197           zneb(i) = seuil_neb
198       ! Calculer la varition de temp. de l'air du a la chaleur sensible      END DO
199       ! transporter par la pluie.  
200       ! Il resterait a rajouter cet effet de la chaleur sensible sur les  
201       ! flux de surface, du a la diff. de temp. entre le 1er niveau et la      !AA Pour plus de securite
202       ! surface.  
203        zalpha_tr = 0.
204       DO i = 1, klon      zfrac_lessi = 0.
205          IF (k<=klevm1) THEN  
206             zmair = (paprs(i,k)-paprs(i,k+1))/rg      !AA----------------------------------------------------------
207             zcpair = rcpd*(1.0+rvtmp2*zq(i))  
208             zcpeau = rcpd*rvtmp2      ! Boucle verticale (du haut vers le bas)
209             zt(i) = ((t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau+zmair*zcpair* &  
210                  zt(i))/(zmair*zcpair+zrfl(i)*dtime*zcpeau)      !IM : klevm1
211             !C        WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))      klevm1 = klev - 1
212          END IF      DO  k = klev, 1, -1
213       END DO  
214           !AA----------------------------------------------------------
215       ! Calculer l'evaporation de la precipitation  
216           DO i = 1, klon
217              zt(i) = t(i, k)
218              zq(i) = q(i, k)
219       IF (evap_prec) THEN         END DO
220          DO i = 1, klon  
221             IF (zrfl(i)>0.) THEN         ! Calculer la varition de temp. de l'air du a la chaleur sensible
222                IF (thermcep) THEN         ! transporter par la pluie.
223                   zdelta = max(0.,sign(1.,rtt-zt(i)))         ! Il resterait a rajouter cet effet de la chaleur sensible sur les
224                   zqs(i) = r2es*foeew(zt(i),zdelta)/pplay(i,k)         ! flux de surface, du a la diff. de temp. entre le 1er niveau et la
225                   zqs(i) = min(0.5,zqs(i))         ! surface.
226                   zcor = 1./(1.-retv*zqs(i))  
227                   zqs(i) = zqs(i)*zcor         DO i = 1, klon
228                ELSE            IF (k<=klevm1) THEN
229                   IF (zt(i)<t_coup) THEN               zmair = (paprs(i, k)-paprs(i, k+1))/rg
230                      zqs(i) = qsats(zt(i))/pplay(i,k)               zcpair = rcpd*(1.0+rvtmp2*zq(i))
231                   ELSE               zcpeau = rcpd*rvtmp2
232                      zqs(i) = qsatl(zt(i))/pplay(i,k)               zt(i) = ((t(i, k + 1) + d_t(i, k + 1)) * zrfl(i) * dtime &
233                   END IF                    * zcpeau + zmair * zcpair* zt(i)) &
234                END IF                    / (zmair * zcpair + zrfl(i) * dtime * zcpeau)
235                zqev = max(0.0,(zqs(i)-zq(i))*zneb(i))            END IF
236                zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &         END DO
237                     (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*rd/rg  
238                zqevt = max(0.0,min(zqevt,zrfl(i)))*rg*dtime/ &         ! Calculer l'evaporation de la precipitation
239                     (paprs(i,k)-paprs(i,k+1))  
240                zqev = min(zqev,zqevt)  
241                zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime  
242           IF (evap_prec) THEN
243                ! pour la glace, on réévapore toute la précip dans la couche du dessous            DO i = 1, klon
244                ! la glace venant de la couche du dessus est simplement dans la couche               IF (zrfl(i)>0.) THEN
245                ! du dessous.                  IF (thermcep) THEN
246                       zdelta = max(0., sign(1., rtt-zt(i)))
247                IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.                     zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
248                       zqs(i) = min(0.5, zqs(i))
249                zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &                     zcor = 1./(1.-retv*zqs(i))
250                     k+1)))*dtime                     zqs(i) = zqs(i)*zcor
251                zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &                  ELSE
252                     k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))                     IF (zt(i)<t_coup) THEN
253                zrfl(i) = zrfln(i)                        zqs(i) = qsats(zt(i))/pplay(i, k)
254             END IF                     ELSE
255          END DO                        zqs(i) = qsatl(zt(i))/pplay(i, k)
256       END IF                     END IF
257                    END IF
258       ! Calculer Qs et L/Cp*dQs/dT:                  zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
259                    zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
260       IF (thermcep) THEN                       (paprs(i, k)-paprs(i, k+1))/pplay(i, k)*zt(i)*rd/rg
261          DO i = 1, klon                  zqevt = max(0.0, min(zqevt, zrfl(i)))*rg*dtime/ &
262             zdelta = max(0.,sign(1.,rtt-zt(i)))                       (paprs(i, k)-paprs(i, k+1))
263             zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta                  zqev = min(zqev, zqevt)
264             zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))                  zrfln(i) = zrfl(i) - zqev*(paprs(i, k)-paprs(i, k+1))/rg/dtime
265             zqs(i) = r2es*foeew(zt(i),zdelta)/pplay(i,k)  
266             zqs(i) = min(0.5,zqs(i))                  ! pour la glace, on réévapore toute la précip dans la
267             zcor = 1./(1.-retv*zqs(i))                  ! couche du dessous la glace venant de la couche du
268             zqs(i) = zqs(i)*zcor                  ! dessus est simplement dans la couche du dessous.
269             zdqs(i) = foede(zt(i),zdelta,zcvm5,zqs(i),zcor)  
270          END DO                  IF (zt(i)<t_coup .AND. reevap_ice) zrfln(i) = 0.
271       ELSE  
272          DO i = 1, klon                  zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
273             IF (zt(i)<t_coup) THEN                       k+1)))*dtime
274                zqs(i) = qsats(zt(i))/pplay(i,k)                  zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i, k)-paprs(i, &
275                zdqs(i) = dqsats(zt(i),zqs(i))                       k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
276             ELSE                  zrfl(i) = zrfln(i)
277                zqs(i) = qsatl(zt(i))/pplay(i,k)               END IF
278                zdqs(i) = dqsatl(zt(i),zqs(i))            END DO
279             END IF         END IF
280          END DO  
281       END IF         ! Calculer Qs et L/Cp*dQs/dT:
282    
283       ! Determiner la condensation partielle et calculer la quantite         IF (thermcep) THEN
284       ! de l'eau condensee:            DO i = 1, klon
285                 zdelta = max(0., sign(1., rtt-zt(i)))
286       IF (cpartiel) THEN               zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
287                 zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
288          !        print*,'Dans partiel k=',k               zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
289                 zqs(i) = min(0.5, zqs(i))
290          !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau               zcor = 1./(1.-retv*zqs(i))
291          !   nuageuse a partir des PDF de Sandrine Bony.               zqs(i) = zqs(i)*zcor
292          !   rneb  : fraction nuageuse               zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
293          !   zqn   : eau totale dans le nuage            END DO
294          !   zcond : eau condensee moyenne dans la maille.         ELSE
295          !           on prend en compte le réchauffement qui diminue la partie condensee            DO i = 1, klon
296                 IF (zt(i)<t_coup) THEN
297          !   Version avec les raqts                  zqs(i) = qsats(zt(i))/pplay(i, k)
298                    zdqs(i) = dqsats(zt(i), zqs(i))
299          IF (iflag_pdf==0) THEN               ELSE
300                    zqs(i) = qsatl(zt(i))/pplay(i, k)
301             DO i = 1, klon                  zdqs(i) = dqsatl(zt(i), zqs(i))
302                zdelq = min(ratqs(i,k),0.99)*zq(i)               END IF
303                rneb(i,k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)            END DO
304                zqn(i) = (zq(i)+zdelq+zqs(i))/2.0         END IF
305             END DO  
306           ! Determiner la condensation partielle et calculer la quantite
307          ELSE         ! de l'eau condensee:
308    
309             !   Version avec les nouvelles PDFs.         IF (cpartiel) THEN
310             DO i = 1, klon  
311                IF (zq(i)<1.E-15) THEN            !        print*, 'Dans partiel k=', k
312                   !C Lionel GUEZ                print*,'ZQ(',i,',',k,')=',zq(i)  
313                   zq(i) = 1.E-15            !   Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
314                END IF            !   nuageuse a partir des PDF de Sandrine Bony.
315             END DO            !   rneb  : fraction nuageuse
316             DO i = 1, klon            !   zqn   : eau totale dans le nuage
317                zpdf_sig(i) = ratqs(i,k)*zq(i)            !   zcond : eau condensee moyenne dans la maille.
318                zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))  
319                zpdf_delta(i) = log(zq(i)/zqs(i))            !           on prend en compte le réchauffement qui diminue
320                zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))            !           la partie condensee
321                zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))  
322                zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)            !   Version avec les raqts
323                zpdf_e1(i) = sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))  
324                zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))            IF (iflag_pdf==0) THEN
325                zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)  
326                zpdf_e2(i) = sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))               DO i = 1, klon
327                zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))                  zdelq = min(ratqs(i, k), 0.99)*zq(i)
328                IF (zpdf_e1(i)<1.E-10) THEN                  rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
329                   rneb(i,k) = 0.                  zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
330                   zqn(i) = zqs(i)               END DO
331                ELSE  
332                   rneb(i,k) = 0.5*zpdf_e1(i)            ELSE
333                   zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)  
334                END IF               !   Version avec les nouvelles PDFs.
335                 DO i = 1, klon
336             END DO                  IF (zq(i)<1.E-15) THEN
337                       zq(i) = 1.E-15
338                    END IF
339          END IF               END DO
340          ! iflag_pdf                                                             DO i = 1, klon
341          DO i = 1, klon                  zpdf_sig(i) = ratqs(i, k)*zq(i)
342             IF (rneb(i,k)<=0.0) zqn(i) = 0.0                  zpdf_k(i) = -sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
343             IF (rneb(i,k)>=1.0) zqn(i) = zq(i)                  zpdf_delta(i) = log(zq(i)/zqs(i))
344             rneb(i,k) = max(0.0,min(1.0,rneb(i,k)))                  zpdf_a(i) = zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
345             !           zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))                  zpdf_b(i) = zpdf_k(i)/(2.*sqrt(2.))
346             !  On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par                  zpdf_e1(i) = zpdf_a(i) - zpdf_b(i)
347             !  la convection.                  zpdf_e1(i) = sign(min(abs(zpdf_e1(i)), 5.), zpdf_e1(i))
348             !  ATTENTION !!! Il va falloir verifier tout ca.                  zpdf_e1(i) = 1. - nr_erf(zpdf_e1(i))
349             zcond(i) = max(0.0,zqn(i)-zqs(i))*rneb(i,k)                  zpdf_e2(i) = zpdf_a(i) + zpdf_b(i)
350             !           print*,'ZDQS ',zdqs(i)                  zpdf_e2(i) = sign(min(abs(zpdf_e2(i)), 5.), zpdf_e2(i))
351             !--Olivier                  zpdf_e2(i) = 1. - nr_erf(zpdf_e2(i))
352             rhcl(i,k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)                  IF (zpdf_e1(i)<1.E-10) THEN
353             IF (rneb(i,k)<=0.0) rhcl(i,k) = zq(i)/zqs(i)                     rneb(i, k) = 0.
354             IF (rneb(i,k)>=1.0) rhcl(i,k) = 1.0                     zqn(i) = zqs(i)
355             !--fin                  ELSE
356          END DO                     rneb(i, k) = 0.5*zpdf_e1(i)
357       ELSE                     zqn(i) = zq(i)*zpdf_e2(i)/zpdf_e1(i)
358          DO i = 1, klon                  END IF
359             IF (zq(i)>zqs(i)) THEN  
360                rneb(i,k) = 1.0               END DO
361             ELSE  
362                rneb(i,k) = 0.0  
363             END IF            END IF
364             zcond(i) = max(0.0,zq(i)-zqs(i))/(1.+zdqs(i))            ! iflag_pdf                                              
365          END DO            DO i = 1, klon
366       END IF               IF (rneb(i, k)<=0.0) zqn(i) = 0.0
367                 IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
368       DO i = 1, klon               rneb(i, k) = max(0.0, min(1.0, rneb(i, k)))
369          zq(i) = zq(i) - zcond(i)               !  On ne divise pas par 1+zdqs pour forcer a avoir l'eau
370          !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD               !  predite par la convection.  ATTENTION !!! Il va
371          zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))               !  falloir verifier tout ca.
372       END DO               zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)
373                 !           print*, 'ZDQS ', zdqs(i)
374       ! Partager l'eau condensee en precipitation et eau liquide nuageuse               !--Olivier
375                 rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
376       DO i = 1, klon               IF (rneb(i, k)<=0.0) rhcl(i, k) = zq(i)/zqs(i)
377          IF (rneb(i,k)>0.0) THEN               IF (rneb(i, k)>=1.0) rhcl(i, k) = 1.0
378             zoliq(i) = zcond(i)               !--fin
379             zrho(i) = pplay(i,k)/zt(i)/rd            END DO
380             zdz(i) = (paprs(i,k)-paprs(i,k+1))/(zrho(i)*rg)         ELSE
381             zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)            DO i = 1, klon
382             zfice(i) = min(max(zfice(i),0.0),1.0)               IF (zq(i)>zqs(i)) THEN
383             zfice(i) = zfice(i)**nexpo                  rneb(i, k) = 1.0
384             zneb(i) = max(rneb(i,k),seuil_neb)               ELSE
385             radliq(i,k) = zoliq(i)/float(ninter+1)                  rneb(i, k) = 0.0
386          END IF               END IF
387       END DO               zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
388              END DO
389       DO n = 1, ninter         END IF
390          DO i = 1, klon  
391             IF (rneb(i,k)>0.0) THEN         DO i = 1, klon
392                zrhol(i) = zrho(i)*zoliq(i)/zneb(i)            zq(i) = zq(i) - zcond(i)
393              !         zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
394                IF (ptconv(i,k)) THEN            zt(i) = zt(i) + zcond(i)*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
395                   zcl(i) = cld_lc_con         END DO
396                   zct(i) = 1./cld_tau_con  
397                ELSE         ! Partager l'eau condensee en precipitation et eau liquide nuageuse
398                   zcl(i) = cld_lc_lsc  
399                   zct(i) = 1./cld_tau_lsc         DO i = 1, klon
400                END IF            IF (rneb(i, k)>0.0) THEN
401                !  quantité d'eau à élminier.               zoliq(i) = zcond(i)
402                zchau(i) = zct(i)*dtime/float(ninter)*zoliq(i)* &               zrho(i) = pplay(i, k)/zt(i)/rd
403                     (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))               zdz(i) = (paprs(i, k)-paprs(i, k+1))/(zrho(i)*rg)
404                !  meme chose pour la glace.               zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
405                IF (ptconv(i,k)) THEN               zfice(i) = min(max(zfice(i), 0.0), 1.0)
406                   zfroi(i) = dtime/float(ninter)/zdz(i)*zoliq(i)* &               zfice(i) = zfice(i)**nexpo
407                        fallvc(zrhol(i))*zfice(i)               zneb(i) = max(rneb(i, k), seuil_neb)
408                ELSE               radliq(i, k) = zoliq(i)/real(ninter+1)
409                   zfroi(i) = dtime/float(ninter)/zdz(i)*zoliq(i)* &            END IF
410                        fallvs(zrhol(i))*zfice(i)         END DO
411                END IF  
412                ztot(i) = zchau(i) + zfroi(i)         DO n = 1, ninter
413                IF (zneb(i)==seuil_neb) ztot(i) = 0.0            DO i = 1, klon
414                ztot(i) = min(max(ztot(i),0.0),zoliq(i))               IF (rneb(i, k)>0.0) THEN
415                zoliq(i) = max(zoliq(i)-ztot(i),0.0)                  zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
416                radliq(i,k) = radliq(i,k) + zoliq(i)/float(ninter+1)  
417             END IF                  IF (ptconv(i, k)) THEN
418          END DO                     zcl(i) = cld_lc_con
419       END DO                     zct(i) = 1./cld_tau_con
420                    ELSE
421       DO i = 1, klon                     zcl(i) = cld_lc_lsc
422          IF (rneb(i,k)>0.0) THEN                     zct(i) = 1./cld_tau_lsc
423             d_ql(i,k) = zoliq(i)                  END IF
424             zrfl(i) = zrfl(i) + max(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i &                  !  quantité d'eau à élminier.
425                  ,k+1))/(rg*dtime)                  zchau(i) = zct(i)*dtime/real(ninter)*zoliq(i)* &
426          END IF                       (1.0-exp(-(zoliq(i)/zneb(i)/zcl(i))**2))*(1.-zfice(i))
427          IF (zt(i)<rtt) THEN                  !  meme chose pour la glace.
428             psfl(i,k) = zrfl(i)                  IF (ptconv(i, k)) THEN
429          ELSE                     zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &
430             prfl(i,k) = zrfl(i)                          fallvc(zrhol(i))*zfice(i)
431          END IF                  ELSE
432       END DO                     zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)* &
433                            fallvs(zrhol(i))*zfice(i)
434       ! Calculer les tendances de q et de t:                  END IF
435                    ztot(i) = zchau(i) + zfroi(i)
436       DO i = 1, klon                  IF (zneb(i)==seuil_neb) ztot(i) = 0.0
437          d_q(i,k) = zq(i) - q(i,k)                  ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
438          d_t(i,k) = zt(i) - t(i,k)                  zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
439       END DO                  radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
440                 END IF
441       !AA--------------- Calcul du lessivage stratiforme  -------------            END DO
442           END DO
443       DO i = 1, klon  
444          zprec_cond(i) = max(zcond(i)-zoliq(i),0.0)* &         DO i = 1, klon
445               (paprs(i,k)-paprs(i,k+1))/rg            IF (rneb(i, k)>0.0) THEN
446          IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN               d_ql(i, k) = zoliq(i)
447             !AA lessivage nucleation LMD5 dans la couche elle-meme               zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.) &
448             IF (t(i,k)>=ztglace) THEN                    * (paprs(i, k) - paprs(i, k + 1)) / (rg * dtime)
449                zalpha_tr = a_tr_sca(3)            END IF
450             ELSE            IF (zt(i)<rtt) THEN
451                zalpha_tr = a_tr_sca(4)               psfl(i, k) = zrfl(i)
452             END IF            ELSE
453             zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))               prfl(i, k) = zrfl(i)
454             pfrac_nucl(i,k) = pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi)            END IF
455             frac_nucl(i,k) = 1. - zneb(i)*zfrac_lessi         END DO
456    
457             ! nucleation avec un facteur -1 au lieu de -0.5         ! Calculer les tendances de q et de t:
458             zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))  
459             pfrac_1nucl(i,k) = pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi)         DO i = 1, klon
460          END IF            d_q(i, k) = zq(i) - q(i, k)
461              d_t(i, k) = zt(i) - t(i, k)
462           END DO
463       END DO  
464       !AA Lessivage par impaction dans les couches en-dessous         !AA--------------- Calcul du lessivage stratiforme  -------------
465       ! boucle sur i                                          
466       DO kk = k - 1, 1, -1         DO i = 1, klon
467          DO i = 1, klon            zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)* &
468             IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN                 (paprs(i, k)-paprs(i, k+1))/rg
469                IF (t(i,kk)>=ztglace) THEN            IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
470                   zalpha_tr = a_tr_sca(1)               !AA lessivage nucleation LMD5 dans la couche elle-meme
471                ELSE               IF (t(i, k)>=ztglace) THEN
472                   zalpha_tr = a_tr_sca(2)                  zalpha_tr = a_tr_sca(3)
473                END IF               ELSE
474                zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))                  zalpha_tr = a_tr_sca(4)
475                pfrac_impa(i,kk) = pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi)               END IF
476                frac_impa(i,kk) = 1. - zneb(i)*zfrac_lessi               zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
477             END IF               pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
478          END DO               frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
479       END DO  
480                 ! nucleation avec un facteur -1 au lieu de -0.5
481       !AA----------------------------------------------------------               zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
482       !                     FIN DE BOUCLE SUR K               pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
483    end DO            END IF
484    
485    !AA-----------------------------------------------------------  
486           END DO
487    ! Pluie ou neige au sol selon la temperature de la 1ere couche         !AA Lessivage par impaction dans les couches en-dessous
488           ! boucle sur i                                        
489    DO i = 1, klon         DO kk = k - 1, 1, -1
490       IF ((t(i,1)+d_t(i,1))<rtt) THEN            DO i = 1, klon
491          snow(i) = zrfl(i)               IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
492          zlh_solid(i) = rlstt - rlvtt                  IF (t(i, kk)>=ztglace) THEN
493       ELSE                     zalpha_tr = a_tr_sca(1)
494          rain(i) = zrfl(i)                  ELSE
495          zlh_solid(i) = 0.                     zalpha_tr = a_tr_sca(2)
496       END IF                  END IF
497    END DO                  zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
498                    pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
499    ! For energy conservation : when snow is present, the solification                  frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
500    ! latent heat is considered.               END IF
501    DO k = 1, klev            END DO
502       DO i = 1, klon         END DO
503          zcpair = rcpd*(1.0+rvtmp2*(q(i,k)+d_q(i,k)))  
504          zmair = (paprs(i,k)-paprs(i,k+1))/rg         !AA----------------------------------------------------------
505          zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime         !                     FIN DE BOUCLE SUR K
506          d_t(i,k) = d_t(i,k) + zlh_solid(i)*zm_solid/(zcpair*zmair)      end DO
507       END DO  
508    END DO      !AA-----------------------------------------------------------
509    
510        ! Pluie ou neige au sol selon la temperature de la 1ere couche
511    
512        DO i = 1, klon
513           IF ((t(i, 1)+d_t(i, 1))<rtt) THEN
514              snow(i) = zrfl(i)
515              zlh_solid(i) = rlstt - rlvtt
516           ELSE
517              rain(i) = zrfl(i)
518              zlh_solid(i) = 0.
519           END IF
520        END DO
521    
522        ! For energy conservation: when snow is present, the solification
523        ! latent heat is considered.
524        DO k = 1, klev
525           DO i = 1, klon
526              zcpair = rcpd*(1.0+rvtmp2*(q(i, k)+d_q(i, k)))
527              zmair = (paprs(i, k)-paprs(i, k+1))/rg
528              zm_solid = (prfl(i, k)-prfl(i, k+1)+psfl(i, k)-psfl(i, k+1))*dtime
529              d_t(i, k) = d_t(i, k) + zlh_solid(i)*zm_solid/(zcpair*zmair)
530           END DO
531        END DO
532    
533  END SUBROUTINE fisrtilp    END SUBROUTINE fisrtilp
534    
535    end module fisrtilp_m

Legend:
Removed from v.52  
changed lines
  Added in v.76

  ViewVC Help
Powered by ViewVC 1.1.21