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

Last change on this file since 5157 was 5153, checked in by abarral, 3 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
RevLine 
[1403]1! $Id: fisrtilp_tr.F90 5153 2024-07-31 16:20:03Z abarral $
[524]2
3
[1992]4SUBROUTINE fisrtilp_tr(dtime, paprs, pplay, t, q, ratqs, d_t, d_q, d_ql, &
[5144]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
[1992]7  ! properties; aeropt.F)
8
9  USE dimphy
[5112]10  USE lmdz_print_control, ONLY: lunout
[5144]11  USE lmdz_yoethf
[5153]12
[5144]13  USE lmdz_yomcst
[5143]14
[1992]15  IMPLICIT NONE
[5153]16 INCLUDE "FCTTRE.h"
[1992]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)
[5144]28  REAL paprs(klon, klev + 1) ! pression a inter-couche
[1992]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)
[5144]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)
[1992]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
[5144]62  PARAMETER (seuil_neb = 0.001)
[1992]63  REAL ct ! inverse du temps pour qu'un nuage precipite
[5144]64  PARAMETER (ct = 1. / 1800.)
[1992]65  REAL cl ! seuil de precipitation
[5144]66  PARAMETER (cl = 2.6E-4)
[1992]67  ! cc      PARAMETER (cl=2.3e-4)
68  ! cc      PARAMETER (cl=2.0e-4)
69  INTEGER ninter ! sous-intervals pour la precipitation
[5144]70  PARAMETER (ninter = 5)
[1992]71  LOGICAL evap_prec ! evaporation de la pluie
[5144]72  PARAMETER (evap_prec = .TRUE.)
[1992]73  REAL coef_eva
[5144]74  PARAMETER (coef_eva = 2.0E-05)
[1992]75  LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
76  REAL ratqs(klon, klev) ! determine la largeur de distribution de vapeur
[5144]77  PARAMETER (calcrat = .TRUE.)
[1992]78  REAL zx_min, rat_max
[5144]79  PARAMETER (zx_min = 1.0, rat_max = 0.01)
[1992]80  REAL zx_max, rat_min
[5144]81  PARAMETER (zx_max = 0.1, rat_min = 0.3)
[1992]82  REAL zx
83
84  LOGICAL cpartiel ! condensation partielle
[5144]85  PARAMETER (cpartiel = .TRUE.)
[1992]86  REAL t_coup
[5144]87  PARAMETER (t_coup = 234.0)
[1992]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
[5144]126  DATA appel1er/.TRUE./
127  fallv(zzz) = 3.29 / 2.0 * ((zzz)**0.16)
[1992]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
[5144]137    IF (abs(dtime / real(ninter) - 360.0)>0.001) THEN
[1992]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
[524]153      DO i = 1, klon
[1992]154        pfrac_nucl(i, k) = 1.
155        pfrac_1nucl(i, k) = 1.
156        pfrac_impa(i, k) = 1.
157      END DO
158    END DO
[524]159
[1992]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
[524]226      DO i = 1, klon
[1992]227        IF (zrfl(i)>0.) THEN
228          IF (thermcep) THEN
[5144]229            zdelta = max(0., sign(1., rtt - zt(i)))
230            zqs(i) = r2es * foeew(zt(i), zdelta) / pplay(i, k)
[1992]231            zqs(i) = min(0.5, zqs(i))
[5144]232            zcor = 1. / (1. - retv * zqs(i))
233            zqs(i) = zqs(i) * zcor
[1992]234          ELSE
235            IF (zt(i)<t_coup) THEN
[5144]236              zqs(i) = qsats(zt(i)) / pplay(i, k)
[1992]237            ELSE
[5144]238              zqs(i) = qsatl(zt(i)) / pplay(i, k)
[1992]239            END IF
240          END IF
[5144]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))
[1992]246          zqev = min(zqev, zqevt)
[5144]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))
[1992]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
[524]260      DO i = 1, klon
[5144]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)
[1992]265        zqs(i) = min(0.5, zqs(i))
[5144]266        zcor = 1. / (1. - retv * zqs(i))
267        zqs(i) = zqs(i) * zcor
[1992]268        zdqs(i) = foede(zt(i), zdelta, zcvm5, zqs(i), zcor)
269      END DO
270    ELSE
[524]271      DO i = 1, klon
[1992]272        IF (zt(i)<t_coup) THEN
[5144]273          zqs(i) = qsats(zt(i)) / pplay(i, k)
[1992]274          zdqs(i) = dqsats(zt(i), zqs(i))
275        ELSE
[5144]276          zqs(i) = qsatl(zt(i)) / pplay(i, k)
[1992]277          zdqs(i) = dqsatl(zt(i), zqs(i))
278        END IF
279      END DO
280    END IF
[524]281
[1992]282    ! Determiner la condensation partielle et calculer la quantite
283    ! de l'eau condensee:
[524]284
[1992]285    IF (cpartiel) THEN
[524]286      DO i = 1, klon
[1992]287
[5144]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))
[1992]295
296        ! --Olivier
[5144]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
[1992]300        ! --fin
301
302      END DO
303    ELSE
[524]304      DO i = 1, klon
[1992]305        IF (zq(i)>zqs(i)) THEN
306          rneb(i, k) = 1.0
307        ELSE
308          rneb(i, k) = 0.0
309        END IF
[5144]310        zcond(i) = max(0.0, zq(i) - zqs(i)) / (1. + zdqs(i))
[1992]311      END DO
312    END IF
313
314    DO i = 1, klon
315      zq(i) = zq(i) - zcond(i)
[5144]316      zt(i) = zt(i) + zcond(i) * rlvtt / rcpd
[1992]317    END DO
318
319    ! Partager l'eau condensee en precipitation et eau liquide nuageuse
320
321    DO i = 1, klon
[5144]322      IF (rneb(i, k)>0.0) THEN
[1992]323        zoliq(i) = zcond(i)
[5144]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)
[1992]328        zfice(i) = zfice(i)**nexpo
[5144]329        zneb(i) = max(rneb(i, k), seuil_neb)
330        radliq(i, k) = zoliq(i) / real(ninter + 1)
[1992]331      END IF
332    END DO
333
334    DO n = 1, ninter
[524]335      DO i = 1, klon
[5144]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)
[1992]342          ztot(i) = zchau(i) + zfroi(i)
343          IF (zneb(i)==seuil_neb) ztot(i) = 0.0
[5144]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)
[1992]347        END IF
348      END DO
349    END DO
350
351    DO i = 1, klon
[5144]352      IF (rneb(i, k)>0.0) THEN
[1992]353        d_ql(i, k) = zoliq(i)
[5144]354        zrfl(i) = zrfl(i) + max(zcond(i) - zoliq(i), 0.0) * (paprs(i, k) - paprs(i, k &
355                + 1)) / (rg * dtime)
[1992]356      END IF
357      IF (zt(i)<rtt) THEN
358        psfl(i, k) = zrfl(i)
[524]359      ELSE
[1992]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
[5144]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
[1992]378        ! AA lessivage nucleation LMD5 dans la couche elle-meme
[5144]379        IF (t(i, k)>=ztglace) THEN
[1992]380          zalpha_tr = a_tr_sca(3)
381        ELSE
382          zalpha_tr = a_tr_sca(4)
383        END IF
[5144]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
[1992]387
388        ! nucleation avec un facteur -1 au lieu de -0.5
[5144]389        zfrac_lessi = 1. - exp(-zprec_cond(i) / zneb(i))
390        pfrac_1nucl(i, k) = pfrac_1nucl(i, k) * (1. - zneb(i) * zfrac_lessi)
[1992]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
[524]397      DO i = 1, klon
[5144]398        IF (rneb(i, k)>0.0 .AND. zprec_cond(i)>0.) THEN
399          IF (t(i, kk)>=ztglace) THEN
[1992]400            zalpha_tr = a_tr_sca(1)
401          ELSE
402            zalpha_tr = a_tr_sca(2)
403          END IF
[5144]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
[1992]407        END IF
408      END DO
409    END DO
[524]410
[1992]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
[5144]420    IF ((t(i, 1) + d_t(i, 1))<rtt) THEN
[1992]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.