source: dynamico_lmdz/aquaplanet/LMDZ5/libf/phylmd/fisrtilp_tr.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

YM

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