source: LMDZ6/trunk/libf/phylmd/fisrtilp_tr.f90 @ 5279

Last change on this file since 5279 was 5274, checked in by abarral, 9 months ago

Replace yomcst.h by existing module

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.8 KB
Line 
1
2! $Id: fisrtilp_tr.f90 5274 2024-10-25 13:41:23Z abarral $
3
4
5SUBROUTINE fisrtilp_tr(dtime, paprs, pplay, t, q, ratqs, d_t, d_q, d_ql, &
6    rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
7    frac_nucl, prfl, psfl, rhcl) ! relative humidity in clear sky (needed for aer optical
8  ! properties; aeropt.F)
9
10
11  USE dimphy
12  USE print_control_mod, ONLY: lunout
13  USE yomcst_mod_h, ONLY: RPI, RCLUM, RHPLA, RKBOL, RNAVO                   &
14          , RDAY, REA, REPSM, RSIYEA, RSIDAY, ROMEGA                  &
15          , R_ecc, R_peri, R_incl                                      &
16          , RA, RG, R1SA                                         &
17          , RSIGMA                                                     &
18          , R, RMD, RMV, RD, RV, RCPD                    &
19          , RMO3, RMCO2, RMC, RMCH4, RMN2O, RMCFC11, RMCFC12        &
20          , RCPV, RCVD, RCVV, RKAPPA, RETV, eps_w                    &
21          , RCW, RCS                                                 &
22          , RLVTT, RLSTT, RLMLT, RTT, RATM                           &
23          , RESTT, RALPW, RBETW, RGAMW, RALPS, RBETS, RGAMS            &
24          , RALPD, RBETD, RGAMD
25IMPLICIT NONE
26  ! ======================================================================
27  ! Auteur(s): Z.X. Li (LMD/CNRS)
28  ! Date: le 20 mars 1995
29  ! Objet: condensation et precipitation stratiforme.
30  ! schema de nuage
31  ! ======================================================================
32  ! ======================================================================
33
34
35  ! Arguments:
36
37  REAL dtime ! intervalle du temps (s)
38  REAL paprs(klon, klev+1) ! pression a inter-couche
39  REAL pplay(klon, klev) ! pression au milieu de couche
40  REAL t(klon, klev) ! temperature (K)
41  REAL q(klon, klev) ! humidite specifique (kg/kg)
42  REAL d_t(klon, klev) ! incrementation de la temperature (K)
43  REAL d_q(klon, klev) ! incrementation de la vapeur d'eau
44  REAL d_ql(klon, klev) ! incrementation de l'eau liquide
45  REAL rneb(klon, klev) ! fraction nuageuse
46  REAL radliq(klon, klev) ! eau liquide utilisee dans rayonnements
47  REAL rain(klon) ! pluies (mm/s)
48  REAL snow(klon) ! neige (mm/s)
49  REAL prfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
50  REAL psfl(klon, klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
51
52  ! jq   For aerosol opt properties needed (see aeropt.F)
53  REAL rhcl(klon, klev)
54
55  ! AA
56  ! Coeffients de fraction lessivee : pour OFF-LINE
57
58  REAL pfrac_nucl(klon, klev)
59  REAL pfrac_1nucl(klon, klev)
60  REAL pfrac_impa(klon, klev)
61
62  ! Fraction d'aerosols lessivee par impaction et par nucleation
63  ! POur ON-LINE
64
65  REAL frac_impa(klon, klev)
66  REAL frac_nucl(klon, klev)
67  ! AA
68
69  ! Options du programme:
70
71  REAL seuil_neb ! un nuage existe vraiment au-dela
72  PARAMETER (seuil_neb=0.001)
73  REAL ct ! inverse du temps pour qu'un nuage precipite
74  PARAMETER (ct=1./1800.)
75  REAL cl ! seuil de precipitation
76  PARAMETER (cl=2.6E-4)
77  ! cc      PARAMETER (cl=2.3e-4)
78  ! cc      PARAMETER (cl=2.0e-4)
79  INTEGER ninter ! sous-intervals pour la precipitation
80  PARAMETER (ninter=5)
81  LOGICAL evap_prec ! evaporation de la pluie
82  PARAMETER (evap_prec=.TRUE.)
83  REAL coef_eva
84  PARAMETER (coef_eva=2.0E-05)
85  LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
86  REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur
87  PARAMETER (calcrat=.TRUE.)
88  REAL zx_min, rat_max
89  PARAMETER (zx_min=1.0, rat_max=0.01)
90  REAL zx_max, rat_min
91  PARAMETER (zx_max=0.1, rat_min=0.3)
92  REAL zx
93
94  LOGICAL cpartiel ! condensation partielle
95  PARAMETER (cpartiel=.TRUE.)
96  REAL t_coup
97  PARAMETER (t_coup=234.0)
98
99  ! Variables locales:
100
101  INTEGER i, k, n, kk
102  REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
103  REAL zrfl(klon), zrfln(klon), zqev, zqevt
104  REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
105  REAL ztglace, zt(klon)
106  INTEGER nexpo ! exponentiel pour glace/eau
107  REAL zdz(klon), zrho(klon), ztot(klon), zrhol(klon)
108  REAL zchau(klon), zfroi(klon), zfice(klon), zneb(klon)
109
110  LOGICAL appel1er
111  SAVE appel1er
112  !$OMP THREADPRIVATE(appel1er)
113
114  ! ---------------------------------------------------------------
115
116  ! AA Variables traceurs:
117  ! AA  Provisoire !!! Parametres alpha du lessivage
118  ! AA  A priori on a 4 scavenging # possibles
119
120  REAL a_tr_sca(4)
121  SAVE a_tr_sca
122  !$OMP THREADPRIVATE(a_tr_sca)
123
124  ! Variables intermediaires
125
126  REAL zalpha_tr
127  REAL zfrac_lessi
128  REAL zprec_cond(klon)
129  ! AA
130  ! ---------------------------------------------------------------
131
132  ! Fonctions en ligne:
133
134  REAL fallv ! vitesse de chute pour crystaux de glace
135  REAL zzz
136  include "YOETHF.h"
137  include "FCTTRE.h"
138  fallv(zzz) = 3.29/2.0*((zzz)**0.16)
139  ! cc      fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
140  ! cc      fallv (zzz) = 3.29 * ((zzz)**0.16)
141
142  DATA appel1er/.TRUE./
143
144  IF (appel1er) THEN
145
146    WRITE (lunout, *) 'fisrtilp, calcrat:', calcrat
147    WRITE (lunout, *) 'fisrtilp, ninter:', ninter
148    WRITE (lunout, *) 'fisrtilp, evap_prec:', evap_prec
149    WRITE (lunout, *) 'fisrtilp, cpartiel:', cpartiel
150    IF (abs(dtime/real(ninter)-360.0)>0.001) THEN
151      WRITE (lunout, *) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
152      WRITE (lunout, *) 'Je prefere un sous-intervalle de 6 minutes'
153      CALL abort
154    END IF
155    appel1er = .FALSE.
156
157    ! AA initialiation provisoire
158    a_tr_sca(1) = -0.5
159    a_tr_sca(2) = -0.5
160    a_tr_sca(3) = -0.5
161    a_tr_sca(4) = -0.5
162
163    ! AA Initialisation a 1 des coefs des fractions lessivees
164
165    DO k = 1, klev
166      DO i = 1, klon
167        pfrac_nucl(i, k) = 1.
168        pfrac_1nucl(i, k) = 1.
169        pfrac_impa(i, k) = 1.
170      END DO
171    END DO
172
173  END IF !  test sur appel1er
174
175  ! MAf Initialisation a 0 de zoliq
176  DO i = 1, klon
177    zoliq(i) = 0.
178  END DO
179  ! Determiner les nuages froids par leur temperature
180
181  ztglace = rtt - 15.0
182  nexpo = 6
183  ! cc      nexpo = 1
184
185  ! Initialiser les sorties:
186
187  DO k = 1, klev + 1
188    DO i = 1, klon
189      prfl(i, k) = 0.0
190      psfl(i, k) = 0.0
191    END DO
192  END DO
193
194  DO k = 1, klev
195    DO i = 1, klon
196      d_t(i, k) = 0.0
197      d_q(i, k) = 0.0
198      d_ql(i, k) = 0.0
199      rneb(i, k) = 0.0
200      radliq(i, k) = 0.0
201      frac_nucl(i, k) = 1.
202      frac_impa(i, k) = 1.
203    END DO
204  END DO
205  DO i = 1, klon
206    rain(i) = 0.0
207    snow(i) = 0.0
208  END DO
209
210  ! Initialiser le flux de precipitation a zero
211
212  DO i = 1, klon
213    zrfl(i) = 0.0
214    zneb(i) = seuil_neb
215  END DO
216
217
218  ! AA Pour plus de securite
219
220  zalpha_tr = 0.
221  zfrac_lessi = 0.
222
223  ! AA----------------------------------------------------------
224
225  ! Boucle verticale (du haut vers le bas)
226
227  DO k = klev, 1, -1
228
229    ! AA----------------------------------------------------------
230
231    DO i = 1, klon
232      zt(i) = t(i, k)
233      zq(i) = q(i, k)
234    END DO
235
236    ! Calculer l'evaporation de la precipitation
237
238    IF (evap_prec) THEN
239      DO i = 1, klon
240        IF (zrfl(i)>0.) THEN
241          IF (thermcep) THEN
242            zdelta = max(0., sign(1.,rtt-zt(i)))
243            zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
244            zqs(i) = min(0.5, zqs(i))
245            zcor = 1./(1.-retv*zqs(i))
246            zqs(i) = zqs(i)*zcor
247          ELSE
248            IF (zt(i)<t_coup) THEN
249              zqs(i) = qsats(zt(i))/pplay(i, k)
250            ELSE
251              zqs(i) = qsatl(zt(i))/pplay(i, k)
252            END IF
253          END IF
254          zqev = max(0.0, (zqs(i)-zq(i))*zneb(i))
255          zqevt = coef_eva*(1.0-zq(i)/zqs(i))*sqrt(zrfl(i))* &
256            (paprs(i,k)-paprs(i,k+1))/pplay(i, k)*zt(i)*rd/rg
257          zqevt = max(0.0, min(zqevt,zrfl(i)))*rg*dtime/ &
258            (paprs(i,k)-paprs(i,k+1))
259          zqev = min(zqev, zqevt)
260          zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))/rg/dtime
261          zq(i) = zq(i) - (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
262            k+1)))*dtime
263          zt(i) = zt(i) + (zrfln(i)-zrfl(i))*(rg/(paprs(i,k)-paprs(i, &
264            k+1)))*dtime*rlvtt/rcpd/(1.0+rvtmp2*zq(i))
265          zrfl(i) = zrfln(i)
266        END IF
267      END DO
268    END IF
269
270    ! Calculer Qs et L/Cp*dQs/dT:
271
272    IF (thermcep) THEN
273      DO i = 1, klon
274        zdelta = max(0., sign(1.,rtt-zt(i)))
275        zcvm5 = r5les*rlvtt*(1.-zdelta) + r5ies*rlstt*zdelta
276        zcvm5 = zcvm5/rcpd/(1.0+rvtmp2*zq(i))
277        zqs(i) = r2es*foeew(zt(i), zdelta)/pplay(i, k)
278        zqs(i) = min(0.5, zqs(i))
279        zcor = 1./(1.-retv*zqs(i))
280        zqs(i) = zqs(i)*zcor
281        zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
282      END DO
283    ELSE
284      DO i = 1, klon
285        IF (zt(i)<t_coup) THEN
286          zqs(i) = qsats(zt(i))/pplay(i, k)
287          zdqs(i) = dqsats(zt(i), zqs(i))
288        ELSE
289          zqs(i) = qsatl(zt(i))/pplay(i, k)
290          zdqs(i) = dqsatl(zt(i), zqs(i))
291        END IF
292      END DO
293    END IF
294
295    ! Determiner la condensation partielle et calculer la quantite
296    ! de l'eau condensee:
297
298    IF (cpartiel) THEN
299      DO i = 1, klon
300
301        zdelq = ratqs(i, k)*zq(i)
302        rneb(i, k) = (zq(i)+zdelq-zqs(i))/(2.0*zdelq)
303        zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
304        IF (rneb(i,k)<=0.0) zqn(i) = 0.0
305        IF (rneb(i,k)>=1.0) zqn(i) = zq(i)
306        rneb(i, k) = max(0.0, min(1.0,rneb(i,k)))
307        zcond(i) = max(0.0, zqn(i)-zqs(i))*rneb(i, k)/(1.+zdqs(i))
308
309        ! --Olivier
310        rhcl(i, k) = (zqs(i)+zq(i)-zdelq)/2./zqs(i)
311        IF (rneb(i,k)<=0.0) rhcl(i, k) = zq(i)/zqs(i)
312        IF (rneb(i,k)>=1.0) rhcl(i, k) = 1.0
313        ! --fin
314
315      END DO
316    ELSE
317      DO i = 1, klon
318        IF (zq(i)>zqs(i)) THEN
319          rneb(i, k) = 1.0
320        ELSE
321          rneb(i, k) = 0.0
322        END IF
323        zcond(i) = max(0.0, zq(i)-zqs(i))/(1.+zdqs(i))
324      END DO
325    END IF
326
327    DO i = 1, klon
328      zq(i) = zq(i) - zcond(i)
329      zt(i) = zt(i) + zcond(i)*rlvtt/rcpd
330    END DO
331
332    ! Partager l'eau condensee en precipitation et eau liquide nuageuse
333
334    DO i = 1, klon
335      IF (rneb(i,k)>0.0) THEN
336        zoliq(i) = zcond(i)
337        zrho(i) = pplay(i, k)/zt(i)/rd
338        zdz(i) = (paprs(i,k)-paprs(i,k+1))/(zrho(i)*rg)
339        zfice(i) = 1.0 - (zt(i)-ztglace)/(273.13-ztglace)
340        zfice(i) = min(max(zfice(i),0.0), 1.0)
341        zfice(i) = zfice(i)**nexpo
342        zneb(i) = max(rneb(i,k), seuil_neb)
343        radliq(i, k) = zoliq(i)/real(ninter+1)
344      END IF
345    END DO
346
347    DO n = 1, ninter
348      DO i = 1, klon
349        IF (rneb(i,k)>0.0) THEN
350          zchau(i) = ct*dtime/real(ninter)*zoliq(i)* &
351            (1.0-exp(-(zoliq(i)/zneb(i)/cl)**2))*(1.-zfice(i))
352          zrhol(i) = zrho(i)*zoliq(i)/zneb(i)
353          zfroi(i) = dtime/real(ninter)/zdz(i)*zoliq(i)*fallv(zrhol(i))* &
354            zfice(i)
355          ztot(i) = zchau(i) + zfroi(i)
356          IF (zneb(i)==seuil_neb) ztot(i) = 0.0
357          ztot(i) = min(max(ztot(i),0.0), zoliq(i))
358          zoliq(i) = max(zoliq(i)-ztot(i), 0.0)
359          radliq(i, k) = radliq(i, k) + zoliq(i)/real(ninter+1)
360        END IF
361      END DO
362    END DO
363
364    DO i = 1, klon
365      IF (rneb(i,k)>0.0) THEN
366        d_ql(i, k) = zoliq(i)
367        zrfl(i) = zrfl(i) + max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k &
368          +1))/(rg*dtime)
369      END IF
370      IF (zt(i)<rtt) THEN
371        psfl(i, k) = zrfl(i)
372      ELSE
373        prfl(i, k) = zrfl(i)
374      END IF
375    END DO
376
377    ! Calculer les tendances de q et de t:
378
379    DO i = 1, klon
380      d_q(i, k) = zq(i) - q(i, k)
381      d_t(i, k) = zt(i) - t(i, k)
382    END DO
383
384    ! AA--------------- Calcul du lessivage stratiforme  -------------
385
386    DO i = 1, klon
387
388      zprec_cond(i) = max(zcond(i)-zoliq(i), 0.0)*(paprs(i,k)-paprs(i,k+1))/ &
389        rg
390      IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN
391        ! AA lessivage nucleation LMD5 dans la couche elle-meme
392        IF (t(i,k)>=ztglace) THEN
393          zalpha_tr = a_tr_sca(3)
394        ELSE
395          zalpha_tr = a_tr_sca(4)
396        END IF
397        zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
398        pfrac_nucl(i, k) = pfrac_nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
399        frac_nucl(i, k) = 1. - zneb(i)*zfrac_lessi
400
401        ! nucleation avec un facteur -1 au lieu de -0.5
402        zfrac_lessi = 1. - exp(-zprec_cond(i)/zneb(i))
403        pfrac_1nucl(i, k) = pfrac_1nucl(i, k)*(1.-zneb(i)*zfrac_lessi)
404      END IF
405
406    END DO ! boucle sur i
407
408    ! AA Lessivage par impaction dans les couches en-dessous
409    DO kk = k - 1, 1, -1
410      DO i = 1, klon
411        IF (rneb(i,k)>0.0 .AND. zprec_cond(i)>0.) THEN
412          IF (t(i,kk)>=ztglace) THEN
413            zalpha_tr = a_tr_sca(1)
414          ELSE
415            zalpha_tr = a_tr_sca(2)
416          END IF
417          zfrac_lessi = 1. - exp(zalpha_tr*zprec_cond(i)/zneb(i))
418          pfrac_impa(i, kk) = pfrac_impa(i, kk)*(1.-zneb(i)*zfrac_lessi)
419          frac_impa(i, kk) = 1. - zneb(i)*zfrac_lessi
420        END IF
421      END DO
422    END DO
423
424    ! AA----------------------------------------------------------
425    ! FIN DE BOUCLE SUR K
426  END DO
427
428  ! AA-----------------------------------------------------------
429
430  ! Pluie ou neige au sol selon la temperature de la 1ere couche
431
432  DO i = 1, klon
433    IF ((t(i,1)+d_t(i,1))<rtt) THEN
434      snow(i) = zrfl(i)
435    ELSE
436      rain(i) = zrfl(i)
437    END IF
438  END DO
439
440  RETURN
441END SUBROUTINE fisrtilp_tr
Note: See TracBrowser for help on using the repository browser.