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

Last change on this file since 5151 was 5144, checked in by abarral, 7 weeks ago

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