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

Line Branch Exec Source
1
!
2
! $Id $
3
!
4
SUBROUTINE cvltr_scav(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
5
     sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,    &
6
     pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM,     &
7
     paprs,it,tr,upd,dnd,inb,icb,                   &
8
     ccntrAA_3d,ccntrENV_3d,coefcoli_3d,            &
9
     dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
10
     qPa,qMel,qTrdi,dtrcvMA,Mint,                   &
11
     zmfd1a,zmfphi2,zmfdam)
12
  !
13
  USE IOIPSL
14
  USE dimphy
15
  USE infotrac_phy, ONLY : nbtr
16
  IMPLICIT NONE
17
  !=====================================================================
18
  ! Objet : convection des traceurs / KE
19
  ! Auteurs: M-A Filiberti and J-Y Grandpeix
20
  ! modifiee par R Pilon : lessivage des traceurs / KE
21
  !=====================================================================
22
23
  include "YOMCST.h"
24
  include "YOECUMF.h"
25
  include "conema3.h"
26
  include "chem.h"
27
28
  ! Entree
29
  REAL,INTENT(IN)                           :: pdtime
30
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: da
31
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
32
  ! RomP
33
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam ! matrices pour simplifier
34
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2    ! l'ecriture des tendances
35
  !
36
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: mpIN
37
  REAL,DIMENSION(klon,klev+1),INTENT(IN)    :: paprs  ! pression aux 1/2 couches (bas en haut)
38
  INTEGER,INTENT(IN)                        :: it     ! numero du traceur
39
  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr     ! q de traceur (bas en haut)
40
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: upd   ! saturated updraft mass flux
41
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: dnd   ! saturated downdraft mass flux
42
  !
43
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA   ! masses precipitantes de l'asc adiab
44
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM   ! masses precipitantes des melanges
45
  !JE  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxrIN   ! vprecip: eau
46
  REAL,DIMENSION(klon,klev+1),INTENT(IN)      :: pmflxrIN   ! vprecip: eau
47
  !JE  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxsIN   ! vprecip: neige
48
  REAL,DIMENSION(klon,klev+1),INTENT(IN)      :: pmflxsIN   ! vprecip: neige
49
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: ev         ! evaporation cv30_routine
50
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: epIN
51
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: te
52
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij        ! fraction dair de lenv
53
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wght_cvfd  ! weights of the layers feeding convection
54
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij       ! contenu en eau condensée spécifique/conc deau condensée massique
55
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm    ! eau condensee precipitee dans mel masse dair sat
56
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm    ! eau condensee precipitee dans aa masse dair sat
57
58
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw        ! contenu en eau condensée dans lasc adiab
59
  REAL,DIMENSION(klon),INTENT(IN)           :: sigd
60
  INTEGER,DIMENSION(klon),INTENT(IN)        :: icb,inb
61
  !
62
  REAL,DIMENSION(klon,klev),INTENT(IN) :: ccntrAA_3d
63
  REAL,DIMENSION(klon,klev),INTENT(IN) :: ccntrENV_3d
64
  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefcoli_3d
65
  !
66
  ! Sortie
67
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcv     ! tendance totale (bas en haut)
68
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcvMA ! M-A Filiberti
69
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: trsptd    ! tendance du transport
70
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrSscav  ! tendance du lessivage courant sat
71
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrsat    ! tendance trsp+sat scav
72
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrUscav  ! tendance du lessivage courant unsat
73
  !
74
  ! Variables locales
75
  INTEGER                           :: i,j,k
76
  REAL,DIMENSION(klon,klev)         :: dxpres   ! difference de pression entre niveau (j+1) et (j)
77
  REAL                              :: pdtimeRG ! pas de temps * gravite
78
  REAL,DIMENSION(klon,nbtr)         :: qfeed     ! tracer concentration feeding convection
79
  ! variables pour les courants satures
80
  REAL,DIMENSION(klon,klev,klev)    :: zmd
81
  REAL,DIMENSION(klon,klev,klev)    :: za
82
  REAL,DIMENSION(klon,klev,nbtr)    :: zmfd,zmfa
83
  REAL,DIMENSION(klon,klev,nbtr)    :: zmfp,zmfu
84
85
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfd1a
86
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfdam
87
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfphi2
88
89
  ! RomP ! les variables sont nettoyees des valeurs aberrantes
90
  REAL,DIMENSION(klon,klev)         :: Pa, Pm  ! pluie AA et mélanges, var temporaire
91
  REAL,DIMENSION(klon,klev)         :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
92
  REAL,DIMENSION(klon,klev)         :: mp            ! flux de masse
93
  REAL,DIMENSION(klon,klev)         :: ep            ! fraction d'eau convertie en precipitation
94
  REAL,DIMENSION(klon,klev)         :: evap          ! evaporation : variable temporaire
95
  REAL,DIMENSION(klon,klev)         :: rho    !environmental density
96
97
  REAL,DIMENSION(klon,klev)         :: kappa ! denominateur du au calcul de la matrice
98
  ! pour obtenir qd et qp
99
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qTrdi ! traceurs descente air insature transport MA
100
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qDi  ! traceurs descente insaturees
101
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPr  ! traceurs colonne precipitante
102
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPa  ! traceurs dans les precip issues lasc. adiab.
103
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qMel ! traceurs dans les precip issues des melanges
104
  REAL,DIMENSION(klon,klev,nbtr)                :: qMeltmp ! variable temporaire
105
  REAL,DIMENSION(klon,klev,nbtr)                :: qpmMint
106
  REAL,DIMENSION(klon,klev),INTENT(OUT)         :: Mint
107
  ! tendances
108
  REAL                              :: tdcvMA           ! terme de transport de traceur (schema Marie Angele)
109
  REAL                              :: trsptrac         ! terme de transport de traceur par l'air
110
  REAL                              :: scavtrac         ! terme de lessivage courant sature
111
  REAL                              :: uscavtrac        ! terme de lessivage courant insature
112
  ! impaction
113
!!!       Correction apres discussion Romain P. / Olivier B.
114
!!!  REAL,PARAMETER                    :: rdrop=2.5e-3     ! rayon des gouttes d'eau
115
  REAL,PARAMETER                    :: rdrop=1.e-3     ! rayon des gouttes d'eau
116
!!!
117
  REAL,DIMENSION(klon,klev)         :: imp              ! coefficient d'impaction
118
  !
119
  LOGICAL,DIMENSION(klon,klev) :: NO_precip
120
  ! var tmp tests
121
  REAL                           :: conserv
122
  real                           :: conservMA
123
124
!jyg<
125
!!  ! ======================================================
126
!!  ! calcul de l'impaction
127
!!  ! ======================================================
128
!!
129
!!  ! impaction sur la surface de la colonne de la descente insaturee
130
!!  ! On prend la moyenne des precip entre le niveau i+1 et i
131
!!  ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
132
!!  ! 1000kg/m3= densite de l'eau
133
!!  ! 0.75e-3 = 3/4 /1000
134
!!  ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
135
!!!!  ! on le neglige ici pour simplifier le code
136
!!
137
!!  DO j=1,klev-1
138
!!     DO i=1,klon
139
!!        imp(i,j) = coefcoli_3d(i,j)*0.75e-3/rdrop *&
140
!!             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
141
!!     ENDDO
142
!!  ENDDO
143
!>jyg
144
  !
145
  ! initialisation pour flux de traceurs, td et autre
146
  !
147
  trsptrac = 0.
148
  scavtrac = 0.
149
  uscavtrac = 0.
150
  qfeed(:,it) = 0.              !RL
151
  DO j=1,klev
152
     DO i=1,klon
153
        zmfd(i,j,it)=0.
154
        zmfa(i,j,it)=0.
155
        zmfu(i,j,it)=0.
156
        zmfp(i,j,it)=0.
157
        zmfphi2(i,j,it)=0.
158
        zmfd1a(i,j,it)=0.
159
        zmfdam(i,j,it)=0.
160
        qDi(i,j,it)=0.
161
        qPr(i,j,it)=0.
162
        qPa(i,j,it)=0.
163
        qMel(i,j,it)=0.
164
        qMeltmp(i,j,it)=0.
165
        qTrdi(i,j,it)=0.
166
        kappa(i,j)=0.
167
        trsptd(i,j,it)=0.
168
        dtrsat(i,j,it)=0.
169
        dtrSscav(i,j,it)=0.
170
        dtrUscav(i,j,it)=0.
171
        dtrcv(i,j,it)=0.
172
        dtrcvMA(i,j,it)=0.
173
        evap(i,j)=0.
174
        dxpres(i,j)=0.
175
        qpmMint(i,j,it)=0.
176
        Mint(i,j)=0.
177
     END DO
178
  END DO
179
180
  ! suppression des valeurs très faibles (~1e-320)
181
  ! multiplication de levaporation pour lavoir par unite de temps
182
  ! et par unite de surface de la maille
183
  ! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
184
  DO j=1,klev
185
     DO i=1,klon
186
        IF(ev(i,j).lt.1.e-16) THEN
187
           evap(i,j)=0.
188
        ELSE
189
           evap(i,j)=ev(i,j)*sigd(i)
190
        ENDIF
191
     END DO
192
  END DO
193
194
  DO j=1,klev
195
     DO i=1,klon
196
        IF(j.LT.klev) THEN
197
           IF(epIN(i,j).LT.1.e-32) THEN
198
              ep(i,j)=0.
199
           ELSE
200
              ep(i,j)=epIN(i,j)
201
           ENDIF
202
        ELSE
203
           ep(i,j)=epmax
204
        ENDIF
205
        IF(mpIN(i,j).LT.1.e-32) THEN
206
           mp(i,j)=0.
207
        ELSE
208
           mp(i,j)=mpIN(i,j)
209
        ENDIF
210
        IF(pmflxsIN(i,j).LT.1.e-32) THEN
211
           pmflxs(i,j)=0.
212
        ELSE
213
           pmflxs(i,j)=pmflxsIN(i,j)
214
        ENDIF
215
        IF(pmflxrIN(i,j).LT.1.e-32) THEN
216
           pmflxr(i,j)=0.
217
        ELSE
218
           pmflxr(i,j)=pmflxrIN(i,j)
219
        ENDIF
220
        IF(wdtrainA(i,j).LT.1.e-32) THEN
221
           Pa(i,j)=0.
222
        ELSE
223
           Pa(i,j)=wdtrainA(i,j)
224
        ENDIF
225
        IF(wdtrainM(i,j).LT.1.e-32) THEN
226
           Pm(i,j)=0.
227
        ELSE
228
           Pm(i,j)=wdtrainM(i,j)
229
        ENDIF
230
     END DO
231
  END DO
232
233
  !==========================================
234
  DO j = klev-1,1,-1
235
     DO i = 1,klon
236
        NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).LT.1.e-10&
237
             .AND.Pa(i,j).LT.1.e-10.AND.Pm(i,j).LT.1.e-10
238
     END DO
239
  END DO
240
241
!jyg<
242
  ! ======================================================
243
  ! calcul de l'impaction
244
  ! ======================================================
245
246
  ! impaction sur la surface de la colonne de la descente insaturee
247
  ! On prend la moyenne des precip entre le niveau i+1 et i
248
  ! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
249
  ! 1000kg/m3= densite de l'eau
250
  ! 0.75e-3 = 3/4 /1000
251
  ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
252
  ! on le neglige ici pour simplifier le code
253
254
  DO j=1,klev-1
255
     DO i=1,klon
256
        imp(i,j) = coefcoli_3d(i,j)*0.75e-3/rdrop *&
257
             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
258
     ENDDO
259
  ENDDO
260
!>jyg
261
  ! =========================================
262
  ! calcul des tendances liees au downdraft
263
  ! =========================================
264
  !cdir collapse
265
  DO k=1,klev
266
     DO j=1,klev
267
        DO i=1,klon
268
           zmd(i,j,k)=0.
269
           za (i,j,k)=0.
270
        END DO
271
     END DO
272
  END DO
273
  ! calcul de la matrice d echange
274
  ! matrice de distribution de la masse entrainee en k
275
  ! commmentaire RomP : mp > 0
276
  DO k=1,klev-1
277
     DO i=1,klon
278
        zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1))   ! ~ mk(k)
279
     END DO
280
  END DO
281
  DO k=2,klev
282
     DO j=k-1,1,-1
283
        DO i=1,klon
284
           IF(mp(i,j+1).GT.1.e-10) THEN
285
              zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1)) !det ~ mk(j)=mk(j+1)*mp(i,j)/mp(i,j+1)
286
           ENDIF
287
        END DO
288
     END DO
289
  END DO
290
  DO k=1,klev
291
     DO j=1,klev-1
292
        DO i=1,klon
293
           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
294
        END DO
295
     END DO
296
  END DO
297
!!!!! quantite  de traceur dans la descente d'air insaturee  :   4 juin 2012
298
  DO k=1,klev
299
     DO j=1,klev-1
300
        DO i=1,klon
301
           IF(mp(i,j+1).GT.1.e-10) THEN
302
              qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
303
           ELSE
304
              qTrdi(i,j,it)=0.!tr(i,j,it)
305
           ENDIF
306
        ENDDO
307
     ENDDO
308
  ENDDO
309
!!!!!
310
  !
311
  ! rajout du terme lie a l ascendance induite
312
  !
313
  DO j=2,klev
314
     DO i=1,klon
315
        za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
316
     END DO
317
  END DO
318
  !
319
  ! tendance courants insatures  ! sans lessivage ancien schema
320
  !
321
  DO k=1,klev
322
     DO j=1,klev
323
        DO i=1,klon
324
           zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
325
        END DO
326
     END DO
327
  END DO
328
  !
329
  ! =========================================
330
  ! calcul des tendances liees aux courants satures   j <-> z ; k <-> z'
331
  ! =========================================
332
  !RL
333
  !  Feeding concentrations
334
  DO j=1,klev
335
     DO i=1,klon
336
        qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
337
     END DO
338
  END DO
339
  !RL
340
  !
341
  DO j=1,klev
342
     DO i=1,klon
343
        !RL
344
        !!        zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it))                     ! da
345
        zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it))                     ! da
346
        !RL
347
     END DO
348
  END DO
349
  !
350
  DO k=1,klev
351
     DO j=1,klev
352
        DO i=1,klon
353
           zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it))  ! phi
354
        END DO
355
     END DO
356
  END DO
357
  ! RomP ajout des matrices liees au lessivage
358
  DO j=1,klev
359
     DO i=1,klon
360
        zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it)   ! da1
361
        zmfdam(i,j,it)=dam(i,j)*tr(i,1,it)   ! dam
362
     END DO
363
  END DO
364
  DO k=1,klev
365
     DO j=1,klev
366
        DO i=1,klon
367
           zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it)  ! psi
368
        END DO
369
     END DO
370
  END DO
371
  DO j=1,klev-1
372
     DO i=1,klon
373
        zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
374
     END DO
375
  END DO
376
  DO j=2,klev
377
     DO i=1,klon
378
        zmfu(i,j,it)=zmfu(i,j,it)+min(0.,upd(i,j)+dnd(i,j))*(tr(i,j,it)-tr(i,j-1,it))
379
     END DO
380
  END DO
381
  ! ===================================================
382
  ! calcul des tendances liees aux courants insatures
383
  ! ===================================================
384
  !  pression
385
  DO k=1, klev
386
     DO i=1, klon
387
        dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
388
     ENDDO
389
  ENDDO
390
  pdtimeRG=pdtime*RG
391
392
  ! q_pa et q_pm    traceurs issues des courants satures se retrouvant dans les precipitations
393
  DO j=1,klev
394
     DO i=1,klon
395
        IF(j.GE.icb(i).AND.j.LE.inb(i)) THEN
396
           IF(clw(i,j).GT.1.e-16) THEN
397
              !JE           qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
398
              qPa(i,j,it)=ccntrAA_3d(i,j)*tr(i,1,it)/clw(i,j)
399
           ELSE
400
              qPa(i,j,it)=0.
401
           ENDIF
402
        ENDIF
403
     END DO
404
  END DO
405
406
  ! calcul de q_pm en 2 parties :
407
  ! 1) calcul de sa valeur pour un niveau z' donne
408
  ! 2) integration sur la verticale sur z'
409
  DO j=1,klev
410
     DO k=1,j-1
411
        DO i=1,klon
412
           IF(k.GE.icb(i).AND.k.LE.inb(i).AND.&
413
                j.LE.inb(i)) THEN
414
              IF(elij(i,k,j).GT.1.e-16) THEN
415
                 !JE             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
416
                 !JE                         *(1.-sij(i,k,j))  +ccntrENV_coef&
417
                 !JE                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
418
                 qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_3d(i,k)*tr(i,1,it)&
419
                      *(1.-sij(i,k,j))  +ccntrENV_3d(i,k)&
420
                      *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
421
              ELSE
422
                 qMeltmp(i,j,it)=0.
423
              ENDIF
424
              qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
425
              Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
426
           ENDIF ! end if dans nuage
427
        END DO
428
     END DO
429
  END DO
430
431
  DO j=1,klev
432
     DO i=1,klon
433
        IF(Mint(i,j).GT.1.e-16) THEN
434
           qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
435
        ELSE
436
           qMel(i,j,it)=0.
437
        ENDIF
438
     END DO
439
  END DO
440
441
  ! calcul de q_d et q_p    traceurs de la descente precipitante
442
  DO j=klev-1,1,-1
443
     DO i=1,klon
444
        IF(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN  ! detrainement
445
           kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
446
                (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
447
                + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
448
449
        ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN! entrainement
450
           IF(j.eq.1) THEN
451
              kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
452
                   (-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
453
                   + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
454
           ELSE
455
              kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
456
                   (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
457
                   + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
458
           ENDIF
459
        ELSE
460
           kappa(i,j)=1.
461
        ENDIF
462
     ENDDO
463
  ENDDO
464
465
  DO j=klev-1,1,-1
466
     DO i=1,klon
467
        IF (abs(kappa(i,j)).LT.1.e-25) THEN    !si denominateur nul (il peut y avoir des mp!=0)
468
           kappa(i,j)=1.
469
           IF(j.eq.1) THEN
470
              qDi(i,j,it)=qDi(i,j+1,it) !orig tr(i,j,it)   ! mp(1)=0 donc tout vient de la couche supérieure
471
           ELSEIF(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN
472
              qDi(i,j,it)=qDi(i,j+1,it)
473
           ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN! entrainement
474
              qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j))
475
           ELSE  ! si  mp (i)=0 et mp(j+1)=0
476
              qDi(i,j,it)=tr(i,j,it) ! orig 0.
477
           ENDIF
478
479
           IF(NO_precip(i,j)) THEN
480
              qPr(i,j,it)=0.
481
           ELSE
482
              qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
483
                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
484
                   +imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/&
485
                   (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
486
           ENDIF
487
        ELSE   !     denominateur non nul
488
           kappa(i,j)=1./kappa(i,j)
489
           ! calcul de qd et qp
490
           !!jyg  (20130119) correction pour le sommet du nuage
491
           !!     if(j.GE.inb(i)) THEN       !au-dessus du nuage, sommet inclu
492
           if(j.GT.inb(i)) THEN       !au-dessus du nuage
493
              qDi(i,j,it)=tr(i,j,it)   ! pas de descente => environnement = descente insaturee
494
              qPr(i,j,it)=0.
495
496
              !  vvv premiere couche du modele ou mp(1)=0  ! det tout le temps  vvv
497
           ELSEIF(j.eq.1) THEN
498
              if(mp(i,2).GT.1.e-10) THEN !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
499
                 if(NO_precip(i,j)) THEN !pas de precip en (i)
500
                    qDi(i,j,it)=qDi(i,j+1,it)
501
                    qPr(i,j,it)=0.
502
                 ELSE
503
                    qDi(i,j,it)=kappa(i,j)*(&
504
                         (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
505
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
506
                         (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
507
                         (-mp(i,j+1)*qDi(i,j+1,it)))
508
509
                    qPr(i,j,it)=kappa(i,j)*(&
510
                         (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
511
                         ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
512
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
513
                         +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
514
                 ENDIF
515
516
              ELSE !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
517
                 qDi(i,j,it)=tr(i,j,it) ! orig 0.
518
                 if(NO_precip(i,j)) THEN
519
                    qPr(i,j,it)=0.
520
                 ELSE
521
                    qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
522
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
523
                         +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
524
                         (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
525
                 ENDIF
526
527
              ENDIF  !mp(2) nul ou non
528
529
              ! vvv  (j!=1.AND.j.LT.inb(i))  en-dessous du sommet nuage   vvv
530
           ELSE
531
              !------------------------------------------------------------- detrainement
532
              if(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN !mp(i,j).GT.1.e-10) THEN
533
                 if(NO_precip(i,j)) THEN
534
                    qDi(i,j,it)=qDi(i,j+1,it)
535
                    qPr(i,j,it)=0.
536
                 ELSE
537
                    qDi(i,j,it)=kappa(i,j)*(&
538
                         (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
539
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
540
                         (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
541
                         (-mp(i,j+1)*qDi(i,j+1,it)))
542
                    !
543
                    qPr(i,j,it)=kappa(i,j)*(&
544
                         (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
545
                         ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
546
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
547
                         +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
548
                 ENDIF !precip
549
                 !------------------------------------------------------------- entrainement
550
              ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN
551
                 if(NO_precip(i,j)) THEN
552
                    qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j))
553
                    qPr(i,j,it)=0.
554
                 ELSE
555
                    qDi(i,j,it)=kappa(i,j)*(&
556
                         (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
557
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
558
                         (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
559
                         (-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it)))
560
                    !
561
                    qPr(i,j,it)=kappa(i,j)*(&
562
                         (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*&
563
                         ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
564
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
565
                         +(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*&
566
                         (imp(i,j)/RG*dxpres(i,j)))
567
                 ENDIF !precip
568
                 !------------------------------------------------------------- ENDIF ! ent/det
569
              ELSE  !mp nul
570
                 qDi(i,j,it)=tr(i,j,it) ! orig 0.
571
                 if(NO_precip(i,j)) THEN
572
                    qPr(i,j,it)=0.
573
                 ELSE
574
                    qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
575
                         Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
576
                         +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
577
                         (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
578
                 ENDIF
579
              ENDIF ! mp nul ou non
580
           ENDIF ! condition sur j
581
        ENDIF ! kappa
582
     ENDDO
583
  ENDDO
584
585
  !! print test descente insaturee
586
  !  DO j=klev,1,-1
587
  !   DO i=1,klon
588
  !     if(it.eq.3) THEN
589
  !   write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') j,&
590
  !!         'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
591
  !          'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
592
  !          'Pm',Pm(i,j),'Mint',Mint(i,j),&
593
  !!          'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
594
  !        'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
595
  !!          'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
596
  !!         'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
597
  !     ENDIF
598
  !   ENDDO
599
  !  ENDDO
600
601
602
  ! ===================================================
603
  ! calcul final des tendances
604
  ! ===================================================
605
606
  DO k=klev-1,1,-1
607
     DO i=1, klon
608
        ! transport
609
        tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)   ! double comptage des downdraft insatures
610
        trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
611
        ! lessivage courants satures
612
        !JE     scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
613
        !JE               -zmfphi2(i,k,it)*ccntrENV_coef&
614
        !JE               -zmfdam(i,k,it)*ccntrAA_coef
615
        scavtrac=-ccntrAA_3d(i,k)*zmfd1a(i,k,it)&
616
             -zmfphi2(i,k,it)*ccntrENV_3d(i,k)&
617
             -zmfdam(i,k,it)*ccntrAA_3d(i,k)
618
        ! lessivage courants insatures
619
        if(k.LE.inb(i).AND.k.GT.1) THEN   ! tendances dans le nuage
620
           !------------------------------------------------------------- detrainement
621
           if(mp(i,k+1).GT.mp(i,k).AND.mp(i,k+1).GT.1.e-10) THEN
622
              uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
623
                   + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
624
              !
625
              !        if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
626
              !                    (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
627
              !                    mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
628
              !                    'mp',mp(i,k)
629
              !------------------------------------------------------------- entrainement
630
           ELSEIF(mp(i,k).GT.mp(i,k+1).AND.mp(i,k).GT.1.e-10) THEN
631
              uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
632
              !
633
              !         if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
634
              !=!------------------------------------------------------------- end ent/det
635
           ELSE !        mp(i,k+1)=0. et mp(i,k)=0.    pluie directement sur l environnement
636
637
              if(NO_precip(i,k)) THEN
638
                 uscavtrac=0.
639
                 !         if(it.eq.3) write(*,'(I2,1X,a,e20.12,82X,a,e20.12)')k,' no P ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
640
              ELSE
641
                 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
642
                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
643
                 !!JE adds
644
                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
645
                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'tr',tr(i,k,it)
646
                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
647
                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'qPr',qPr(i,k,it)
648
                 !         if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
649
                 !!Je end
650
651
              ENDIF
652
           ENDIF  ! mp/det/ent
653
           !------------------------------------------------------------- premiere couche
654
        ELSEIF(k.eq.1) THEN  !                                      mp(1)=0.
655
           if(mp(i,2).GT.1.e-10) THEN  !detrainement
656
              uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
657
              !                   + mp(i,2)*(0.-tr(i,k,it))
658
              !
659
              !       if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
660
              !                  (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
661
              !                  mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
662
              !                   'mp',mp(i,k)
663
           ELSE   ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
664
              if(NO_precip(i,1)) THEN
665
                 uscavtrac=0.
666
              ELSE
667
                 uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
668
              ENDIF
669
              !         if(it.eq.3) write(*,'(I2,1X,a,2X,e20.12,82X,a,e20.12)')k,'1 P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k)
670
           ENDIF
671
672
        ELSE   ! k > INB  au-dessus du nuage
673
           uscavtrac=0.
674
        ENDIF
675
676
        ! =====    tendances finales  ======
677
        trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k)              ! td transport sans eau dans courants satures
678
        dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k)            ! td du lessivage dans courants satures
679
        dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k)           ! td courant insat
680
        dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k)   ! td courant sat
681
        dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it)         td conv
682
!!!!!!
683
        dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
684
     ENDDO
685
  ENDDO
686
687
  ! test de conservation du traceur
688
  !print*,"_____________________________________________________________"
689
  !print*,"                                                             "
690
  !      conserv=0.
691
  !      conservMA=0.
692
  !      DO k= klev-1,1,-1
693
  !        DO i=1, klon
694
  !         conserv=conserv+dtrcv(i,k,it)*   &
695
  !        (paprs(i,k)-paprs(i,k+1))/RG
696
  !         conservMA=conservMA+dtrcvMA(i,k,it)*   &
697
  !        (paprs(i,k)-paprs(i,k+1))/RG
698
  !
699
  !      if(it.eq.3)  write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
700
  !         'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
701
  !         ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,'         conservMA ',conservMA,'conserv ',conserv
702
  !!
703
  !        ENDDO
704
  !      ENDDO
705
  !       if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
706
707
END SUBROUTINE cvltr_scav