source: LMDZ5/branches/testing/libf/phylmd/fisrtilp_tr.F90 @ 3632

Last change on this file since 3632 was 2408, checked in by Laurent Fairhead, 9 years ago

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