GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: phylmd/fisrtilp_tr.F90 Lines: 0 148 0.0 %
Date: 2023-06-30 12:56:34 Branches: 0 86 0.0 %

Line Branch Exec Source
1
2
! $Id: fisrtilp_tr.F90 2346 2015-08-21 15:13:46Z emillour $
3
4
5
SUBROUTINE 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
  USE print_control_mod, ONLY: lunout
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
154
      DO i = 1, klon
155
        pfrac_nucl(i, k) = 1.
156
        pfrac_1nucl(i, k) = 1.
157
        pfrac_impa(i, k) = 1.
158
      END DO
159
    END DO
160
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
227
      DO i = 1, klon
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
261
      DO i = 1, klon
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
272
      DO i = 1, klon
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
282
283
    ! Determiner la condensation partielle et calculer la quantite
284
    ! de l'eau condensee:
285
286
    IF (cpartiel) THEN
287
      DO i = 1, klon
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
305
      DO i = 1, klon
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
336
      DO i = 1, klon
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)
360
      ELSE
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
398
      DO i = 1, klon
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
411
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
429
END SUBROUTINE fisrtilp_tr