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

Line Branch Exec Source
1
!
2
! $Id $
3
!
4
SUBROUTINE cvltr(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
5
!!           sigd,sij,clw,elij,epmlmMm,eplaMm,              &   !RL
6
           sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm,    &     !RL
7
           pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM,     &
8
           paprs,it,tr,upd,dnd,inb,icb,                   &
9
           dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
10
           qPa,qMel,qTrdi,dtrcvMA,Mint,                   &
11
           zmfd1a,zmfphi2,zmfdam)
12
  USE IOIPSL
13
  USE dimphy
14
  USE infotrac_phy, ONLY : nbtr
15
  IMPLICIT NONE
16
!=====================================================================
17
! Objet : convection des traceurs / KE
18
! Auteurs: M-A Filiberti and J-Y Grandpeix
19
! modifiee par R Pilon : lessivage des traceurs / KE
20
!=====================================================================
21
22
  include "YOMCST.h"
23
  include "YOECUMF.h"
24
  include "conema3.h"
25
26
! Entree
27
  REAL,INTENT(IN)                           :: pdtime
28
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: da
29
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
30
! RomP
31
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam ! matrices pour simplifier
32
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2    ! l'ecriture des tendances
33
!
34
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: mpIN
35
  REAL,DIMENSION(klon,klev+1),INTENT(IN)    :: paprs  ! pression aux 1/2 couches (bas en haut)
36
!  REAL,DIMENSION(klon,klev),INTENT(IN)    :: pplay ! pression aux 1/2 couches (bas en haut)
37
  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr     ! q de traceur (bas en haut)
38
  INTEGER,INTENT(IN)                        :: it
39
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: upd   ! saturated updraft mass flux
40
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: dnd   ! saturated downdraft mass flux
41
!
42
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA   ! masses precipitantes de l'asc adiab
43
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM   ! masses precipitantes des melanges
44
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxrIN   ! vprecip: eau
45
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxsIN   ! vprecip: neige
46
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: ev         ! evaporation cv30_routine
47
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: epIN
48
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: te
49
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij        ! fraction dair de lenv
50
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wght_cvfd  ! weights of the layers feeding convection
51
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij       ! contenu en eau condens�e sp�cifique/conc deau condens�e massique
52
  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm    ! eau condensee precipitee dans mel masse dair sat
53
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm    ! eau condensee precipitee dans aa masse dair sat
54
55
  REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw        ! contenu en eau condens�e dans lasc adiab
56
  REAL,DIMENSION(klon),INTENT(IN)           :: sigd
57
  INTEGER,DIMENSION(klon),INTENT(IN)        :: icb,inb
58
! Sortie
59
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcv     ! tendance totale (bas en haut)
60
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcvMA ! M-A Filiberti
61
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: trsptd    ! tendance du transport
62
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrSscav  ! tendance du lessivage courant sat
63
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrsat    ! tendance trsp+sat scav
64
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrUscav  ! tendance du lessivage courant unsat
65
!
66
! Variables locales
67
  INTEGER                           :: i,j,k
68
  REAL,DIMENSION(klon,klev)         :: dxpres   ! difference de pression entre niveau (j+1) et (j)
69
  REAL                              :: pdtimeRG ! pas de temps * gravite
70
! variables pour les courants satures
71
  REAL,DIMENSION(klon,klev,klev)    :: zmd
72
  REAL,DIMENSION(klon,klev,klev)    :: za
73
  REAL,DIMENSION(klon,klev,nbtr)    :: zmfd,zmfa
74
  REAL,DIMENSION(klon,klev,nbtr)    :: zmfp,zmfu
75
  REAL,DIMENSION(klon,nbtr)         :: qfeed     ! tracer concentration feeding convection
76
77
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfd1a
78
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfdam
79
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfphi2
80
81
! RomP ! les variables sont nettoyees des valeurs aberrantes
82
  REAL,DIMENSION(klon,klev)         :: Pa, Pm  ! pluie AA et m�langes, var temporaire
83
  REAL,DIMENSION(klon,klev)         :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
84
  REAL,DIMENSION(klon,klev)         :: mp            ! flux de masse
85
  REAL,DIMENSION(klon,klev)         :: ep            ! fraction d'eau convertie en precipitation
86
  REAL,DIMENSION(klon,klev)         :: evap          ! evaporation : variable temporaire
87
  REAL,DIMENSION(klon,klev)         :: rho    !environmental density
88
89
  REAL,DIMENSION(klon,klev)         :: kappa ! denominateur du au calcul de la matrice
90
                                             ! pour obtenir qd et qp
91
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qTrdi ! traceurs descente air insature transport MA
92
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qDi  ! traceurs descente insaturees
93
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPr  ! traceurs colonne precipitante
94
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPa  ! traceurs dans les precip issues lasc. adiab.
95
  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qMel ! traceurs dans les precip issues des melanges
96
  REAL,DIMENSION(klon,klev,nbtr)                :: qMeltmp ! variable temporaire
97
  REAL,DIMENSION(klon,klev,nbtr)                :: qpmMint
98
  REAL,DIMENSION(klon,klev),INTENT(OUT)         :: Mint
99
! tendances
100
  REAL                              :: tdcvMA           ! terme de transport de traceur (schema Marie Angele)
101
  REAL                              :: trsptrac         ! terme de transport de traceur par l'air
102
  REAL                              :: scavtrac         ! terme de lessivage courant sature
103
  REAL                              :: uscavtrac        ! terme de lessivage courant insature
104
! impaction
105
!!!       Correction apres discussion Romain P. / Olivier B.
106
!!!  REAL,PARAMETER                    :: rdrop=2.5e-3     ! rayon des gouttes d'eau
107
  REAL,PARAMETER                    :: rdrop=1.e-3     ! rayon des gouttes d'eau
108
!!!
109
  REAL,DIMENSION(klon,klev)         :: imp              ! coefficient d'impaction
110
! parametres lessivage
111
  REAL                    :: ccntrAA_coef        ! \alpha_a : fract aerosols de l'AA convertis en CCN
112
  REAL                    :: ccntrENV_coef       ! \beta_m  : fract aerosols de l'env convertis en CCN
113
  REAL                    :: coefcoli            ! coefficient de collision des gouttes sur les aerosols
114
!
115
  LOGICAL,DIMENSION(klon,klev) :: NO_precip
116
!  LOGICAL                      :: scavON
117
! var tmp tests
118
  REAL                           :: conserv
119
  real                           :: conservMA
120
121
! coefficient lessivage
122
   ccntrAA_coef     = 0.
123
   ccntrENV_coef    = 0.
124
   coefcoli         = 0.
125
126
!$OMP MASTER
127
  call getin('ccntrAA_coef',ccntrAA_coef)
128
  call getin('ccntrENV_coef',ccntrENV_coef)
129
  call getin('coefcoli',coefcoli)
130
!$OMP END MASTER
131
!$OMP BARRIER
132
  print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
133
134
!  scavON=.TRUE.
135
!  if(scavON) then
136
!   ccntrAA_coef     = 1.
137
!   ccntrENV_coef    = 1.
138
!   coefcoli         = 1.
139
!  else
140
!   ccntrAA_coef     = 0.
141
!   ccntrENV_coef    = 0.
142
!   coefcoli         = 0.
143
!  endif
144
145
! ======================================================
146
! calcul de l'impaction
147
! ======================================================
148
!initialisation
149
  do j=1,klev
150
   do i=1,klon
151
     imp(i,j)=0.
152
   enddo
153
  enddo
154
! impaction sur la surface de la colonne de la descente insaturee
155
! On prend la moyenne des precip entre le niveau i+1 et i
156
! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
157
!  1000kg/m3= densit� de l'eau
158
! 0.75e-3 = 3/4 /1000
159
! Par la suite, I est tout le temps multipli� par sig_d pour avoir l'impaction sur la surface de la maille
160
! on le n�glige ici pour simplifier le code
161
  do j=1,klev-1
162
   do i=1,klon
163
    imp(i,j) = coefcoli*0.75e-3/rdrop *&
164
             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
165
!    rho(i,j)=pplay(i,j)/(rd*te(i,j))
166
   enddo
167
  enddo
168
!
169
! initialisation pour flux de traceurs, td et autre
170
   trsptrac = 0.
171
   scavtrac = 0.
172
   uscavtrac = 0.
173
   qfeed(:,it) = 0.              !RL
174
  DO j=1,klev
175
   DO i=1,klon
176
    zmfd(i,j,it)=0.
177
    zmfa(i,j,it)=0.
178
    zmfu(i,j,it)=0.
179
    zmfp(i,j,it)=0.
180
    zmfphi2(i,j,it)=0.
181
    zmfd1a(i,j,it)=0.
182
    zmfdam(i,j,it)=0.
183
    qDi(i,j,it)=0.
184
    qPr(i,j,it)=0.
185
    qPa(i,j,it)=0.
186
    qMel(i,j,it)=0.
187
    qMeltmp(i,j,it)=0.
188
    qTrdi(i,j,it)=0.
189
    kappa(i,j)=0.
190
    trsptd(i,j,it)=0.
191
    dtrsat(i,j,it)=0.
192
    dtrSscav(i,j,it)=0.
193
    dtrUscav(i,j,it)=0.
194
    dtrcv(i,j,it)=0.
195
    dtrcvMA(i,j,it)=0.
196
    evap(i,j)=0.
197
    dxpres(i,j)=0.
198
    qpmMint(i,j,it)=0.
199
    Mint(i,j)=0.
200
   END DO
201
  END DO
202
203
! suppression des valeurs tr�s faibles (~1e-320)
204
! multiplication de levaporation pour lavoir par unite de temps
205
! et par unite de surface de la maille
206
! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
207
  DO j=1,klev
208
   DO i=1,klon
209
    if(ev(i,j).lt.1.e-16) then
210
     evap(i,j)=0.
211
    else
212
     evap(i,j)=ev(i,j)*sigd(i)
213
    endif
214
   END DO
215
  END DO
216
217
  DO j=1,klev
218
   DO i=1,klon
219
   if(j.lt.klev) then
220
    if(epIN(i,j).lt.1.e-32) then
221
     ep(i,j)=0.
222
    else
223
     ep(i,j)=epIN(i,j)
224
    endif
225
   else
226
    ep(i,j)=epmax
227
   endif
228
    if(mpIN(i,j).lt.1.e-32) then
229
     mp(i,j)=0.
230
    else
231
     mp(i,j)=mpIN(i,j)
232
    endif
233
    if(pmflxsIN(i,j).lt.1.e-32) then
234
     pmflxs(i,j)=0.
235
    else
236
     pmflxs(i,j)=pmflxsIN(i,j)
237
    endif
238
    if(pmflxrIN(i,j).lt.1.e-32) then
239
     pmflxr(i,j)=0.
240
    else
241
     pmflxr(i,j)=pmflxrIN(i,j)
242
    endif
243
    if(wdtrainA(i,j).lt.1.e-32) then
244
     Pa(i,j)=0.
245
    else
246
     Pa(i,j)=wdtrainA(i,j)
247
    endif
248
    if(wdtrainM(i,j).lt.1.e-32) then
249
     Pm(i,j)=0.
250
    else
251
     Pm(i,j)=wdtrainM(i,j)
252
    endif
253
   END DO
254
  END DO
255
256
!==========================================
257
  DO j = klev-1,1,-1
258
   DO i = 1,klon
259
     NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
260
                    .and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
261
   END DO
262
  END DO
263
264
! =========================================
265
! calcul des tendances liees au downdraft
266
! =========================================
267
!cdir collapse
268
  DO k=1,klev
269
   DO j=1,klev
270
    DO i=1,klon
271
     zmd(i,j,k)=0.
272
     za (i,j,k)=0.
273
    END DO
274
   END DO
275
  END DO
276
! calcul de la matrice d echange
277
! matrice de distribution de la masse entrainee en k
278
! commmentaire RomP : mp > 0
279
  DO k=1,klev-1
280
     DO i=1,klon
281
        zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1))   ! ~ mk(k)
282
     END DO
283
  END DO
284
  DO k=2,klev
285
     DO j=k-1,1,-1
286
        DO i=1,klon
287
           if(mp(i,j+1).gt.1.e-10) then
288
              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)
289
           ENDif
290
        END DO
291
     END DO
292
  END DO
293
  DO k=1,klev
294
     DO j=1,klev-1
295
        DO i=1,klon
296
           za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
297
        END DO
298
     END DO
299
  END DO
300
!!!!! quantite  de traceur dans la descente d'air insaturee  :   4 juin 2012
301
  DO k=1,klev
302
     DO j=1,klev-1
303
        DO i=1,klon
304
        if(mp(i,j+1).gt.1.e-10) then
305
          qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
306
        else
307
          qTrdi(i,j,it)=0.!tr(i,j,it)
308
        endif
309
        ENDDO
310
     ENDDO
311
  ENDDO
312
!!!!!
313
!
314
! rajout du terme lie a l ascendance induite
315
!
316
  DO j=2,klev
317
     DO i=1,klon
318
        za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
319
     END DO
320
  END DO
321
!
322
! tendance courants insatures  ! sans lessivage ancien schema
323
!
324
  DO k=1,klev
325
     DO j=1,klev
326
        DO i=1,klon
327
           zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
328
        END DO
329
     END DO
330
  END DO
331
!
332
! =========================================
333
! calcul des tendances liees aux courants satures   j <-> z ; k <-> z'
334
! =========================================
335
!
336
!RL
337
!  Feeding concentrations
338
  DO j=1,klev
339
     DO i=1,klon
340
        qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
341
     END DO
342
  END DO
343
!RL
344
!
345
  DO j=1,klev
346
     DO i=1,klon
347
!RL
348
!!        zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it))                     ! da
349
        zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it))                     ! da
350
!RL
351
     END DO
352
  END DO
353
!
354
  DO k=1,klev
355
     DO j=1,klev
356
        DO i=1,klon
357
           zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it))  ! phi
358
        END DO
359
     END DO
360
  END DO
361
! RomP ajout des matrices liees au lessivage
362
  DO j=1,klev
363
     DO i=1,klon
364
        zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it)   ! da1
365
        zmfdam(i,j,it)=dam(i,j)*tr(i,1,it)   ! dam
366
     END DO
367
  END DO
368
  DO k=1,klev
369
     DO j=1,klev
370
        DO i=1,klon
371
          zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it)  ! psi
372
        END DO
373
     END DO
374
  END DO
375
  DO j=1,klev-1
376
     DO i=1,klon
377
        zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
378
     END DO
379
  END DO
380
  DO j=2,klev
381
     DO i=1,klon
382
        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))
383
     END DO
384
  END DO
385
! ===================================================
386
! calcul des tendances liees aux courants insatures
387
! ===================================================
388
!  pression
389
  DO k=1, klev
390
     DO i=1, klon
391
        dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
392
     ENDDO
393
  ENDDO
394
  pdtimeRG=pdtime*RG
395
396
! q_pa et q_pm    traceurs issues des courants satures se retrouvant dans les precipitations
397
  DO j=1,klev
398
     DO i=1,klon
399
        if(j.ge.icb(i).and.j.le.inb(i)) then
400
          if(clw(i,j).gt.1.e-16) then
401
           qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
402
          else
403
           qPa(i,j,it)=0.
404
          endif
405
        endif
406
     END DO
407
  END DO
408
409
! calcul de q_pm en 2 parties :
410
! 1) calcul de sa valeur pour un niveau z' donne
411
! 2) integration sur la verticale sur z'
412
     DO j=1,klev
413
      DO k=1,j-1
414
        DO i=1,klon
415
        if(k.ge.icb(i).and.k.le.inb(i).and.&
416
           j.le.inb(i)) then
417
            if(elij(i,k,j).gt.1.e-16) then
418
             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
419
                         *(1.-sij(i,k,j))  +ccntrENV_coef&
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
     scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
613
               -zmfphi2(i,k,it)*ccntrENV_coef&
614
               -zmfdam(i,k,it)*ccntrAA_coef
615
! lessivage courants insatures
616
   if(k.le.inb(i).and.k.gt.1) then   ! tendances dans le nuage
617
!------------------------------------------------------------- detrainement
618
      if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
619
        uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
620
                   + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
621
!
622
!        if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
623
!                    (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
624
!                    mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
625
!                    'mp',mp(i,k)
626
!------------------------------------------------------------- entrainement
627
      elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
628
         uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
629
!
630
!         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)
631
!=!------------------------------------------------------------- end ent/det
632
      else !        mp(i,k+1)=0. et mp(i,k)=0.    pluie directement sur l environnement
633
634
        if(NO_precip(i,k)) then
635
          uscavtrac=0.
636
!         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)
637
        else
638
          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
639
!         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)
640
        endif
641
      endif  ! mp/det/ent
642
!------------------------------------------------------------- premiere couche
643
   elseif(k.eq.1) then  !                                      mp(1)=0.
644
      if(mp(i,2).gt.1.e-10) then  !detrainement
645
         uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
646
!                   + mp(i,2)*(0.-tr(i,k,it))
647
!
648
!       if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
649
!                  (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
650
!                  mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
651
!                   'mp',mp(i,k)
652
      else   ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
653
        if(NO_precip(i,1)) then
654
          uscavtrac=0.
655
        else
656
          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
657
        endif
658
!         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)
659
      endif
660
661
   else   ! k > INB  au-dessus du nuage
662
    uscavtrac=0.
663
   endif
664
665
! =====    tendances finales  ======
666
     trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k)              ! td transport sans eau dans courants satures
667
     dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k)            ! td du lessivage dans courants satures
668
     dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k)           ! td courant insat
669
     dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k)   ! td courant sat
670
     dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it)         td conv
671
!!!!!!
672
     dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
673
     ENDDO
674
  ENDDO
675
676
! test de conservation du traceur
677
!print*,"_____________________________________________________________"
678
!print*,"                                                             "
679
!      conserv=0.
680
!      conservMA=0.
681
!      DO k= klev-1,1,-1
682
!        DO i=1, klon
683
!         conserv=conserv+dtrcv(i,k,it)*   &
684
!        (paprs(i,k)-paprs(i,k+1))/RG
685
!         conservMA=conservMA+dtrcvMA(i,k,it)*   &
686
!        (paprs(i,k)-paprs(i,k+1))/RG
687
!
688
!      if(it.eq.3)  write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
689
!         'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
690
!         ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,'         conservMA ',conservMA,'conserv ',conserv
691
!!
692
!        ENDDO
693
!      ENDDO
694
!       if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
695
696
END SUBROUTINE cvltr