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

Last change on this file since 5353 was 5285, checked in by abarral, 6 weeks ago

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