Ignore:
Timestamp:
Apr 25, 2013, 5:27:27 PM (12 years ago)
Author:
Laurent Fairhead
Message:

Version testing basée sur r1745

Testing release based on r1745

Location:
LMDZ5/branches/testing
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/branches/testing

  • LMDZ5/branches/testing/libf/phylmd/cvltr.F90

    r1279 r1750  
    22! $Id $
    33!
    4 SUBROUTINE cvltr(pdtime,da, phi, mp,paprs,pplay,x,upd,dnd,dx)
     4SUBROUTINE cvltr(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
     5           sigd,sij,clw,elij,epmlmMm,eplaMm,              &
     6           pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM,     &
     7           paprs,it,tr,upd,dnd,inb,icb,                   &
     8           dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
     9           qPa,qMel,qTrdi,dtrcvMA,Mint,                   &
     10           zmfd1a,zmfphi2,zmfdam)
     11  USE IOIPSL
    512  USE dimphy
     13  USE infotrac, ONLY : nbtr,tname
    614  IMPLICIT NONE
    715!=====================================================================
    816! Objet : convection des traceurs / KE
    917! Auteurs: M-A Filiberti and J-Y Grandpeix
     18! modifiee par R Pilon : lessivage des traceurs / KE
    1019!=====================================================================
    1120
    1221  include "YOMCST.h"
    13   include "YOECUMF.h"
     22  include "YOECUMF.h"
     23  include "conema3.h"
    1424
    1525! Entree
    … …  
    1727  REAL,DIMENSION(klon,klev),INTENT(IN)      :: da
    1828  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
    19   REAL,DIMENSION(klon,klev),INTENT(IN)      :: mp
    20   REAL,DIMENSION(klon,klev+1),INTENT(IN)    :: paprs ! pression aux 1/2 couches (bas en haut)
    21   REAL,DIMENSION(klon,klev),INTENT(IN)      :: pplay ! pression pour le milieu de chaque couche
    22   REAL,DIMENSION(klon,klev),INTENT(IN)      :: x     ! q de traceur (bas en haut)
     29! RomP
     30  REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam ! matrices pour simplifier
     31  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2    ! l'ecriture des tendances
     32!
     33  REAL,DIMENSION(klon,klev),INTENT(IN)      :: mpIN
     34  REAL,DIMENSION(klon,klev+1),INTENT(IN)    :: paprs  ! pression aux 1/2 couches (bas en haut)
     35!  REAL,DIMENSION(klon,klev),INTENT(IN)    :: pplay ! pression aux 1/2 couches (bas en haut)
     36  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr     ! q de traceur (bas en haut)
     37  INTEGER,INTENT(IN)                        :: it
    2338  REAL,DIMENSION(klon,klev),INTENT(IN)      :: upd   ! saturated updraft mass flux
    2439  REAL,DIMENSION(klon,klev),INTENT(IN)      :: dnd   ! saturated downdraft mass flux
    25 
     40!
     41  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA   ! masses precipitantes de l'asc adiab
     42  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM   ! masses precipitantes des melanges
     43  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxrIN   ! vprecip: eau
     44  REAL,DIMENSION(klon,klev),INTENT(IN)      :: pmflxsIN   ! vprecip: neige
     45  REAL,DIMENSION(klon,klev),INTENT(IN)      :: ev         ! evaporation cv30_routine
     46  REAL,DIMENSION(klon,klev),INTENT(IN)      :: epIN
     47  REAL,DIMENSION(klon,klev),INTENT(IN)      :: te
     48  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij        ! fraction dair de lenv
     49  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij       ! contenu en eau condensée spécifique/conc deau condensée massique
     50  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm    ! eau condensee precipitee dans mel masse dair sat
     51  REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm    ! eau condensee precipitee dans aa masse dair sat
     52
     53  REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw        ! contenu en eau condensée dans lasc adiab
     54  REAL,DIMENSION(klon),INTENT(IN)           :: sigd
     55  INTEGER,DIMENSION(klon),INTENT(IN)        :: icb,inb
    2656! Sortie
    27   REAL,DIMENSION(klon,klev),INTENT(OUT) :: dx ! tendance de traceur  (bas en haut)
    28 
    29 ! Variables locales     
    30 ! REAL,DIMENSION(klon,klev)       :: zed
    31   REAL,DIMENSION(klon,klev,klev)  :: zmd
    32   REAL,DIMENSION(klon,klev,klev)  :: za
    33   REAL,DIMENSION(klon,klev)       :: zmfd,zmfa
    34   REAL,DIMENSION(klon,klev)       :: zmfp,zmfu
    35   INTEGER                         :: i,k,j
    36   REAL                            :: pdtimeRG
     57  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcv     ! tendance totale (bas en haut)
     58  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrcvMA ! M-A Filiberti
     59  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: trsptd    ! tendance du transport
     60  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrSscav  ! tendance du lessivage courant sat
     61  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrsat    ! tendance trsp+sat scav
     62  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)     :: dtrUscav  ! tendance du lessivage courant unsat
     63!
     64! Variables locales
     65  INTEGER                           :: i,j,k
     66  REAL,DIMENSION(klon,klev)         :: dxpres   ! difference de pression entre niveau (j+1) et (j)
     67  REAL                              :: pdtimeRG ! pas de temps * gravite
     68! variables pour les courants satures
     69  REAL,DIMENSION(klon,klev,klev)    :: zmd
     70  REAL,DIMENSION(klon,klev,klev)    :: za
     71  REAL,DIMENSION(klon,klev,nbtr)    :: zmfd,zmfa
     72  REAL,DIMENSION(klon,klev,nbtr)    :: zmfp,zmfu
     73
     74  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfd1a
     75  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfdam
     76  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: zmfphi2
     77
     78! RomP ! les variables sont nettoyees des valeurs aberrantes
     79  REAL,DIMENSION(klon,klev)         :: Pa, Pm  ! pluie AA et mélanges, var temporaire
     80  REAL,DIMENSION(klon,klev)         :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
     81  REAL,DIMENSION(klon,klev)         :: mp            ! flux de masse
     82  REAL,DIMENSION(klon,klev)         :: ep            ! fraction d'eau convertie en precipitation
     83  REAL,DIMENSION(klon,klev)         :: evap          ! evaporation : variable temporaire
     84  REAL,DIMENSION(klon,klev)         :: rho    !environmental density
     85
     86  REAL,DIMENSION(klon,klev)         :: kappa ! denominateur du au calcul de la matrice
     87                                             ! pour obtenir qd et qp
     88  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qTrdi ! traceurs descente air insature transport MA
     89  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qDi  ! traceurs descente insaturees
     90  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPr  ! traceurs colonne precipitante
     91  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qPa  ! traceurs dans les precip issues lasc. adiab.
     92  REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT)    :: qMel ! traceurs dans les precip issues des melanges
     93  REAL,DIMENSION(klon,klev,nbtr)                :: qMeltmp ! variable temporaire
     94  REAL,DIMENSION(klon,klev,nbtr)                :: qpmMint
     95  REAL,DIMENSION(klon,klev),INTENT(OUT)         :: Mint
     96! tendances
     97  REAL                              :: tdcvMA           ! terme de transport de traceur (schema Marie Angele)
     98  REAL                              :: trsptrac         ! terme de transport de traceur par l'air
     99  REAL                              :: scavtrac         ! terme de lessivage courant sature
     100  REAL                              :: uscavtrac        ! terme de lessivage courant insature
     101! impaction
     102!!!       Correction apres discussion Romain P. / Olivier B.
     103!!!  REAL,PARAMETER                    :: rdrop=2.5e-3     ! rayon des gouttes d'eau
     104  REAL,PARAMETER                    :: rdrop=1.e-3     ! rayon des gouttes d'eau
     105!!!
     106  REAL,DIMENSION(klon,klev)         :: imp              ! coefficient d'impaction
     107! parametres lessivage
     108  REAL                    :: ccntrAA_coef        ! \alpha_a : fract aerosols de l'AA convertis en CCN
     109  REAL                    :: ccntrENV_coef       ! \beta_m  : fract aerosols de l'env convertis en CCN
     110  REAL                    :: coefcoli            ! coefficient de collision des gouttes sur les aerosols
     111!
     112  LOGICAL,DIMENSION(klon,klev) :: NO_precip
     113!  LOGICAL                      :: scavON
     114! var tmp tests
     115  REAL                           :: conserv
     116  real                           :: conservMA
     117
     118! coefficient lessivage
     119   ccntrAA_coef     = 0.
     120   ccntrENV_coef    = 0.
     121   coefcoli         = 0.
     122
     123  call getin('ccntrAA_coef',ccntrAA_coef)
     124  call getin('ccntrENV_coef',ccntrENV_coef)
     125  call getin('coefcoli',coefcoli)
     126  print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
     127
     128!  scavON=.TRUE.
     129!  if(scavON) then
     130!   ccntrAA_coef     = 1.
     131!   ccntrENV_coef    = 1.
     132!   coefcoli         = 1.
     133!  else
     134!   ccntrAA_coef     = 0.
     135!   ccntrENV_coef    = 0.
     136!   coefcoli         = 0.
     137!  endif
     138
     139! ======================================================
     140! calcul de l'impaction
     141! ======================================================
     142!initialisation
     143  do j=1,klev
     144   do i=1,klon
     145     imp(i,j)=0.
     146   enddo
     147  enddo
     148! impaction sur la surface de la colonne de la descente insaturee
     149! On prend la moyenne des precip entre le niveau i+1 et i
     150! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l)
     151!  1000kg/m3= densité de l'eau
     152! 0.75e-3 = 3/4 /1000
     153! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
     154! on le néglige ici pour simplifier le code
     155  do j=1,klev-1
     156   do i=1,klon
     157    imp(i,j) = coefcoli*0.75e-3/rdrop *&
     158             0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
     159!    rho(i,j)=pplay(i,j)/(rd*te(i,j))
     160   enddo
     161  enddo
     162!
     163! initialisation pour flux de traceurs, td et autre
     164   trsptrac = 0.
     165   scavtrac = 0.
     166   uscavtrac = 0.
     167
     168  DO j=1,klev
     169   DO i=1,klon
     170    zmfd(i,j,it)=0.
     171    zmfa(i,j,it)=0.
     172    zmfu(i,j,it)=0.
     173    zmfp(i,j,it)=0.
     174    zmfphi2(i,j,it)=0.
     175    zmfd1a(i,j,it)=0.
     176    zmfdam(i,j,it)=0.
     177    qDi(i,j,it)=0.
     178    qPr(i,j,it)=0.
     179    qPa(i,j,it)=0.
     180    qMel(i,j,it)=0.
     181    qMeltmp(i,j,it)=0.
     182    qTrdi(i,j,it)=0.
     183    kappa(i,j)=0.
     184    trsptd(i,j,it)=0.
     185    dtrsat(i,j,it)=0.
     186    dtrSscav(i,j,it)=0.
     187    dtrUscav(i,j,it)=0.
     188    dtrcv(i,j,it)=0.
     189    dtrcvMA(i,j,it)=0.
     190    evap(i,j)=0.
     191    dxpres(i,j)=0.
     192    qpmMint(i,j,it)=0.
     193    Mint(i,j)=0.
     194   END DO
     195  END DO
     196
     197! suppression des valeurs très faibles (~1e-320)
     198! multiplication de levaporation pour lavoir par unite de temps
     199! et par unite de surface de la maille
     200! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
     201  DO j=1,klev
     202   DO i=1,klon
     203    if(ev(i,j).lt.1.e-16) then
     204     evap(i,j)=0.
     205    else
     206     evap(i,j)=ev(i,j)*sigd(i)
     207    endif
     208   END DO
     209  END DO
     210
     211  DO j=1,klev
     212   DO i=1,klon
     213   if(j.lt.klev) then
     214    if(epIN(i,j).lt.1.e-32) then
     215     ep(i,j)=0.
     216    else
     217     ep(i,j)=epIN(i,j)
     218    endif
     219   else
     220    ep(i,j)=epmax
     221   endif
     222    if(mpIN(i,j).lt.1.e-32) then
     223     mp(i,j)=0.
     224    else
     225     mp(i,j)=mpIN(i,j)
     226    endif
     227    if(pmflxsIN(i,j).lt.1.e-32) then
     228     pmflxs(i,j)=0.
     229    else
     230     pmflxs(i,j)=pmflxsIN(i,j)
     231    endif
     232    if(pmflxrIN(i,j).lt.1.e-32) then
     233     pmflxr(i,j)=0.
     234    else
     235     pmflxr(i,j)=pmflxrIN(i,j)
     236    endif
     237    if(wdtrainA(i,j).lt.1.e-32) then
     238     Pa(i,j)=0.
     239    else
     240     Pa(i,j)=wdtrainA(i,j)
     241    endif
     242    if(wdtrainM(i,j).lt.1.e-32) then
     243     Pm(i,j)=0.
     244    else
     245     Pm(i,j)=wdtrainM(i,j)
     246    endif
     247   END DO
     248  END DO
     249
     250!==========================================
     251  DO j = klev-1,1,-1
     252   DO i = 1,klon
     253     NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
     254                    .and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
     255   END DO
     256  END DO
    37257
    38258! =========================================
    … …  
    40260! =========================================
    41261!cdir collapse
    42   DO j=1,klev
    43   DO i=1,klon
    44 !   zed(i,j)=0.
    45     zmfd(i,j)=0.
    46     zmfa(i,j)=0.
    47     zmfu(i,j)=0.
    48     zmfp(i,j)=0.
    49   END DO
    50   END DO
    51 !cdir collapse
    52262  DO k=1,klev
    53   DO j=1,klev
    54   DO i=1,klon
    55     zmd(i,j,k)=0.
    56     za (i,j,k)=0.
    57   END DO
    58   END DO
    59   END DO
    60 ! entrainement
    61 ! DO k=1,klev-1
    62 !    DO i=1,klon
    63 !       zed(i,k)=max(0.,mp(i,k)-mp(i,k+1))
    64 !    END DO
    65 ! END DO
    66 
     263   DO j=1,klev
     264    DO i=1,klon
     265     zmd(i,j,k)=0.
     266     za (i,j,k)=0.
     267    END DO
     268   END DO
     269  END DO
    67270! calcul de la matrice d echange
    68271! matrice de distribution de la masse entrainee en k
    69 
     272! commmentaire RomP : mp > 0
    70273  DO k=1,klev-1
    71274     DO i=1,klon
    72         zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1))
     275        zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1))   ! ~ mk(k)
    73276     END DO
    74277  END DO
    … …  
    76279     DO j=k-1,1,-1
    77280        DO i=1,klon
    78            if(mp(i,j+1).ne.0) then
    79               zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1))
     281           if(mp(i,j+1).gt.1.e-10) then
     282              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)
    80283           ENDif
    81284        END DO
    … …  
    89292     END DO
    90293  END DO
     294!!!!! quantite  de traceur dans la descente d'air insaturee  :   4 juin 2012
     295  DO k=1,klev
     296     DO j=1,klev-1
     297        DO i=1,klon
     298        if(mp(i,j+1).gt.1.e-10) then
     299          qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
     300        else
     301          qTrdi(i,j,it)=0.!tr(i,j,it)
     302        endif
     303        ENDDO
     304     ENDDO
     305  ENDDO
     306!!!!!
    91307!
    92308! rajout du terme lie a l ascendance induite
    … …  
    98314  END DO
    99315!
    100 ! tendances
    101 !           
     316! tendance courants insatures  ! sans lessivage ancien schema
     317!
    102318  DO k=1,klev
    103319     DO j=1,klev
    104320        DO i=1,klon
    105            zmfd(i,j)=zmfd(i,j)+za(i,j,k)*(x(i,k)-x(i,j))
     321           zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
    106322        END DO
    107323     END DO
    … …  
    109325!
    110326! =========================================
    111 ! calcul des tendances liees aux flux satures
     327! calcul des tendances liees aux courants satures   j <-> z ; k <-> z'
    112328! =========================================
    113329  DO j=1,klev
    114330     DO i=1,klon
    115         zmfa(i,j)=da(i,j)*(x(i,1)-x(i,j))
     331        zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it))                     ! da
    116332     END DO
    117333  END DO
    … …  
    119335     DO j=1,klev
    120336        DO i=1,klon
    121            zmfp(i,j)=zmfp(i,j)+phi(i,j,k)*(x(i,k)-x(i,j))
     337           zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it))  ! phi
     338        END DO
     339     END DO
     340  END DO
     341! RomP ajout des matrices liees au lessivage
     342  DO j=1,klev
     343     DO i=1,klon
     344        zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it)   ! da1
     345        zmfdam(i,j,it)=dam(i,j)*tr(i,1,it)   ! dam
     346     END DO
     347  END DO
     348  DO k=1,klev
     349     DO j=1,klev
     350        DO i=1,klon
     351          zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it)  ! psi
    122352        END DO
    123353     END DO
    … …  
    125355  DO j=1,klev-1
    126356     DO i=1,klon
    127         zmfu(i,j)=max(0.,upd(i,j+1)+dnd(i,j+1))*(x(i,j+1)-x(i,j))
     357        zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
    128358     END DO
    129359  END DO
    130360  DO j=2,klev
    131361     DO i=1,klon
    132         zmfu(i,j)=zmfu(i,j)+min(0.,upd(i,j)+dnd(i,j))*(x(i,j)-x(i,j-1))
    133      END DO
    134   END DO
    135 
    136 ! =========================================
    137 ! calcul final des tendances
    138 ! =========================================
     362        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))
     363     END DO
     364  END DO
     365! ===================================================
     366! calcul des tendances liees aux courants insatures
     367! ===================================================
     368!  pression 
    139369  DO k=1, klev
    140370     DO i=1, klon
    141         dx(i,k)=paprs(i,k)-paprs(i,k+1)
     371        dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
    142372     ENDDO
    143373  ENDDO
    144374  pdtimeRG=pdtime*RG
    145 !cdir collapse
    146   DO k=1, klev
    147      DO i=1, klon
    148         dx(i,k)=(zmfd(i,k)+zmfu(i,k)       &
    149                 +zmfa(i,k)+zmfp(i,k))*pdtimeRG/dx(i,k)
    150         !          print*,'dx',k,dx(i,k)
     375
     376! q_pa et q_pm    traceurs issues des courants satures se retrouvant dans les precipitations
     377  DO j=1,klev
     378     DO i=1,klon
     379        if(j.ge.icb(i).and.j.le.inb(i)) then
     380          if(clw(i,j).gt.1.e-16) then
     381           qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
     382          else
     383           qPa(i,j,it)=0.
     384          endif
     385        endif
     386     END DO
     387  END DO
     388
     389! calcul de q_pm en 2 parties :
     390! 1) calcul de sa valeur pour un niveau z' donne
     391! 2) integration sur la verticale sur z'
     392     DO j=1,klev
     393      DO k=1,j-1
     394        DO i=1,klon
     395        if(k.ge.icb(i).and.k.le.inb(i).and.&
     396           j.le.inb(i)) then
     397            if(elij(i,k,j).gt.1.e-16) then
     398             qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)&
     399                         *(1.-sij(i,k,j))  +ccntrENV_coef&
     400                         *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j)
     401            else
     402             qMeltmp(i,j,it)=0.
     403            endif
     404          qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
     405          Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
     406        endif ! end if dans nuage
     407        END DO
     408     END DO
     409  END DO
     410
     411     DO j=1,klev
     412        DO i=1,klon
     413          if(Mint(i,j).gt.1.e-16) then
     414           qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
     415          else
     416           qMel(i,j,it)=0.
     417          endif
     418     END DO
     419  END DO
     420
     421! calcul de q_d et q_p    traceurs de la descente precipitante
     422   DO j=klev-1,1,-1
     423    DO i=1,klon
     424     if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then  ! detrainement
     425       kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
     426                (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
     427                + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
     428             
     429     elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
     430       if(j.eq.1) then
     431        kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
     432                 (-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
     433                 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
     434       else
     435        kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
     436                 (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
     437                 + (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j)))
     438       endif
     439      else
     440        kappa(i,j)=1.
     441      endif
     442    ENDDO
     443   ENDDO
     444
     445  DO j=klev-1,1,-1
     446   DO i=1,klon
     447    if (abs(kappa(i,j)).lt.1.e-25) then    !si denominateur nul (il peut y avoir des mp!=0)
     448     kappa(i,j)=1.
     449     if(j.eq.1) then
     450       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
     451     elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then
     452       qDi(i,j,it)=qDi(i,j+1,it)
     453     elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
     454       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))
     455     else  ! si  mp (i)=0 et mp(j+1)=0
     456       qDi(i,j,it)=tr(i,j,it) ! orig 0.
     457     endif
     458
     459      if(NO_precip(i,j)) then
     460       qPr(i,j,it)=0.
     461      else
     462       qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     463                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
     464                   +imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/&
     465                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
     466      endif
     467    else   !     denominateur non nul
     468     kappa(i,j)=1./kappa(i,j)     
     469! calcul de qd et qp
     470!!jyg  (20130119) correction pour le sommet du nuage
     471!!     if(j.ge.inb(i)) then       !au-dessus du nuage, sommet inclu
     472     if(j.gt.inb(i)) then       !au-dessus du nuage
     473       qDi(i,j,it)=tr(i,j,it)   ! pas de descente => environnement = descente insaturee
     474       qPr(i,j,it)=0.
     475
     476!  vvv premiere couche du modele ou mp(1)=0  ! det tout le temps  vvv
     477     elseif(j.eq.1) then
     478      if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
     479       if(NO_precip(i,j)) then !pas de precip en (i)
     480        qDi(i,j,it)=qDi(i,j+1,it)
     481        qPr(i,j,it)=0.
     482       else
     483        qDi(i,j,it)=kappa(i,j)*(&
     484            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     485            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
     486            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
     487            (-mp(i,j+1)*qDi(i,j+1,it)))
     488
     489        qPr(i,j,it)=kappa(i,j)*(&
     490            (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
     491            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     492            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
     493            +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
     494       endif
     495
     496      else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
     497        qDi(i,j,it)=tr(i,j,it) ! orig 0.
     498        if(NO_precip(i,j)) then
     499         qPr(i,j,it)=0.
     500        else
     501         qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     502                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
     503                   +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
     504                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
     505        endif
     506
     507      endif  !mp(2) nul ou non
     508
     509! vvv  (j!=1.and.j.lt.inb(i))  en-dessous du sommet nuage   vvv
     510     else
     511!------------------------------------------------------------- detrainement
     512      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
     513         if(NO_precip(i,j)) then
     514          qDi(i,j,it)=qDi(i,j+1,it)
     515          qPr(i,j,it)=0.
     516         else
     517          qDi(i,j,it)=kappa(i,j)*(&
     518            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     519            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
     520            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
     521            (-mp(i,j+1)*qDi(i,j+1,it)))
     522!
     523          qPr(i,j,it)=kappa(i,j)*(&
     524            (-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*&
     525            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     526            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
     527            +(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j)))
     528         endif !precip
     529!------------------------------------------------------------- entrainement
     530      elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then
     531         if(NO_precip(i,j)) then
     532          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))
     533          qPr(i,j,it)=0.
     534         else
     535          qDi(i,j,it)=kappa(i,j)*(&
     536            (-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     537            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +&
     538            (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
     539            (-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it)))
     540!
     541          qPr(i,j,it)=kappa(i,j)*(&
     542            (-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*&
     543            ((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     544            Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))&
     545            +(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*&
     546            (imp(i,j)/RG*dxpres(i,j)))
     547         endif !precip
     548!------------------------------------------------------------- endif ! ent/det
     549      else  !mp nul
     550        qDi(i,j,it)=tr(i,j,it) ! orig 0.
     551        if(NO_precip(i,j)) then
     552         qPr(i,j,it)=0.
     553        else
     554         qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
     555                   Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)&
     556                   +imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/&
     557                  (pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
     558        endif
     559      endif ! mp nul ou non
     560     endif ! condition sur j
     561    endif ! kappa
     562   ENDDO
     563  ENDDO
     564
     565!! print test descente insaturee
     566!  DO j=klev,1,-1
     567!   DO i=1,klon
     568!     if(it.eq.3) then
     569!   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,&
     570!!         'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
     571!          'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
     572!          'Pm',Pm(i,j),'Mint',Mint(i,j),&
     573!!          'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
     574!        'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
     575!!          'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
     576!!         'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
     577!     endif
     578!   ENDDO
     579!  ENDDO
     580
     581
     582! ===================================================
     583! calcul final des tendances
     584! ===================================================
     585
     586  DO k=klev-1,1,-1
     587    DO i=1, klon
     588! transport
     589     tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)   ! double comptage des downdraft insatures
     590     trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
     591! lessivage courants satures
     592     scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
     593               -zmfphi2(i,k,it)*ccntrENV_coef&
     594               -zmfdam(i,k,it)*ccntrAA_coef
     595! lessivage courants insatures
     596   if(k.le.inb(i).and.k.gt.1) then   ! tendances dans le nuage
     597!------------------------------------------------------------- detrainement
     598      if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
     599        uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
     600                   + mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
     601!
     602!        if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
     603!                    (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
     604!                    mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
     605!                    'mp',mp(i,k)
     606!------------------------------------------------------------- entrainement
     607      elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
     608         uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
     609!
     610!         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)
     611!=!------------------------------------------------------------- end ent/det
     612      else !        mp(i,k+1)=0. et mp(i,k)=0.    pluie directement sur l environnement
     613
     614        if(NO_precip(i,k)) then
     615          uscavtrac=0.
     616!         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)
     617        else
     618          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
     619!         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)
     620        endif
     621      endif  ! mp/det/ent
     622!------------------------------------------------------------- premiere couche
     623   elseif(k.eq.1) then  !                                      mp(1)=0.
     624      if(mp(i,2).gt.1.e-10) then  !detrainement
     625         uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
     626!                   + mp(i,2)*(0.-tr(i,k,it))
     627!
     628!       if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
     629!                  (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+&
     630!                  mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),&
     631!                   'mp',mp(i,k)
     632      else   ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
     633        if(NO_precip(i,1)) then
     634          uscavtrac=0.
     635        else
     636          uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG
     637        endif
     638!         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)
     639      endif
     640
     641   else   ! k > INB  au-dessus du nuage
     642    uscavtrac=0.
     643   endif
     644
     645! =====    tendances finales  ======
     646     trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k)              ! td transport sans eau dans courants satures
     647     dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k)            ! td du lessivage dans courants satures
     648     dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k)           ! td courant insat
     649     dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k)   ! td courant sat
     650     dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it)         td conv
     651!!!!!!
     652     dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
    151653     ENDDO
    152654  ENDDO
    153655
    154656! test de conservation du traceur
     657!print*,"_____________________________________________________________"
     658!print*,"                                                             "
    155659!      conserv=0.
    156 !      DO k=1, klev
     660!      conservMA=0.
     661!      DO k= klev-1,1,-1
    157662!        DO i=1, klon
    158 !         conserv=conserv+dx(i,k)*   &
     663!         conserv=conserv+dtrcv(i,k,it)*   &
    159664!        (paprs(i,k)-paprs(i,k+1))/RG
     665!         conservMA=conservMA+dtrcvMA(i,k,it)*   &
     666!        (paprs(i,k)-paprs(i,k+1))/RG
     667!
     668!      if(it.eq.3)  write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
     669!         'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
     670!         ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,'         conservMA ',conservMA,'conserv ',conserv
     671!!
    160672!        ENDDO
    161673!      ENDDO
    162 !      print *,'conserv',conserv
    163      
     674!       if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
     675
    164676END SUBROUTINE cvltr
Note: See TracChangeset for help on using the changeset viewer.