source: LMDZ5/trunk/libf/phylmd/fisrtilp_tr.F90 @ 2260

Last change on this file since 2260 was 1992, checked in by lguez, 11 years ago

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

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