source: LMDZ6/branches/Amaury_dev/libf/phylmd/fisrtilp_tr.F90 @ 5218

Last change on this file since 5218 was 5153, checked in by abarral, 5 months ago

Revert FCTTRE to INCLUDE to assess impact of inlining

  • 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.4 KB
Line 
1! $Id: fisrtilp_tr.F90 5153 2024-07-31 16:20:03Z abarral $
2
3
4SUBROUTINE fisrtilp_tr(dtime, paprs, pplay, t, q, ratqs, d_t, d_q, d_ql, &
5        rneb, radliq, rain, snow, pfrac_impa, pfrac_nucl, pfrac_1nucl, frac_impa, &
6        frac_nucl, prfl, psfl, rhcl) ! relative humidity in clear sky (needed for aer optical
7  ! properties; aeropt.F)
8
9  USE dimphy
10  USE lmdz_print_control, ONLY: lunout
11  USE lmdz_yoethf
12
13  USE lmdz_yomcst
14
15  IMPLICIT NONE
16 INCLUDE "FCTTRE.h"
17  ! ======================================================================
18  ! Auteur(s): Z.X. Li (LMD/CNRS)
19  ! Date: le 20 mars 1995
20  ! Objet: condensation et precipitation stratiforme.
21  ! schema de nuage
22  ! ======================================================================
23  ! ======================================================================
24
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  DATA appel1er/.TRUE./
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  IF (appel1er) THEN
132
133    WRITE (lunout, *) 'fisrtilp, calcrat:', calcrat
134    WRITE (lunout, *) 'fisrtilp, ninter:', ninter
135    WRITE (lunout, *) 'fisrtilp, evap_prec:', evap_prec
136    WRITE (lunout, *) 'fisrtilp, cpartiel:', cpartiel
137    IF (abs(dtime / real(ninter) - 360.0)>0.001) THEN
138      WRITE (lunout, *) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
139      WRITE (lunout, *) 'Je prefere un sous-intervalle de 6 minutes'
140      CALL abort
141    END IF
142    appel1er = .FALSE.
143
144    ! AA initialiation provisoire
145    a_tr_sca(1) = -0.5
146    a_tr_sca(2) = -0.5
147    a_tr_sca(3) = -0.5
148    a_tr_sca(4) = -0.5
149
150    ! AA Initialisation a 1 des coefs des fractions lessivees
151
152    DO k = 1, klev
153      DO i = 1, klon
154        pfrac_nucl(i, k) = 1.
155        pfrac_1nucl(i, k) = 1.
156        pfrac_impa(i, k) = 1.
157      END DO
158    END DO
159
160  END IF !  test sur appel1er
161
162  ! MAf Initialisation a 0 de zoliq
163  DO i = 1, klon
164    zoliq(i) = 0.
165  END DO
166  ! Determiner les nuages froids par leur temperature
167
168  ztglace = rtt - 15.0
169  nexpo = 6
170  ! cc      nexpo = 1
171
172  ! Initialiser les sorties:
173
174  DO k = 1, klev + 1
175    DO i = 1, klon
176      prfl(i, k) = 0.0
177      psfl(i, k) = 0.0
178    END DO
179  END DO
180
181  DO k = 1, klev
182    DO i = 1, klon
183      d_t(i, k) = 0.0
184      d_q(i, k) = 0.0
185      d_ql(i, k) = 0.0
186      rneb(i, k) = 0.0
187      radliq(i, k) = 0.0
188      frac_nucl(i, k) = 1.
189      frac_impa(i, k) = 1.
190    END DO
191  END DO
192  DO i = 1, klon
193    rain(i) = 0.0
194    snow(i) = 0.0
195  END DO
196
197  ! Initialiser le flux de precipitation a zero
198
199  DO i = 1, klon
200    zrfl(i) = 0.0
201    zneb(i) = seuil_neb
202  END DO
203
204
205  ! AA Pour plus de securite
206
207  zalpha_tr = 0.
208  zfrac_lessi = 0.
209
210  ! AA----------------------------------------------------------
211
212  ! Boucle verticale (du haut vers le bas)
213
214  DO k = klev, 1, -1
215
216    ! AA----------------------------------------------------------
217
218    DO i = 1, klon
219      zt(i) = t(i, k)
220      zq(i) = q(i, k)
221    END DO
222
223    ! Calculer l'evaporation de la precipitation
224
225    IF (evap_prec) THEN
226      DO i = 1, klon
227        IF (zrfl(i)>0.) THEN
228          IF (thermcep) THEN
229            zdelta = max(0., sign(1., rtt - zt(i)))
230            zqs(i) = r2es * foeew(zt(i), zdelta) / pplay(i, k)
231            zqs(i) = min(0.5, zqs(i))
232            zcor = 1. / (1. - retv * zqs(i))
233            zqs(i) = zqs(i) * zcor
234          ELSE
235            IF (zt(i)<t_coup) THEN
236              zqs(i) = qsats(zt(i)) / pplay(i, k)
237            ELSE
238              zqs(i) = qsatl(zt(i)) / pplay(i, k)
239            END IF
240          END IF
241          zqev = max(0.0, (zqs(i) - zq(i)) * zneb(i))
242          zqevt = coef_eva * (1.0 - zq(i) / zqs(i)) * sqrt(zrfl(i)) * &
243                  (paprs(i, k) - paprs(i, k + 1)) / pplay(i, k) * zt(i) * rd / rg
244          zqevt = max(0.0, min(zqevt, zrfl(i))) * rg * dtime / &
245                  (paprs(i, k) - paprs(i, k + 1))
246          zqev = min(zqev, zqevt)
247          zrfln(i) = zrfl(i) - zqev * (paprs(i, k) - paprs(i, k + 1)) / rg / dtime
248          zq(i) = zq(i) - (zrfln(i) - zrfl(i)) * (rg / (paprs(i, k) - paprs(i, &
249                  k + 1))) * dtime
250          zt(i) = zt(i) + (zrfln(i) - zrfl(i)) * (rg / (paprs(i, k) - paprs(i, &
251                  k + 1))) * dtime * rlvtt / rcpd / (1.0 + rvtmp2 * zq(i))
252          zrfl(i) = zrfln(i)
253        END IF
254      END DO
255    END IF
256
257    ! Calculer Qs et L/Cp*dQs/dT:
258
259    IF (thermcep) THEN
260      DO i = 1, klon
261        zdelta = max(0., sign(1., rtt - zt(i)))
262        zcvm5 = r5les * rlvtt * (1. - zdelta) + r5ies * rlstt * zdelta
263        zcvm5 = zcvm5 / rcpd / (1.0 + rvtmp2 * zq(i))
264        zqs(i) = r2es * foeew(zt(i), zdelta) / pplay(i, k)
265        zqs(i) = min(0.5, zqs(i))
266        zcor = 1. / (1. - retv * zqs(i))
267        zqs(i) = zqs(i) * zcor
268        zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
269      END DO
270    ELSE
271      DO i = 1, klon
272        IF (zt(i)<t_coup) THEN
273          zqs(i) = qsats(zt(i)) / pplay(i, k)
274          zdqs(i) = dqsats(zt(i), zqs(i))
275        ELSE
276          zqs(i) = qsatl(zt(i)) / pplay(i, k)
277          zdqs(i) = dqsatl(zt(i), zqs(i))
278        END IF
279      END DO
280    END IF
281
282    ! Determiner la condensation partielle et calculer la quantite
283    ! de l'eau condensee:
284
285    IF (cpartiel) THEN
286      DO i = 1, klon
287
288        zdelq = ratqs(i, k) * zq(i)
289        rneb(i, k) = (zq(i) + zdelq - zqs(i)) / (2.0 * zdelq)
290        zqn(i) = (zq(i) + zdelq + zqs(i)) / 2.0
291        IF (rneb(i, k)<=0.0) zqn(i) = 0.0
292        IF (rneb(i, k)>=1.0) zqn(i) = zq(i)
293        rneb(i, k) = max(0.0, min(1.0, rneb(i, k)))
294        zcond(i) = max(0.0, zqn(i) - zqs(i)) * rneb(i, k) / (1. + zdqs(i))
295
296        ! --Olivier
297        rhcl(i, k) = (zqs(i) + zq(i) - zdelq) / 2. / zqs(i)
298        IF (rneb(i, k)<=0.0) rhcl(i, k) = zq(i) / zqs(i)
299        IF (rneb(i, k)>=1.0) rhcl(i, k) = 1.0
300        ! --fin
301
302      END DO
303    ELSE
304      DO i = 1, klon
305        IF (zq(i)>zqs(i)) THEN
306          rneb(i, k) = 1.0
307        ELSE
308          rneb(i, k) = 0.0
309        END IF
310        zcond(i) = max(0.0, zq(i) - zqs(i)) / (1. + zdqs(i))
311      END DO
312    END IF
313
314    DO i = 1, klon
315      zq(i) = zq(i) - zcond(i)
316      zt(i) = zt(i) + zcond(i) * rlvtt / rcpd
317    END DO
318
319    ! Partager l'eau condensee en precipitation et eau liquide nuageuse
320
321    DO i = 1, klon
322      IF (rneb(i, k)>0.0) THEN
323        zoliq(i) = zcond(i)
324        zrho(i) = pplay(i, k) / zt(i) / rd
325        zdz(i) = (paprs(i, k) - paprs(i, k + 1)) / (zrho(i) * rg)
326        zfice(i) = 1.0 - (zt(i) - ztglace) / (273.13 - ztglace)
327        zfice(i) = min(max(zfice(i), 0.0), 1.0)
328        zfice(i) = zfice(i)**nexpo
329        zneb(i) = max(rneb(i, k), seuil_neb)
330        radliq(i, k) = zoliq(i) / real(ninter + 1)
331      END IF
332    END DO
333
334    DO n = 1, ninter
335      DO i = 1, klon
336        IF (rneb(i, k)>0.0) THEN
337          zchau(i) = ct * dtime / real(ninter) * zoliq(i) * &
338                  (1.0 - exp(-(zoliq(i) / zneb(i) / cl)**2)) * (1. - zfice(i))
339          zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
340          zfroi(i) = dtime / real(ninter) / zdz(i) * zoliq(i) * fallv(zrhol(i)) * &
341                  zfice(i)
342          ztot(i) = zchau(i) + zfroi(i)
343          IF (zneb(i)==seuil_neb) ztot(i) = 0.0
344          ztot(i) = min(max(ztot(i), 0.0), zoliq(i))
345          zoliq(i) = max(zoliq(i) - ztot(i), 0.0)
346          radliq(i, k) = radliq(i, k) + zoliq(i) / real(ninter + 1)
347        END IF
348      END DO
349    END DO
350
351    DO i = 1, klon
352      IF (rneb(i, k)>0.0) THEN
353        d_ql(i, k) = zoliq(i)
354        zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.0) * (paprs(i, k) - paprs(i, k &
355                + 1)) / (rg * dtime)
356      END IF
357      IF (zt(i)<rtt) THEN
358        psfl(i, k) = zrfl(i)
359      ELSE
360        prfl(i, k) = zrfl(i)
361      END IF
362    END DO
363
364    ! Calculer les tendances de q et de t:
365
366    DO i = 1, klon
367      d_q(i, k) = zq(i) - q(i, k)
368      d_t(i, k) = zt(i) - t(i, k)
369    END DO
370
371    ! AA--------------- Calcul du lessivage stratiforme  -------------
372
373    DO i = 1, klon
374
375      zprec_cond(i) = max(zcond(i) - zoliq(i), 0.0) * (paprs(i, k) - paprs(i, k + 1)) / &
376              rg
377      IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
378        ! AA lessivage nucleation LMD5 dans la couche elle-meme
379        IF (t(i, k)>=ztglace) THEN
380          zalpha_tr = a_tr_sca(3)
381        ELSE
382          zalpha_tr = a_tr_sca(4)
383        END IF
384        zfrac_lessi = 1. - exp(zalpha_tr * zprec_cond(i) / zneb(i))
385        pfrac_nucl(i, k) = pfrac_nucl(i, k) * (1. - zneb(i) * zfrac_lessi)
386        frac_nucl(i, k) = 1. - zneb(i) * zfrac_lessi
387
388        ! nucleation avec un facteur -1 au lieu de -0.5
389        zfrac_lessi = 1. - exp(-zprec_cond(i) / zneb(i))
390        pfrac_1nucl(i, k) = pfrac_1nucl(i, k) * (1. - zneb(i) * zfrac_lessi)
391      END IF
392
393    END DO ! boucle sur i
394
395    ! AA Lessivage par impaction dans les couches en-dessous
396    DO kk = k - 1, 1, -1
397      DO i = 1, klon
398        IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
399          IF (t(i, kk)>=ztglace) THEN
400            zalpha_tr = a_tr_sca(1)
401          ELSE
402            zalpha_tr = a_tr_sca(2)
403          END IF
404          zfrac_lessi = 1. - exp(zalpha_tr * zprec_cond(i) / zneb(i))
405          pfrac_impa(i, kk) = pfrac_impa(i, kk) * (1. - zneb(i) * zfrac_lessi)
406          frac_impa(i, kk) = 1. - zneb(i) * zfrac_lessi
407        END IF
408      END DO
409    END DO
410
411    ! AA----------------------------------------------------------
412    ! FIN DE BOUCLE SUR K
413  END DO
414
415  ! AA-----------------------------------------------------------
416
417  ! Pluie ou neige au sol selon la temperature de la 1ere couche
418
419  DO i = 1, klon
420    IF ((t(i, 1) + d_t(i, 1))<rtt) THEN
421      snow(i) = zrfl(i)
422    ELSE
423      rain(i) = zrfl(i)
424    END IF
425  END DO
426
427END SUBROUTINE fisrtilp_tr
Note: See TracBrowser for help on using the repository browser.