Changeset 1742


Ignore:
Timestamp:
Apr 5, 2013, 1:49:35 PM (11 years ago)
Author:
idelkadi
Message:

1- Inclusion des developpements de la these de Romain Pilon sur le
lessivage des aerosols :

a/ par les pluies convectives (modifs cv30_routines et cv3_routines pour

sortir les champs nécessaires au calcul off-line ; modif cvltr)

b/ par les pluies stratiformes (modifs phytrac et introduction

lsc_scav).

2- Choix entre plusieurs schemas pour les pluies stratiformes, commande
par iflag_lscav.

3- Quelques corrections dans la convection "Nouvelle Physique" pour
assurer la conservation des traceurs (cv3p1_mixing et cva_driver) (travail
de Robin Locatelli).

Location:
LMDZ5/trunk/libf/phylmd
Files:
2 added
22 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cltrac.F90

    r1279 r1742  
    22! $Id $
    33!
    4 SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp,d_tr)
     4SUBROUTINE cltrac(dtime,coef,t,tr,flux,paprs,pplay,delp, &
     5                  d_tr,d_tr_dry,flux_tr_dry)                    !jyg
     6
    57  USE dimphy
    68  IMPLICIT NONE
    … …  
    1315! Arguments:
    1416!-----------
    15 ! dtime....input-R- intervalle du temps (en secondes)
    16 ! coef.....input-R- le coefficient d'echange (m**2/s) l>1
    17 ! t........input-R- temperature (K)
    18 ! tr.......input-R- la q. de traceurs
    19 ! flux.....input-R- le flux de traceurs a la surface
    20 ! paprs....input-R- pression a inter-couche (Pa)
    21 ! pplay....input-R- pression au milieu de couche (Pa)
    22 ! delp.....input-R- epaisseur de couche (Pa)
    23 ! cdrag....input-R- cdrag pour le flux de surface (non active)
    24 ! tr0......input-R- traceurs a la surface ou dans l'ocean (non active)
    25 ! d_tr.....output-R- le changement de tr
    26 ! flux_tr..output-R- flux de tr
     17! dtime.......input-R- intervalle du temps (en secondes)
     18! coef........input-R- le coefficient d'echange (m**2/s) l>1
     19! t...........input-R- temperature (K)
     20! tr..........input-R- la q. de traceurs
     21! flux........input-R- le flux de traceurs a la surface
     22! paprs.......input-R- pression a inter-couche (Pa)
     23! pplay.......input-R- pression au milieu de couche (Pa)
     24! delp........input-R- epaisseur de couche (Pa)
     25! cdrag.......input-R- cdrag pour le flux de surface (non active)
     26! tr0.........input-R- traceurs a la surface ou dans l'ocean (non active)
     27! d_tr........output-R- le changement de tr
     28! d_tr_dry....output-R- le changement de tr du au depot sec (1st layer)
     29! flux_tr_dry.output-R- depot sec
     30!!! flux_tr..output-R- flux de tr
    2731!======================================================================
    2832  include "YOMCST.h"
    … …  
    4044!
    4145  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: d_tr
     46  REAL ,DIMENSION(klon),INTENT(OUT)       :: d_tr_dry          !jyg
     47  REAL ,DIMENSION(klon),INTENT(OUT)       :: flux_tr_dry       !jyg
    4248!  REAL ,DIMENSION(klon,klev),INTENT(OUT) :: flux_tr
    4349!
    … …  
    6672     zx_alf1(i) = (paprs(i,1)-pplay(i,2))/(pplay(i,1)-pplay(i,2))
    6773     zx_alf2(i) = 1.0 - zx_alf1(i)
    68      zx_flux(i) =  -flux(i)*dtime*RG
     74     flux_tr_dry(i) = -flux(i)*dtime                              !jyg
     75     zx_flux(i) =  flux_tr_dry(i)*RG                              !jyg
     76!!     zx_flux(i) =  -flux(i)*dtime*RG                            !jyg
    6977! Pour le moment le flux est prescrit cdrag et zx_coef(1) vaut 0
    7078     cdrag(i) = 0.0
    … …  
    95103     zx_dtr(i,2) = (zx_coef(i,2)-zx_alf2(i)*zx_coef(i,1)) /   &
    96104          zx_buf(i)
     105     d_tr_dry(i) = -zx_flux(i)/zx_buf(i)                          !jyg
    97106  ENDDO
    98107

  • LMDZ5/trunk/libf/phylmd/concvl.F

    r1650 r1742  
    1010     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
    1111     .             qcondc,wd,pmflxr,pmflxs,
    12      .             da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
     12! RomP >>>
     13!!     .             da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
     14     .             da,phi,mp,phi2,d1a,dam,sij,clw,elij,        ! RomP
     15     .             dd_t,dd_q,lalim_conv,wght_th,               ! RomP
     16     .             evap, ep, epmlmMm,eplaMm,                   ! RomP
     17     .             wdtrainA,wdtrainM)                          ! RomP
     18! RomP <<<
    1319***************************************************************
    1420*                                                             *
    … …  
    9197
    9298       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
     99! RomP >>>
     100       real phi2(klon,klev,klev)                                   
     101       real d1a(klon,klev),dam(klon,klev)
     102       real sij(klon,klev,klev),clw(klon,klev),elij(klon,klev,klev)
     103       REAL wdtrainA(klon,klev),wdtrainM(klon,klev)
     104       REAL evap(klon,klev),ep(klon,klev)
     105       REAL epmlmMm(klon,klev,klev),eplaMm(klon,klev)
     106! RomP <<<
    93107       REAL cape(klon),cin(klon),tvp(klon,klev)
    94108       REAL Tconv(klon,klev)
    … …  
    248262         DO i = 1, klon
    249263          cbmf(i) = 0.
    250           plcl(i) = 0.
     264!!          plcl(i) = 0.
    251265          sigd(i) = 0.
    252266         ENDDO
    … …  
    256270      plfc(:)  = 0.
    257271      wbeff(:) = 100.
     272      plcl(:) = 0.
    258273
    259274      DO k = 1, klev+1
    … …  
    339354      if (iflag_con.eq.30) then
    340355
    341       CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
     356      print *, '-> cv_driver'      !jyg
     357      CALL cv_driver(klon,klev,klevp1,ntra,iflag_con,
    342358     :              t,q,qs,u,v,tra,
    343359     $              em_p,em_ph,iflag,
    344360     $              d_t,d_q,d_u,d_v,d_tra,rain,
    345 !!     $              pmflxr,cbmf,work1,work2,           !jyg
    346      $              Vprecip,cbmf,work1,work2,            !jyg
     361     $              Vprecip,cbmf,work1,work2,                  !jyg
    347362     $              kbas,ktop,
    348363     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
    349      $              da,phi,mp)
     364     $              da,phi,mp,phi2,d1a,dam,sij,clw,elij,       !RomP
     365     $              evap,ep,epmlmMm,eplaMm,                    !RomP
     366     $              wdtrainA,wdtrainM)                         !RomP
     367      print *, 'cv_driver ->'      !jyg
    350368c
    351369      DO i = 1,klon
    … …  
    369387     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
    370388     $              asupmaxmin,lalim_conv,
    371 !AC!
    372      $              da,phi)
    373 !AC!
     389!AC!+!RomP
     390     $              da,phi,mp,phi2,d1a,dam,sij,clw,  ! RomP
     391     $              elij,evap,ep,wdtrainA,wdtrainM)  ! RomP
     392!AC!+!RomP
    374393      endif 
    375394C------------------------------------------------------------------

  • LMDZ5/trunk/libf/phylmd/cv30_routines.F

    r1403 r1742  
    18311831     :              ,th,tv,lv,cpn,ep,sigp,clw
    18321832     :              ,m,ment,elij,delt,plcl
    1833      :              ,mp,rp,up,vp,trap,wt,water,evap,b)
     1833     :              ,mp,rp,up,vp,trap,wt,water,evap,b ! RomP-jyg
     1834     :              ,wdtrainA,wdtrainM)               ! 26/08/10  RomP-jyg
    18341835      implicit none
    18351836
    … …  
    18571858      real trap(nloc,na,ntra)
    18581859      real b(nloc,na)
     1860! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
     1861! lascendance adiabatique et des flux melanges Pa et Pm.
     1862! Distinction des wdtrain
     1863! Pa = wdtrainA     Pm = wdtrainM
     1864      real  wdtrainA(nloc,na), wdtrainM(nloc,na)
    18591865
    18601866c local variables
    … …  
    18981904c         enddo
    18991905c        enddo
    1900 
     1906!! RomP >>>
     1907         do i=1,nd
     1908          do il=1,ncum
     1909          wdtrainA(il,i)=0.0     
     1910          wdtrainM(il,i)=0.0     
     1911         enddo
     1912        enddo
     1913!! RomP <<<
    19011914c
    19021915c   ***  check whether ep(inb)=0, if so, skip precipitating    ***
    … …  
    19351948         if (cvflag_grav) then
    19361949          wdtrain(il)=grav*ep(il,i)*m(il,i)*clw(il,i)
     1950          wdtrainA(il,i) = wdtrain(il)/grav     !   Pa  26/08/10   RomP
    19371951         else
    19381952          wdtrain(il)=10.0*ep(il,i)*m(il,i)*clw(il,i)
     1953          wdtrainA(il,i) = wdtrain(il)/10.      !   Pa  26/08/10   RomP
    19391954         endif
    19401955        endif
    … …  
    19421957
    19431958       if(i.gt.1)then
     1959
    19441960        do 320 j=1,i-1
    19451961         do il=1,ncum
    … …  
    19551971         enddo
    19561972320     continue
     1973         do il=1,ncum
     1974           if (cvflag_grav) then
     1975           wdtrainM(il,i) = wdtrain(il)/grav-wdtrainA(il,i)    !   Pm  26/08/10   RomP
     1976           else
     1977           wdtrainM(il,i) = wdtrain(il)/10.-wdtrainA(il,i)    !   Pm  26/08/10   RomP
     1978           endif
     1979         enddo
     1980
    19571981       endif
    19581982
    … …  
    30223046        end
    30233047
     3048!!RomP >>>
    30243049      SUBROUTINE cv30_tracer(nloc,len,ncum,nd,na,
    3025      &                        ment,sij,da,phi)
     3050     &                       ment,sij,da,phi,phi2,d1a,dam,
     3051     &                       ep,VPrecip,elij,clw,epmlmMm,eplaMm,
     3052     &                       icb,inb)
    30263053        implicit none
     3054
     3055#include "cv30param.h"
     3056
    30273057c inputs:
    30283058        integer ncum, nd, na, nloc,len
    30293059        real ment(nloc,na,na),sij(nloc,na,na)
     3060        real clw(nloc,nd),elij(nloc,na,na)
     3061        real ep(nloc,na)
     3062        integer icb(nloc),inb(nloc)
     3063        real VPrecip(nloc,nd+1)
    30303064c ouputs:
    30313065        real da(nloc,na),phi(nloc,na,na)
     3066        real phi2(nloc,na,na)
     3067        real d1a(nloc,na),dam(nloc,na)
     3068        real epmlmMm(nloc,na,na),eplaMm(nloc,na)
     3069! variables pour tracer dans precip de l'AA et des mel
    30323070c local variables:
    30333071        integer i,j,k
    3034 c       
    3035         da(:,:)=0.
    3036 c
     3072        real epm(nloc,na,na)
     3073c
     3074! variables d'Emanuel : du second indice au troisieme
     3075! --->    tab(i,k,j) -> de l origine k a l arrivee j
     3076!  ment, sij, elij
     3077! variables personnelles : du troisieme au second indice
     3078! --->    tab(i,j,k) -> de k a j
     3079! phi, phi2
     3080!
     3081! initialisations
     3082       do j=1,na
     3083        do i=1,ncum
     3084         da(i,j)=0.
     3085         d1a(i,j)=0.
     3086         dam(i,j)=0.
     3087         eplaMm(i,j)=0.
     3088        enddo
     3089       enddo
     3090      do k=1,na
     3091       do j=1,na
     3092        do i=1,ncum
     3093         epm(i,j,k)=0.
     3094         epmlmMm(i,j,k)=0.
     3095         phi(i,j,k)=0.
     3096         phi2(i,j,k)=0.
     3097        enddo
     3098       enddo
     3099      enddo
     3100c
     3101!  fraction deau condensee dans les melanges convertie en precip : epm
     3102! et eau condensée précipitée dans masse d'air saturé : l_m*dM_m/dzdz.dzdz
     3103        do j=1,na
     3104         do k=1,j-1
     3105           do i=1,ncum
     3106            if(k.ge.icb(i).and.k.le.inb(i).and.
     3107     &         j.le.inb(i)) then
     3108             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
     3109             epm(i,j,k)=max(epm(i,j,k),0.0)
     3110            endif
     3111           end do
     3112         end do
     3113        end do
     3114!
     3115        do j=1,na
     3116         do k=1,na
     3117           do i=1,ncum
     3118            if(k.ge.icb(i).and.k.le.inb(i)) then
     3119             eplaMm(i,j)=eplaMm(i,j) + ep(i,j)*clw(i,j)
     3120     &                  *ment(i,j,k)*(1.-sij(i,j,k))
     3121            endif
     3122           end do
     3123         end do
     3124        end do
     3125!
     3126        do j=1,na
     3127         do k=1,j-1
     3128           do i=1,ncum
     3129            if(k.ge.icb(i).and.k.le.inb(i).and.
     3130     &         j.le.inb(i)) then
     3131             epmlmMm(i,j,k)=epm(i,j,k)*elij(i,k,j)*ment(i,k,j)
     3132            endif
     3133           end do
     3134         end do
     3135        end do
     3136
     3137!  matrices pour calculer la tendance des concentrations dans cvltr.F90
    30373138        do j=1,na
    30383139          do k=1,na
    30393140            do i=1,ncum
    3040             da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
    3041             phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
    3042 c            print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
     3141             da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
     3142             phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
     3143             d1a(i,j)=d1a(i,j)+ment(i,k,j)*ep(i,k)
     3144     &              *(1.-sij(i,k,j))
    30433145            end do
    30443146          end do
    30453147        end do
    3046    
     3148
     3149        do j=1,na
     3150          do k=1,j-1
     3151            do i=1,ncum
     3152             dam(i,j)=dam(i,j)+ment(i,k,j)
     3153     &             *epm(i,j,k)*(1.-ep(i,k))*(1.-sij(i,k,j))
     3154             phi2(i,j,k)=phi(i,j,k)*epm(i,j,k)
     3155            end do
     3156          end do
     3157        end do
     3158
    30473159        return
    30483160        end
    3049 
     3161!RomP <<<
    30503162
    30513163      SUBROUTINE cv30_uncompress(nloc,len,ncum,nd,ntra,idcum
    30523164     :         ,iflag
    3053      :         ,precip,VPrecip,sig,w0
     3165     :         ,precip,VPrecip,evap,ep,sig,w0
    30543166     :         ,ft,fq,fu,fv,ftra
    30553167     :         ,inb
    30563168     :         ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
    3057      :         ,da,phi,mp
     3169     :         ,da,phi,mp,phi2,d1a,dam,sij
     3170     :         ,elij,clw,epmlmMm,eplaMm
     3171     :         ,wdtrainA,wdtrainM
    30583172     :         ,iflag1
    3059      :         ,precip1,VPrecip1,sig1,w01
     3173     :         ,precip1,VPrecip1,evap1,ep1,sig1,w01
    30603174     :         ,ft1,fq1,fu1,fv1,ftra1
    30613175     :         ,inb1
    30623176     :         ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
    3063      :         ,da1,phi1,mp1)
     3177     :         ,da1,phi1,mp1,phi21,d1a1,dam1,sij1
     3178     :         ,elij1,clw1,epmlmMm1,eplaMm1
     3179     :         ,wdtrainA1,wdtrainM1)
    30643180      implicit none
    30653181
    … …  
    30723188      integer inb(nloc)
    30733189      real precip(nloc)
    3074       real VPrecip(nloc,nd+1)
     3190      real VPrecip(nloc,nd+1),evap(nloc,nd)
     3191      real ep(nloc,nd)
    30753192      real sig(nloc,nd), w0(nloc,nd)
    30763193      real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
    … …  
    30813198      real wd(nloc),cape(nloc)
    30823199      real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
     3200!RomP >>>
     3201      real phi2(nloc,nd,nd)
     3202      real d1a(nloc,nd),dam(nloc,nd)
     3203      real wdtrainA(nloc,nd), wdtrainM(nloc,nd)
     3204      real sij(nloc,nd,nd)
     3205      real elij(nloc,nd,nd),clw(nloc,nd)
     3206      real epmlmMm(nloc,nd,nd),eplaMm(nloc,nd)
     3207!RomP <<<
    30833208
    30843209c outputs:
    … …  
    30863211      integer inb1(len)
    30873212      real precip1(len)
    3088       real VPrecip1(len,nd+1)
     3213      real VPrecip1(len,nd+1),evap1(len,nd) !<<< RomP
     3214      real ep1(len,nd)                      !<<< RomP
    30893215      real sig1(len,nd), w01(len,nd)
    30903216      real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
    … …  
    30953221      real wd1(nloc),cape1(nloc)
    30963222      real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
     3223!RomP >>>
     3224      real phi21(len,nd,nd)
     3225      real d1a1(len,nd),dam1(len,nd)
     3226      real wdtrainA1(len,nd), wdtrainM1(len,nd)
     3227      real sij1(len,nd,nd)
     3228      real elij1(len,nd,nd),clw1(len,nd)
     3229      real epmlmMm1(len,nd,nd),eplaMm1(len,nd)
     3230!RomP <<<
    30973231
    30983232c local variables:
    … …  
    31103244          do 2010 i=1,ncum
    31113245            VPrecip1(idcum(i),k)=VPrecip(i,k)
     3246            evap1(idcum(i),k)=evap(i,k)   !<<< RomP
    31123247            sig1(idcum(i),k)=sig(i,k)
    31133248            w01(idcum(i),k)=w0(i,k)
    … …  
    31233258            da1(idcum(i),k)=da(i,k)
    31243259            mp1(idcum(i),k)=mp(i,k)
     3260!RomP >>>
     3261            ep1(idcum(i),k)=ep(i,k)
     3262            d1a1(idcum(i),k)=d1a(i,k)
     3263            dam1(idcum(i),k)=dam(i,k)
     3264            clw1(idcum(i),k)=clw(i,k)
     3265            eplaMm1(idcum(i),k)=eplaMm(i,k)
     3266            wdtrainA1(idcum(i),k)=wdtrainA(i,k)
     3267            wdtrainM1(idcum(i),k)=wdtrainM(i,k)
     3268!RomP <<<
    31253269 2010     continue
    31263270 2020   continue
    … …  
    31413285         do k=1,nd
    31423286          do i=1,ncum
     3287            sij1(idcum(i),k,j)=sij(i,k,j)
    31433288            phi1(idcum(i),k,j)=phi(i,k,j)
     3289            phi21(idcum(i),k,j)=phi2(i,k,j)
     3290            elij1(idcum(i),k,j)=elij(i,k,j)
     3291            epmlmMm1(idcum(i),k,j)=epmlmMm(i,k,j)
    31443292          end do
    31453293         end do

  • LMDZ5/trunk/libf/phylmd/cv3_routines.F

    r1650 r1742  
    19501950     :              ,th,tv,lv,cpn,ep,sigp,clw
    19511951     :              ,m,ment,elij,delt,plcl,coef_clos
    1952      o              ,mp,rp,up,vp,trap,wt,water,evap,b,sigd)
     1952     o              ,mp,rp,up,vp,trap,wt,water,evap,b,sigd
     1953     o              ,wdtrainA,wdtrainM)                                ! RomP
    19531954      implicit none
    19541955
    … …  
    19791980      real trap(nloc,na,ntra)
    19801981      real b(nloc,na), sigd(nloc)
     1982! 25/08/10 - RomP---- ajout des masses precipitantes ejectees
     1983! lascendance adiabatique et des flux melanges Pa et Pm.
     1984! Distinction des wdtrain
     1985! Pa = wdtrainA     Pm = wdtrainM
     1986      real  wdtrainA(nloc,na), wdtrainM(nloc,na)
    19811987
    19821988c local variables
    … …  
    20212027!AC!         enddo
    20222028!AC!        enddo
     2029!! RomP >>>
     2030         do i=1,nd
     2031          do il=1,ncum
     2032          wdtrainA(il,i)=0.0     
     2033          wdtrainM(il,i)=0.0     
     2034         enddo
     2035        enddo
     2036!! RomP <<<
    20232037c
    20242038c   ***  check whether ep(inb)=0, if so, skip precipitating    ***
    … …  
    20652079         if (cvflag_grav) then
    20662080          wdtrain(il)=grav*ep(il,i)*m(il,i)*clw(il,i)
     2081          wdtrainA(il,i) = wdtrain(il)/grav     !   Pa   RomP
    20672082         else
    20682083          wdtrain(il)=10.0*ep(il,i)*m(il,i)*clw(il,i)
     2084          wdtrainA(il,i) = wdtrain(il)/10.      !   Pa   RomP
    20692085         endif
    20702086        endif
    … …  
    20792095           if (cvflag_grav) then
    20802096            wdtrain(il)=wdtrain(il)+grav*awat*ment(il,j,i)
     2097            wdtrainM(il,i) = wdtrain(il)/grav-wdtrainA(il,i)    !   Pm  RomP
    20812098           else
    20822099            wdtrain(il)=wdtrain(il)+10.0*awat*ment(il,j,i)
     2100            wdtrainM(il,i) = wdtrain(il)/10.-wdtrainA(il,i)    !   Pm  RomP
    20832101           endif
    20842102          endif
    … …  
    35403558        end
    35413559
    3542 !AC!
     3560!AC! et !RomP >>>
    35433561      SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
    3544      &                        ment,sij,da,phi)
     3562     &                       ment,sigij,da,phi,phi2,d1a,dam,
     3563     &                       ep,Vprecip,elij,clw,icb,inb)
    35453564        implicit none
     3565
     3566#include "cv3param.h"
     3567
    35463568c inputs:
    35473569        integer ncum, nd, na, nloc,len
    3548         real ment(nloc,na,na),sij(nloc,na,na)
     3570        real ment(nloc,na,na),sigij(nloc,na,na)
     3571        real clw(nloc,nd),elij(nloc,na,na)
     3572        real ep(nloc,na)
     3573        integer icb(nloc),inb(nloc)
     3574        real VPrecip(nloc,nd+1)
    35493575c ouputs:
    35503576        real da(nloc,na),phi(nloc,na,na)
     3577        real phi2(nloc,na,na)
     3578        real d1a(nloc,na),dam(nloc,na)
     3579! variables pour tracer dans precip de l'AA et des mel
    35513580c local variables:
    35523581        integer i,j,k
    3553 c       
    3554         da(:,:)=0.
    3555 c
     3582        real epm(nloc,na,na)
     3583c
     3584! variables d'Emanuel : du second indice au troisieme
     3585! --->    tab(i,k,j) -> de l origine k a l arrivee j
     3586!  ment, sigij, elij
     3587! variables personnelles : du troisieme au second indice
     3588! --->    tab(i,j,k) -> de k a j
     3589! phi, phi2
     3590!
     3591! initialisations
     3592c
     3593       da(:,:)=0.
     3594       d1a(:,:)=0.
     3595       dam(:,:)=0.
     3596       epm(:,:,:)=0.
     3597c
     3598!  fraction deau condensee dans les melanges convertie en precip
     3599        do j=1,na
     3600         do k=1,na
     3601           do i=1,ncum
     3602            if(k.ge.icb(i).and.k.le.inb(i).and.
     3603     &         j.ge.k.and.j.le.inb(i)) then
     3604             epm(i,j,k)=1.-(1.-ep(i,j))*clw(i,j)/elij(i,k,j)
     3605             epm(i,j,k)=max(epm(i,j,k),0.0)
     3606            endif
     3607           end do
     3608         end do
     3609        end do
     3610
     3611!  matrices pour calculer la tendance des concentrations dans cvltr.F90
    35563612        do j=1,na
    35573613          do k=1,na
    35583614            do i=1,ncum
    3559             da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
    3560             phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
     3615             da(i,j)=da(i,j)+(1.-sigij(i,k,j))*ment(i,k,j)
     3616             phi(i,j,k)=sigij(i,k,j)*ment(i,k,j)
     3617             d1a(i,j)=d1a(i,j)+ment(i,k,j)*ep(i,k)
     3618     &              *(1.-sigij(i,k,j))
     3619            if(k.le.j) then
     3620             dam(i,j)=dam(i,j)+ment(i,k,j)
     3621     &             *epm(i,k,j)*(1.-ep(i,k))*(1.-sigij(i,k,j))
     3622
     3623             phi2(i,j,k)=phi(i,j,k)*epm(i,j,k)   
     3624            else
     3625             dam(i,j)=0.
     3626             phi2(i,j,k)=0.
     3627            endif
    35613628            end do
    35623629          end do
    35633630        end do
     3631   
    35643632        return
    35653633        end
    3566 !AC!
     3634!AC! et !RomP <<<
    35673635
    35683636      SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum

  • LMDZ5/trunk/libf/phylmd/cv3a_uncompress.F

    r1650 r1742  
    99     :         ,Plim1,Plim2,asupmax,supmax0
    1010     :         ,asupmaxmin
    11 !AC!
    12      :         ,da,phi
    13 !AC!
     11!
     12     :         ,da,phi                               !AC!
     13     :         ,mp,phi2,d1a,dam,sigij                  !RomP
     14     :         ,wdtrainA,wdtrainM,elij,clw           !RomP
     15     :         ,evap,ep                              !RomP
     16!
    1417     o         ,iflag1,kbas1,ktop1
    1518     :         ,precip1,cbmf1,plcl1,plfc1,wbeff1,sig1,w01,ptop21
    … …  
    2124     :         ,Plim11,Plim21,asupmax1,supmax01
    2225     :         ,asupmaxmin1     
    23 !AC!
    24      :         ,da1,phi1  )
    25 !AC!
     26!
     27     o         ,da1,phi1                             !AC!
     28     o         ,mp1,phi21,d1a1,dam1,sigij1             !RomP
     29     o         ,wdtrainA1,wdtrainM1,elij1,clw1       !RomP
     30     o         ,evap1,ep1)                           !RomP
     31!
    2632***************************************************************
    2733*                                                             *
    … …  
    5662      real asupmax(nloc,nd),supmax0(nloc)
    5763      real asupmaxmin(nloc)
    58 !AC!
    59       real da(nloc,nd),phi(nloc,nd,nd)
    60 !AC!
     64!
     65      real da(nloc,nd),phi(nloc,nd,nd)                    !AC!
     66      real mp(nloc,nd)                                    !RomP
     67      real phi2(nloc,nd,nd)                               !RomP
     68      real d1a(nloc,nd),dam(nloc,nd)                      !RomP
     69      real wdtrainA(nloc,nd), wdtrainM(nloc,nd)           !RomP
     70      real sigij(nloc,nd,nd)                                !RomP
     71      real elij(nloc,nd,nd),clw(nloc,nd)                  !RomP
     72      real evap(nloc,nd),ep(nloc,nd)                      !RomP
     73!
    6174c outputs:
    6275      integer iflag1(len),kbas1(len),ktop1(len)
    … …  
    7689      real asupmax1(len,nd),supmax01(len)
    7790      real asupmaxmin1(len)
    78 !AC!
    79       real da1(nloc,nd),phi1(nloc,nd,nd)
    80 !AC!
     91!
     92      real da1(nloc,nd),phi1(nloc,nd,nd)                  !AC!
     93      real mp1(nloc,nd)                                   !RomP
     94      real phi21(nloc,nd,nd)                              !RomP
     95      real d1a1(nloc,nd),dam1(nloc,nd)                    !RomP
     96      real wdtrainA1(len,nd), wdtrainM1(len,nd)           !RomP
     97      real sigij1(len,nd,nd)                                !RomP
     98      real elij1(len,nd,nd),clw1(len,nd)                  !RomP
     99      real evap1(len,nd),ep1(len,nd)                      !RomP
     100!
    81101c
    82102c local variables:
    … …  
    122142            fqd1(idcum(i),k)=fqd(i,k)
    123143            asupmax1(idcum(i),k)=asupmax(i,k)
    124 !AC!
    125             da1(idcum(i),k)=da(i,k)
    126 !AC!
     144!
     145            da1(idcum(i),k)=da(i,k)                       !AC!
     146            mp1(idcum(i),k)      = mp(i,k)                !RomP
     147            d1a1(idcum(i),k)     = d1a(i,k)               !RomP
     148            dam1(idcum(i),k)     = dam(i,k)               !RomP
     149            wdtrainA1(idcum(i),k)= wdtrainA(i,k)          !RomP
     150            wdtrainM1(idcum(i),k)= wdtrainM(i,k)          !RomP
     151            clw1(idcum(i),k)     = clw(i,k)               !RomP
     152            evap1(idcum(i),k)    = evap(i,k)              !RomP
     153            ep1(idcum(i),k)      = ep(i,k)                !RomP
     154!
    127155 2010    continue
    128156 2020   continue
    … …  
    146174         do k=1,nd
    147175          do i=1,ncum
    148             phi1(idcum(i),k,j)=phi(i,k,j)
     176            phi1(idcum(i),k,j)=phi(i,k,j)                 !AC!
     177            phi21(idcum(i),k1,k2)= phi2(idcum(i),k1,k2)   !RomP
     178            sigij1(idcum(i),k1,k2) = sigij(idcum(i),k1,k2)    !RomP
     179            elij1(idcum(i),k1,k2)= elij(idcum(i),k1,k2)   !RomP
    149180          end do
    150181         end do
    … …  
    157188c          do 2200 i=1,ncum
    158189c            ment1(idcum(i),k1,k2) = ment(i,k1,k2)
    159 c            sij1(idcum(i),k1,k2) = sij(i,k1,k2)
     190c            sigij1(idcum(i),k1,k2) = sigij(i,k1,k2)
    160191c2200      enddo
    161192c2210     enddo

  • LMDZ5/trunk/libf/phylmd/cv3p_mixing.F

    r1650 r1742  
    33     :                    ,unk,vnk,hp,tv,tvp,ep,clw,sig
    44     :                    ,ment,qent,hent,uent,vent,nent
    5      :                   ,sij,elij,supmax,ments,qents,traent)
     5     :                   ,sigij,elij,supmax,ments,qents,traent)
    66***************************************************************
    77*                                                             *
    … …  
    3636      real ment(nloc,na,na), qent(nloc,na,na)
    3737      real uent(nloc,na,na), vent(nloc,na,na)
    38       real sij(nloc,na,na), elij(nloc,na,na)
     38      real sigij(nloc,na,na), elij(nloc,na,na)
    3939      real supmax(nloc,na)     ! Highest mixing fraction of mixed updraughts
    4040                               ! with the sign of (h-hp)
    4141      real traent(nloc,nd,nd,ntra)
    4242      real ments(nloc,nd,nd), qents(nloc,nd,nd)
    43       real sigij(nloc,nd,nd)
    4443      real hent(nloc,nd,nd)
    4544      integer nent(nloc,nd)
    … …  
    5756      real Sbef(nloc), Sup(nloc), Smin(nloc)
    5857      real asij(nloc), smax(nloc), scrit(nloc)
     58      real sij(nloc,nd,nd)
    5959      real csum(nloc,nd)
    6060      real awat

  • LMDZ5/trunk/libf/phylmd/cv_driver.F

    r965 r1742  
    99     &                   icb1,inb1,
    1010     &                   delt,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1,
    11      &                   da1,phi1,mp1)
     11     &                   da1,phi1,mp1,phi21,d1a1,dam1,sij1,clw1,elij1,
     12     &                   evap1,ep1,epmlmMm1,eplaMm1,
     13     &                   wdtrainA1,wdtrainM1)
    1214C
    1315      USE dimphy
    … …  
    6769cym#include "dimensions.h"
    6870cym#include "dimphy.h"
    69 
     71c
     72c Input
    7073      integer len
    7174      integer nd
    … …  
    7477      integer iflag_con
    7578      integer ntra
     79      real delt
    7680      real t1(len,nd)
    7781      real q1(len,nd)
    … …  
    7983      real u1(len,nd)
    8084      real v1(len,nd)
     85      real tra1(len,nd,ntra)
    8186      real p1(len,nd)
    8287      real ph1(len,ndp1)
     88c
     89c Output
    8390      integer iflag1(len)
    8491      real ft1(len,nd)
    … …  
    8693      real fu1(len,nd)
    8794      real fv1(len,nd)
     95      real ftra1(len,nd,ntra)
    8896      real precip1(len)
    8997      real cbmf1(len)
    90       real VPrecip1(len,nd+1)
     98      real sig1(klon,klev)
     99      real w01(klon,klev)
     100      real VPrecip1(len,nd+1)
     101      real evap1(len,nd)                    !RomP
     102      real ep1(len,nd)                      !RomP
    91103      real Ma1(len,nd)
    92104      real upwd1(len,nd)
    … …  
    98110      real cape1(len)     
    99111
     112! RomP >>>
     113      real wdtrainA1(len,nd), wdtrainM1(len,nd)
     114      real sij1(len,nd,nd),elij1(len,nd,nd)
    100115      real da1(len,nd),phi1(len,nd,nd),mp1(len,nd)
    101       real da(len,nd),phi(len,nd,nd),mp(len,nd)
    102       real tra1(len,nd,ntra)
    103       real ftra1(len,nd,ntra)
    104 
    105       real delt
     116
     117      real phi21(len,nd,nd)
     118      real d1a1(len,nd), dam1(len,nd)
     119      real epmlmMm1(len,nd,nd),eplaMm1(len,nd)
     120! RomP <<<
     121
    106122
    107123!-------------------------------------------------------------------
    … …  
    243259      real tvp1(klon,klev)
    244260      real clw1(klon,klev)
    245       real sig1(klon,klev)
    246       real w01(klon,klev)
    247261      real th1(klon,klev)
    248262c
    … …  
    277291      real ments(nloc,klev,klev), qents(nloc,klev,klev)
    278292      real sij(nloc,klev,klev), elij(nloc,klev,klev)
     293! RomP >>>
     294      real da(nloc,klev),phi(nloc,klev,klev),mp(nloc,klev)
     295      real epmlmMm(nloc,klev,klev),eplaMm(nloc,klev) 
     296      real phi2(nloc,klev,klev)
     297      real d1a(nloc,klev), dam(nloc,klev)
     298      real wdtrainA(nloc,klev),wdtrainM(nloc,klev)
     299      real sigd(nloc)
     300! RomP <<<
    279301      real qp(nloc,klev), up(nloc,klev), vp(nloc,klev)
    280302      real wt(nloc,klev), water(nloc,klev), evap(nloc,klev)
    … …  
    295317! --- SET CONSTANTS AND PARAMETERS
    296318!-------------------------------------------------------------------
    297 
     319      print *, '-> cv_driver'      !jyg
    298320c -- set simulation flags:
    299321c   (common cvflag)
    … …  
    325347!---------------------------------------------------------------------
    326348
    327       do 20 k=1,nd
    328         do 10 i=1,len
    329          ft1(i,k)=0.0
    330          fq1(i,k)=0.0
    331          fu1(i,k)=0.0
    332          fv1(i,k)=0.0
    333          tvp1(i,k)=0.0
    334          tp1(i,k)=0.0
    335          clw1(i,k)=0.0
     349         ft1(:,:)=0.0
     350         fq1(:,:)=0.0
     351         fu1(:,:)=0.0
     352         fv1(:,:)=0.0
     353         tvp1(:,:)=0.0
     354         tp1(:,:)=0.0
     355         clw1(:,:)=0.0
    336356cym
    337          clw(i,k)=0.0   
    338          gz1(i,k) = 0.
    339          VPrecip1(i,k) = 0.
    340          Ma1(i,k)=0.0
    341          upwd1(i,k)=0.0
    342          dnwd1(i,k)=0.0
    343          dnwd01(i,k)=0.0
    344          qcondc1(i,k)=0.0
    345  10     continue
    346  20   continue
    347 
    348       do 30 j=1,ntra
    349        do 31 k=1,nd
    350         do 32 i=1,len
    351          ftra1(i,k,j)=0.0
    352  32     continue   
    353  31    continue   
    354  30   continue   
    355 
    356       do 60 i=1,len
    357         precip1(i)=0.0
    358         iflag1(i)=0
    359         wd1(i)=0.0
    360         cape1(i)=0.0
    361         VPrecip1(i,nd+1)=0.0
    362  60   continue
     357         clw(:,:)=0.0   
     358         gz1(:,:) = 0.
     359         VPrecip1(:,:) = 0.
     360         Ma1(:,:)=0.0
     361         upwd1(:,:)=0.0
     362         dnwd1(:,:)=0.0
     363         dnwd01(:,:)=0.0
     364         qcondc1(:,:)=0.0
     365
     366         ftra1(:,:,:)=0.0
     367
     368        elij1(:,:,:) = 0.0
     369        sij1(:,:,:) = 0.0
     370
     371        precip1(:)=0.0
     372        iflag1(:)=0
     373        wd1(:)=0.0
     374        cape1(:)=0.0
    363375
    364376      if (iflag_con.eq.30) then
    … …  
    441453 400  continue
    442454
    443 c       print*,'klon, ncum = ',len,ncum
     455         print*,'cv_driver : klon, ncum = ',len,ncum
    444456
    445457      IF (ncum.gt.0) THEN
    … …  
    541553
    542554      if (iflag_con.eq.30) then
    543        CALL cv30_unsat(nloc,ncum,nd,nd,ntra,icb,inb    ! na->nd
     555!RomP >>>
     556       CALL cv30_unsat(nloc,ncum,nd,nd,ntra,icb,inb         ! na->nd
    544557     :               ,t,q,qs,gz,u,v,tra,p,ph
    545558     :               ,th,tv,lv,cpn,ep,sigp,clw
    546559     :               ,m,ment,elij,delt,plcl
    547      o          ,mp,qp,up,vp,trap,wt,water,evap,b)
     560     :          ,mp,qp,up,vp,trap,wt,water,evap,b
     561     o          ,wdtrainA,wdtrainM)
     562!RomP <<<
    548563      endif
    549564
    … …  
    588603
    589604      if (iflag_con.eq.30) then
     605!RomP >>>
    590606       CALL cv30_tracer(nloc,len,ncum,nd,nd,
    591      :                  ment,sij,da,phi)
     607     :                  ment,sij,da,phi,phi2,d1a,dam,
     608     :                  ep,VPrecip,elij,clw,epmlmMm,eplaMm,
     609     :                  icb,inb)
     610!RomP <<<
    592611      endif
    593612
    … …  
    603622       CALL cv30_uncompress(nloc,len,ncum,nd,ntra,idcum
    604623     :          ,iflag
    605      :          ,precip,VPrecip,sig,w0
     624     :          ,precip,VPrecip,evap,ep,sig,w0        !RomP
    606625     :          ,ft,fq,fu,fv,ftra
    607      :          ,inb 
     626     :          ,inb
    608627     :          ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
    609      :          ,da,phi,mp
     628     :          ,da,phi,mp,phi2,d1a,dam,sij           !RomP
     629     :          ,elij,clw,epmlmMm,eplaMm              !RomP 
     630     :          ,wdtrainA,wdtrainM                    !RomP
    610631     o          ,iflag1
    611      o          ,precip1,VPrecip1,sig1,w01
     632     o          ,precip1,VPrecip1,evap1,ep1,sig1,w01     !RomP
    612633     o          ,ft1,fq1,fu1,fv1,ftra1
    613634     o          ,inb1
    614635     o          ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
    615      o          ,da1,phi1,mp1)
     636     o          ,da1,phi1,mp1,phi21,d1a1,dam1,sij1       !RomP
     637     o          ,elij1,clw1,epmlmMm1,eplaMm1             !RomP
     638     o          ,wdtrainA1,wdtrainM1)                    !RomP
    616639      endif
    617640
    … …  
    6326559999  continue
    633656
     657      print *, 'fin cv_driver ->'      !jyg
    634658      return
    635659      end

  • LMDZ5/trunk/libf/phylmd/cva_driver.F

    r1652 r1742  
    2121     &                   Plim11,Plim21,asupmax1,supmax01,asupmaxmin1
    2222     &                   ,lalim_conv,
    23 !AC!
    24      &                   da1,phi1)
    25 !AC!
     23     &                   da1,phi1,mp1,phi21,d1a1,dam1,sigij1,clw1,     ! RomP
     24     &                   elij1,evap1,ep1,                              ! RomP
     25     &                   wdtrainA1,wdtrainM1)                          ! RomP
    2626***************************************************************
    2727*                                                             *
    … …  
    175175c
    176176!AC!
    177       real da1(len,nd),phi1(len,nd,nd)
    178       real da(len,nd),phi(len,nd,nd)
     177!!      real da1(len,nd),phi1(len,nd,nd)
     178!!      real da(len,nd),phi(len,nd,nd)
    179179!AC!
    180180      real ftd1(len,nd)
    … …  
    186186      real asupmaxmin1(len)
    187187      integer lalim_conv(len)
     188! RomP >>>
     189      real wdtrainA1(len,nd), wdtrainM1(len,nd)
     190      real wdtrainA(nloc,klev),wdtrainM(nloc,klev)
     191      real da1(len,nd),phi1(len,nd,nd),mp1(len,nd)
     192      real da(len,nd),phi(len,nd,nd)
     193      real evap1(len,nd),ep1(len,nd)
     194      real sigij1(len,nd,nd),elij1(len,nd,nd)
     195      real phi2(len,nd,nd)
     196      real d1a(len,nd), dam(len,nd)
     197      real phi21(len,nd,nd)
     198      real d1a1(len,nd), dam1(len,nd)
     199! RomP <<<
    188200!-------------------------------------------------------------------
    189201! --- ARGUMENTS
    … …  
    397409      real cin(nloc)
    398410      real m(nloc,klev)
    399       real ment(nloc,klev,klev), sij(nloc,klev,klev)
     411      real ment(nloc,klev,klev), sigij(nloc,klev,klev)
    400412      real qent(nloc,klev,klev)
    401413      real hent(nloc,klev,klev)
    … …  
    505517      call  zilch(cbmf1   ,nword1)
    506518      call  zilch(ptop21  ,nword1)
    507       sigd1=0.
     519      sigd1(:)=0.
    508520      call  zilch(Ma1     ,nword2)
    509521      call  zilch(mip1    ,nword2)
    … …  
    771783     :                       ,unk,vnk,hp,tv,tvp,ep,clw,sig
    772784     :                    ,ment,qent,hent,uent,vent,nent
    773      :                   ,sij,elij,supmax,ments,qents,traent)
     785     :                   ,sigij,elij,supmax,ments,qents,traent)
    774786!        print*, 'cv3p_mixing-> supmax ', (supmax(1,k), k=1,nd)
    775787     
    … …  
    828840     :                       ,ph,t,q,qs,u,v,tra,h,lv,qnk
    829841     :                       ,unk,vnk,hp,tv,tvp,ep,clw,m,sig
    830      o   ,ment,qent,uent,vent,nent,sij,elij,ments,qents,traent)
     842     o   ,ment,qent,uent,vent,nent,sigij,elij,ments,qents,traent)
    831843         CALL zilch(hent,nloc*klev*klev)
    832844        ELSE
    … …  
    842854     :                     ,ph,t,q,qs,u,v,h,lv,qnk
    843855     :                     ,hp,tv,tvp,ep,clw,cbmf
    844      o                     ,m,ment,qent,uent,vent,nent,sij,elij)
     856     o                     ,m,ment,qent,uent,vent,nent,sigij,elij)
    845857      endif
    846858c
    … …  
    865877     :               ,ep,sigp,clw
    866878     :               ,m,ment,elij,delt,plcl,coef_clos
    867      o          ,mp,qp,up,vp,trap,wt,water,evap,b,sigd)
     879     o          ,mp,qp,up,vp,trap,wt,water,evap,b,sigd
     880     o          ,wdtrainA,wdtrainM)   ! RomP
    868881      endif
    869882     
    … …  
    925938
    926939      if (iflag_con.eq.3) then
     940!RomP >>>
    927941       CALL cv3_tracer(nloc,len,ncum,nd,nd,
    928      :                  ment,sij,da,phi)
     942     :                  ment,sigij,da,phi,phi2,d1a,dam,
     943     :                  ep,Vprecip,elij,clw,icb,inb)
     944!RomP <<<
    929945      endif
    930946
    … …  
    947963     :          ,Plim1,Plim2,asupmax,supmax0
    948964     :          ,asupmaxmin
    949 !AC!
    950      :          ,da,phi
    951 !AC!
     965     :          ,da,phi,mp,phi2,d1a,dam,sigij         ! RomP
     966     :          ,wdtrainA,wdtrainM,elij,clw           ! RomP
     967     :          ,evap,ep                              ! RomP
    952968     o          ,iflag1,kbas1,ktop1
    953969     o          ,precip1,cbmf1,plcl1,plfc1,wbeff1,sig1,w01,ptop21
    … …  
    959975     o          ,Plim11,Plim21,asupmax1,supmax01
    960976     o          ,asupmaxmin1
    961 !AC!
    962      o          ,da1,phi1)
    963 !AC!
     977     o          ,da1,phi1,mp1,phi21,d1a1,dam1,sigij1  ! RomP
     978     o          ,wdtrainA1,wdtrainM1,elij1,clw1       ! RomP
     979     o          ,evap1,ep1)                           ! RomP
    964980      endif
    965981

  • LMDZ5/trunk/libf/phylmd/cvltr.F90

    r1279 r1742  
    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

  • LMDZ5/trunk/libf/phylmd/fisrtilp.F90

    r1575 r1742  
    44!
    55SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
    6      d_t, d_q, d_ql, rneb, radliq, rain, snow, &
    7      pfrac_impa, pfrac_nucl, pfrac_1nucl, &
    8      frac_impa, frac_nucl, &
    9      prfl, psfl, rhcl, zqta, fraca, &
     6     d_t, d_q, d_ql, rneb, radliq, rain, snow,          &
     7     pfrac_impa, pfrac_nucl, pfrac_1nucl,               &
     8     frac_impa, frac_nucl, beta,                        &
     9     prfl, psfl, rhcl, zqta, fraca,                     &
    1010     ztv, zpspsk, ztla, zthl, iflag_cldcon)
    1111
    … …  
    124124  REAL zprec_cond(klon)
    125125  !AA
     126! RomP >>> 15 nov 2012
     127  REAL   beta(klon,klev) ! taux de conversion de l'eau cond
     128! RomP <<<
    126129  REAL zmair, zcpair, zcpeau
    127130  !     Pour la conversion eau-neige
    … …  
    171174           pfrac_1nucl(i,k)=1.
    172175           pfrac_impa(i,k)=1.
     176           beta(i,k)=0.  !RomP initialisation
    173177        ENDDO
    174178     ENDDO
    … …  
    549553     DO i = 1,klon
    550554        !
     555        if(zcond(i).gt.zoliq(i)+1.e-10) then
     556         beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
     557        else
     558         beta(i,k) = 0.
     559        endif
    551560        zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
    552561             * (paprs(i,k)-paprs(i,k+1))/RG

  • LMDZ5/trunk/libf/phylmd/ini_histrac.h

    r1570 r1742  
    44  IF (ecrit_tra>0.) THEN
    55!$OMP MASTER
    6      CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian)
     6!!!     CALL ymds2ju(annee_ref, 1, day_ref, 0.0, zjulian)
     7! correction pour l heure initiale                               !jyg
     8!                                                               !jyg
     9      CALL ymds2ju(annee_ref, 1, day_ref, hour, zjulian)         !jyg
     10
    711     CALL histbeg_phy("histrac", itau_phy, zjulian, pdtphys,nhori, nid_tra)
    8      CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",klev, presnivs, nvert)
     12     CALL histvert(nid_tra, "presnivs", "Vertical levels", "Pa",klev, presnivs, nvert,"down")
    913
    1014     zsto = pdtphys
    … …  
    1721          "kg m-2", iim, jj_nb, nhori, klev, 1, klev, nvert, 32, "ave(X)", &
    1822          zsto,zout)
     23! RomP >>>
     24     CALL histdef(nid_tra, "sourceBE", "source 7Be", &
     25          "at/kgA/s", iim, jj_nb, nhori, klev, 1, klev, nvert, 32, "ave(X)", &
     26          zsto,zout)
     27! RomP <<<
    1928
    2029!TRACEURS
    … …  
    3039        IF (lessivage .AND. aerosol(it)) THEN
    3140           CALL histdef(nid_tra, "fl"//tname(iiq),"Flux "//ttext(iiq), &
    32                 "U/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
    33                 "ave(X)", zsto,zout)
    34         END IF
     41                "at/m2/s",iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     42                "ave(X)", zsto,zout)
     43           CALL histdef(nid_tra, "d_tr_ls_"//tname(iiq),      &
     44                "tendance lessivage large scale"// ttext(iiq), "?",&
     45                iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     46                "ave(X)", zsto,zout)
     47           CALL histdef(nid_tra, "d_tr_insc_"//tname(iiq),      &
     48                "tendance lessivage large scale"// ttext(iiq), "?",&
     49                iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     50                "ave(X)", zsto,zout)
     51           CALL histdef(nid_tra, "d_tr_bcscav_"//tname(iiq),      &
     52                "tendance lessivage large scale"// ttext(iiq), "?",&
     53                iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     54                "ave(X)", zsto,zout)
     55           CALL histdef(nid_tra, "d_tr_evls_"//tname(iiq),      &
     56                "tendance lessivage large scale"// ttext(iiq), "?",&
     57                iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
     58                "ave(X)", zsto,zout)
     59!  Tracer concentration in LS precipitation at surface
     60           CALL histdef(nid_tra, "qpr_ls_"//tname(iiq),       &
     61                "concentration in LS precip"// ttext(iiq), "at/kgw", &
     62                iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     63                "ave(X)", zsto,zout)     
     64                 END IF
    3565
    3666! TD THERMIQUES
    … …  
    5080        ENDIF
    5181
     82! RomP >>>
     83        IF (iflag_con.EQ.30) THEN
     84           CALL histdef(nid_tra, "d_tr_cvMA_"//tname(iiq),   &
     85                "tendance convection"// ttext(iiq), "?",&
     86                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     87                "ave(X)", zsto,zout)
     88           CALL histdef(nid_tra, "d_tr_trsp_"//tname(iiq),   &
     89                "tendance transport "// ttext(iiq), "at/kga",   &
     90                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     91                "ave(X)", zsto,zout)
     92           CALL histdef(nid_tra, "d_tr_sscav_"//tname(iiq),   &
     93                "tendance lessivage flux satures "// ttext(iiq), "at/kga",   &
     94                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     95                "ave(X)", zsto,zout)
     96           CALL histdef(nid_tra, "d_tr_sat_"//tname(iiq),   &
     97                "tendance flux satures "// ttext(iiq), "at/kga",  &
     98                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     99                "ave(X)", zsto,zout)
     100           CALL histdef(nid_tra, "d_tr_uscav_"//tname(iiq),  &
     101                "tendance flux insatures "// ttext(iiq), "at/kga", &
     102                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     103                "ave(X)", zsto,zout)
     104           CALL histdef(nid_tra, "tr_pr_"//tname(iiq),  &
     105                "concentration dans precip"// ttext(iiq), "at/kga", &
     106                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     107                "ave(X)", zsto,zout)
     108           CALL histdef(nid_tra, "tr_aa_"//tname(iiq),  &
     109                "concentration precip issu AA"// ttext(iiq), "at/kga", &
     110                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     111                "ave(X)", zsto,zout)
     112           CALL histdef(nid_tra, "tr_mel_"//tname(iiq),  &
     113                "concentration precip issu melange"// ttext(iiq), "at/kga", &
     114                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     115                "ave(X)", zsto,zout)
     116           CALL histdef(nid_tra, "tr_di_"//tname(iiq),  &
     117                "concentration dans descente insaturee"// ttext(iiq), "at/kga", &
     118                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     119                "ave(X)", zsto,zout)
     120           CALL histdef(nid_tra, "tr_trspdi_"//tname(iiq),  &
     121                "conc descente insaturee MA"// ttext(iiq), "at/kga", &
     122                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     123                "ave(X)", zsto,zout)
     124           CALL histdef(nid_tra, "zmfd1a_"//tname(iiq),  &
     125                "zmfd1a"// ttext(iiq), "_", &
     126                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     127                "ave(X)", zsto,zout)
     128           CALL histdef(nid_tra, "zmfphi2_"//tname(iiq),  &
     129                "zmfphi2"// ttext(iiq), "_", &
     130                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     131                "ave(X)", zsto,zout)
     132           CALL histdef(nid_tra, "zmfdam_"//tname(iiq),  &
     133                "zmfdam"// ttext(iiq), "_", &
     134                iim,jj_nb,nhori, klev,1,klev,nvert, 32,    &
     135                "ave(X)", zsto,zout)
     136          ENDIF
     137! RomP <<<
     138           CALL histdef(nid_tra, "dtrdyn_"//tname(iiq),  &
     139                "td dyn tra"// ttext(iiq), "at/kga", &
     140                iim,jj_nb,nhori, klev,1,klev,nvert, 32,  &
     141                "ave(X)", zsto,zout)
     142! TD decroissance radioactive
     143           CALL histdef(nid_tra, "d_tr_dec_"//tname(iiq),   &
     144                "tendance decroi radio "// ttext(iiq), "",  &
     145                iim,jj_nb,nhori, klev,1,klev,nvert, 32,  &
     146                "ave(X)", zsto,zout)
     147
    52148! TD COUCHE-LIMITE
     149      IF (couchelimite) THEN
    53150        CALL histdef(nid_tra, "d_tr_cl_"//tname(iiq),      &
    54151             "tendance couche limite"// ttext(iiq), "?",   &
    55152             iim,jj_nb,nhori, klev,1,klev,nvert, 32,       &
    56153             "ave(X)", zsto,zout)
     154!  Dry deposit (1st layer and surface)
     155        CALL histdef(nid_tra, "d_tr_dry_"//tname(iiq),       &
     156             "tendancy dry deposit"// ttext(iiq), "at/kga/step", &
     157             iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     158             "ave(X)", zsto,zout)     
     159        CALL histdef(nid_tra, "flux_tr_dry_"//tname(iiq),       &
     160             "dry deposit at surf (downward)"// ttext(iiq), "at/m2/step", &
     161             iim,jj_nb,nhori, 1,1,1, -99, 32,                 &
     162             "ave(X)", zsto,zout)     
     163      ENDIF
    57164     ENDDO
     165
     166     CALL histdef(nid_tra, "Mint", "Mint","",         &
     167          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     168          "inst(X)", zout,zout)
     169     CALL histdef(nid_tra, "frac_impa", "frac_impa","",         &
     170          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     171          "inst(X)", zout,zout)
     172     CALL histdef(nid_tra, "frac_nucl", "frac_nucl","",         &
     173          iim,jj_nb,nhori, klev,1,klev,nvert, 32,          &
     174          "inst(X)", zout,zout)
    58175!---------------   
    59176!

  • LMDZ5/trunk/libf/phylmd/init_be.F90

    r1279 r1742  
    11!$Id $
    22
    3 SUBROUTINE init_be(pctsrf,masktr,tautr,vdeptr,scavtr,srcbe)
     3SUBROUTINE init_be(pctsrf,pplay,masktr,tautr,vdeptr,scavtr,srcbe)
     4!!!SUBROUTINE init_be(pctsrf,masktr,tautr,vdeptr,scavtr,srcbe)
    45
    56  USE dimphy
    … …  
    2627!
    2728  REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf !Pourcentage de sol (f(nature du sol))
     29  REAL,DIMENSION(klon,klev), INTENT(IN) :: pplay  ! Pressions en milieu de couches
    2830!
    2931! Output Arguments
    … …  
    3739! Local Variables
    3840!
     41!!!  INTEGER              :: iref      ! numero d'un point oceanique donnant la grille de pression de reference
    3942  REAL,DIMENSION(klon) :: rlatgeo   ! latitudes geomagnetiques de la grille
    4043  REAL                 :: glt       ! latitude du pole geomagnetique
    4144  REAL                 :: glg       ! longitude du pole geomagnetique
    4245  REAL                 :: latgeo,qcos
    43   INTEGER              :: k,i
     46  INTEGER              :: k,i, kref, k2
     47  INTEGER              :: nref
     48  PARAMETER (nref=39)
     49  REAL,DIMENSION(nref) :: pref      ! grille de pression de reference (bas des couches)
     50  DATA pref  /   &
     51      101249.99999999994, 100387.17261011522, 99447.35334189111,  98357.43412194174,   &
     52      97046.47707771382,  95447.1116450629,   93496.85259615642,  91139.46548240296,   &
     53      88326.55568744117,  85019.60710580258,  81192.7404556645,   76836.48366938648,   &
     54      71962.81275769137,  66611.56331321516,  60857.914829743604, 54819.84484441629,   &
     55      48663.06257114699,  42598.95465845692,  36869.104365898806, 31709.927925633147,  &
     56      27296.757208636915, 23682.282929080895, 20766.025578936627, 18336.105961406534,  &
     57      16178.04816768436,  14168.286905562818, 12275.719926478887, 10507.798835225762,  &
     58      8876.585404909414,  7391.283929569539,  6057.514475749798,  4877.165909157005,   &
     59      3848.34936408203,   2965.444753540027,  2219.2391544640013, 1597.15366044666,    &
     60      1083.5531161631498, 660.1311067852655,  306.36072267002805 /
     61!$OMP THREADPRIVATE(pref)
    4462
    4563  WRITE(*,*)'PASSAGE init_be ...'
    4664
    47 ! Source actuellement definie pour klev = 19 et klev >= 39
    48   IF (klev /= 19 .AND. klev<39) CALL abort_gcm("init_be","Source du be7 necessite klev=19 ou klev>=39",1)
    49 !
     65! la source est maintenant définie independemment de la valeur de klev.
     66!!! Source actuellement definie pour klev = 19 et klev >= 39
     67!!  IF (klev /= 19 .AND. klev<39) CALL abort_gcm("init_be","Source du be7 necessite klev=19 ou klev>=39",1)
     68!!!
    5069! Definition des constantes
    5170! -------------------------
    … …  
    5372  vdeptr = 1.E-3
    5473  scavtr = 0.5
     74!!!!!jyg le 13/03/2013; puis 20/03/2013 : pref est maintenant une table.
     75!!!
     76!!!   Recherche d'un point rlat=0., rlon=180.
     77!!      iref=(klon+1)/2
     78!!      DO i = 1,klon
     79!!        IF (abs(rlatd(i)) .LT. 0.15 .AND. cos(rlond(i)) .LT. -0.85) iref=i
     80!!      ENDDO
     81!!!
     82!!!   Grille de pression de reference (= approx de sommets de couches)
     83!!      pref(1) = pplay(iref,1)+0.5*(pplay(iref,1)-pplay(iref,2))
     84!!      DO k = 2,klev
     85!!        pref(k) = 0.5*(pplay(iref,k-1)+pplay(iref,k))
     86!!      ENDDO
     87!!!
    5588
    5689  WRITE(*,*) '-------------- SOURCE DE BERYLLIUM ------------------- '
    … …  
    77110! 3-mettre la source de Be ds la bonne unite (en at/kgA/s)
    78111!
    79   glt=78.5*rpi/180.
    80   glg=69.0*rpi/180.
     112  glt =  78.5*rpi/180.
     113  glg = -69.0*rpi/180.
    81114
    82115  DO i = 1,klon
    83116     qcos=sin(glt)*sin(rlatd(i))
    84      qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)+glg)
     117!!jyg
     118!!     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)+glg)
     119     qcos=qcos+cos(glt)*cos(rlatd(i))*cos(rlond(i)-glg)
     120!!jyg end
    85121     IF ( qcos .LT. -1.) qcos = -1.
    86122     IF ( qcos .GT. 1.)  qcos = 1.
    … …  
    88124  ENDDO
    89125
    90 !===========================
    91 !  Cas 19 niveaux verticaux
    92 !===========================
    93   IF (klev.eq.19) then
     126!!!===========================
     127!!!  Cas 19 niveaux verticaux
     128!!!===========================
     129!!  IF (klev.eq.19) then
     130!!     DO k = 1,klev
     131!!        DO i = 1,klon
     132!!!!!jyg le 13/03/2013
     133!!!
     134!!! k est le niveau dans la grille locale
     135!!! Determination du niveau kref dans la grille de refernce
     136!!      kref = 1
     137!!      DO k2 = 1,klev
     138!!        IF (pref(k2) .GT. pplay(i,k)) kref=k2
     139!!      ENDDO
     140!!!!!
     141!!           latgeo=(180./rpi)*abs(rlatgeo(i))
     142!!           IF ( kref .EQ. 1 ) THEN
     143!!              IF (latgeo.GE.50.0) srcbe(i,k)=0.1
     144!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09
     145!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.07
     146!!           END IF
     147!!           IF ( kref .EQ. 2 ) THEN
     148!!              IF (latgeo.GE.50.0) srcbe(i,k)=0.12
     149!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.1
     150!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.09
     151!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07
     152!!           END IF
     153!!           IF ( kref .EQ. 3 ) THEN
     154!!              IF (latgeo.GE.50.0) srcbe(i,k)=0.14
     155!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12
     156!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.1
     157!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.09
     158!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07
     159!!           END IF
     160!!           IF ( kref .EQ. 4 ) THEN
     161!!              IF (latgeo.GE.50.0) srcbe(i,k)=0.175
     162!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
     163!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14
     164!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12
     165!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1
     166!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09
     167!!           END IF
     168!!           IF ( kref .EQ. 5 ) THEN
     169!!              IF (latgeo.GE.50.0) srcbe(i,k)=0.28
     170!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.26
     171!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.23
     172!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.175
     173!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.14
     174!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12
     175!!           END IF
     176!!           IF ( kref .EQ. 6 ) THEN
     177!!              IF (latgeo.GE.50.0) srcbe(i,k)=0.56
     178!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.49
     179!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.42
     180!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.28
     181!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.26
     182!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.245
     183!!           END IF
     184!!           IF ( kref .EQ. 7 ) THEN
     185!!              IF (latgeo.GE.50.0) srcbe(i,k)=1.05
     186!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.875
     187!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.7
     188!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.52
     189!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.44
     190!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.385
     191!!           END IF
     192!!           IF ( kref .EQ. 8 ) THEN
     193!!              IF (latgeo.GE.50.0) srcbe(i,k)=2.
     194!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.8
     195!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.5
     196!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=1.
     197!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.8
     198!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.75
     199!!           END IF
     200!!           IF ( kref .EQ. 9 ) THEN
     201!!              IF (latgeo.GE.50.0) srcbe(i,k)=4.
     202!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3.5
     203!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=3.
     204!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=2.5
     205!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.8
     206!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4
     207!!           END IF
     208!!           IF ( kref .EQ. 10 ) THEN
     209!!              IF (latgeo.GE.50.0) srcbe(i,k)=8.5
     210!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8.
     211!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=7.
     212!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=4.5
     213!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.5
     214!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.
     215!!           END IF
     216!!           IF ( kref .EQ. 11 ) THEN
     217!!              IF (latgeo.GE.50.0) srcbe(i,k)=17.
     218!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=15.
     219!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=11.
     220!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=8.
     221!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=5.
     222!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.
     223!!           END IF
     224!!           IF ( kref .EQ. 12 ) THEN
     225!!              IF (latgeo.GE.50.0) srcbe(i,k)=25.
     226!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22.
     227!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.
     228!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=11.
     229!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.5
     230!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
     231!!           END IF
     232!!           IF ( kref .EQ. 13 ) THEN
     233!!              IF (latgeo.GE.60.0) srcbe(i,k)=33.
     234!!              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=32.
     235!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=30.
     236!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=22.
     237!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
     238!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.
     239!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.
     240!!           END IF
     241!!           IF ( kref .EQ. 14 ) THEN
     242!!              IF (latgeo.GE.60.0) srcbe(i,k)=48.
     243!!              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
     244!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=36.
     245!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=26.
     246!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=17.5
     247!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
     248!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
     249!!           END IF
     250!!           IF ( kref .EQ. 15 ) THEN
     251!!              IF (latgeo.GE.70.0) srcbe(i,k)=58.
     252!!              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=57.
     253!!              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
     254!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=38.
     255!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25.
     256!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
     257!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
     258!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
     259!!           END IF
     260!!           IF ( kref .EQ. 16 ) THEN
     261!!              IF (latgeo.GE.70.0) srcbe(i,k)=70.
     262!!              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=65.
     263!!              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
     264!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32.
     265!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=20.
     266!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=13.
     267!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=9.
     268!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.5
     269!!           END IF
     270!!           IF ( kref .GE. 17 ) THEN
     271!!              IF (latgeo.GE.70.0) srcbe(i,k)=80.
     272!!              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70.
     273!!              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
     274!!              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=27.
     275!!              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.5
     276!!              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=12.
     277!!              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8.
     278!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
     279!!           END IF
     280!!        END DO
     281!!     END DO
     282!!  END IF ! fin de 19 niveaux verticaux
     283!!
     284!!!!!!  IF (klev .ge. 39) then
    94285     DO k = 1,klev
    95286        DO i = 1,klon
     287!!!jyg le 13/03/2013
     288!
     289! k est le niveau dans la grille locale
     290! Determination du niveau kref dans la grille de refernce
     291      kref = 1
     292      DO k2 = 1,nref
     293        IF (pref(k2) .GT. pplay(i,k)) kref=k2
     294      ENDDO
     295!!!
    96296           latgeo=(180./rpi)*abs(rlatgeo(i))
    97            IF ( k .EQ. 1 ) THEN
     297           IF ( kref .LE. 4 ) THEN
     298              IF (latgeo.GE.50.0) srcbe(i,k)=0.07
     299           END IF
     300           IF ( kref .EQ. 5 ) THEN
    98301              IF (latgeo.GE.50.0) srcbe(i,k)=0.1
    99               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09
    100               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.07
    101            END IF
    102            IF ( k .EQ. 2 ) THEN
    103               IF (latgeo.GE.50.0) srcbe(i,k)=0.12
    104               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.1
    105               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.09
     302              IF (latgeo.GE.20.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09
    106303              IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07
    107304           END IF
    108            IF ( k .EQ. 3 ) THEN
     305           IF ( kref .EQ. 6 ) THEN
    109306              IF (latgeo.GE.50.0) srcbe(i,k)=0.14
    110307              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12
    … …  
    113310              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07
    114311           END IF
    115            IF ( k .EQ. 4 ) THEN
    116               IF (latgeo.GE.50.0) srcbe(i,k)=0.175
    117               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
    118               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.14
    119               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.12
    120               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.1
    121               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09
    122            END IF
    123            IF ( k .EQ. 5 ) THEN
    124               IF (latgeo.GE.50.0) srcbe(i,k)=0.28
    125               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.26
    126               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.23
    127               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.175
    128               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.14
    129               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12
    130            END IF
    131            IF ( k .EQ. 6 ) THEN
    132               IF (latgeo.GE.50.0) srcbe(i,k)=0.56
    133               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.49
    134               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.42
    135               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.28
    136               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.26
    137               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.245
    138            END IF
    139            IF ( k .EQ. 7 ) THEN
    140               IF (latgeo.GE.50.0) srcbe(i,k)=1.05
    141               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.875
    142               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.7
    143               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.52
    144               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.44
    145               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.385
    146            END IF
    147            IF ( k .EQ. 8 ) THEN
    148               IF (latgeo.GE.50.0) srcbe(i,k)=2.
    149               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.8
    150               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=1.5
    151               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=1.
    152               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.8
    153               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.75
    154            END IF
    155            IF ( k .EQ. 9 ) THEN
    156               IF (latgeo.GE.50.0) srcbe(i,k)=4.
    157               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3.5
    158               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=3.
    159               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=2.5
    160               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=1.8
    161               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4
    162            END IF
    163            IF ( k .EQ. 10 ) THEN
    164               IF (latgeo.GE.50.0) srcbe(i,k)=8.5
    165               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8.
    166               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=7.
    167               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=4.5
    168               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=3.5
    169               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.
    170            END IF
    171            IF ( k .EQ. 11 ) THEN
    172               IF (latgeo.GE.50.0) srcbe(i,k)=17.
    173               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=15.
    174               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=11.
    175               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=8.
    176               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=5.
    177               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.
    178            END IF
    179            IF ( k .EQ. 12 ) THEN
    180               IF (latgeo.GE.50.0) srcbe(i,k)=25.
    181               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22.
    182               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.
    183               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=11.
    184               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=7.5
    185               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
    186            END IF
    187            IF ( k .EQ. 13 ) THEN
    188               IF (latgeo.GE.60.0) srcbe(i,k)=33.
    189               IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=32.
    190               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=30.
    191               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=22.
    192               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
    193               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.
    194               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.
    195            END IF
    196            IF ( k .EQ. 14 ) THEN
    197               IF (latgeo.GE.60.0) srcbe(i,k)=48.
    198               IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
    199               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=36.
    200               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=26.
    201               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=17.5
    202               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
    203               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
    204            END IF
    205            IF ( k .EQ. 15 ) THEN
    206               IF (latgeo.GE.70.0) srcbe(i,k)=58.
    207               IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=57.
    208               IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
    209               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=38.
    210               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=25.
    211               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=15.
    212               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=12.5
    213               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
    214            END IF
    215            IF ( k .EQ. 16 ) THEN
    216               IF (latgeo.GE.70.0) srcbe(i,k)=70.
    217               IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=65.
    218               IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=50.
    219               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=32.
    220               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=20.
    221               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=13.
    222               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=9.
    223               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.5
    224            END IF
    225            IF ( k .GE. 17 ) THEN
    226               IF (latgeo.GE.70.0) srcbe(i,k)=80.
    227               IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70.
    228               IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=45.
    229               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=27.
    230               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=17.5
    231               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=12.
    232               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=8.
    233               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
    234            END IF
    235         END DO
    236      END DO
    237   END IF ! fin de 19 niveaux verticaux
    238 
    239 !================================
    240 !  Cas 39 niveaux verticaux
    241 !================================
    242   IF (klev .ge. 39) then
    243      DO k = 1,klev
    244         DO i = 1,klon
    245            latgeo=(180./rpi)*abs(rlatgeo(i))
    246            IF ( k .LE. 4 ) THEN
    247               IF (latgeo.GE.50.0) srcbe(i,k)=0.07
    248            END IF
    249            IF ( k .EQ. 5 ) THEN
    250               IF (latgeo.GE.50.0) srcbe(i,k)=0.1
    251               IF (latgeo.GE.20.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.09
    252               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.07
    253            END IF
    254            IF ( k .EQ. 6 ) THEN
    255               IF (latgeo.GE.50.0) srcbe(i,k)=0.14
    256               IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.12
    257               IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.1
    258               IF (latgeo.GE.20.0 .and. latgeo.LT.30.0) srcbe(i,k)=0.09
    259               IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=0.07
    260            END IF
    261            IF ( k .EQ. 7 ) THEN
     312           IF ( kref .EQ. 7 ) THEN
    262313              IF (latgeo.GE.50.0) srcbe(i,k)=0.16
    263314              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
    … …  
    267318              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.09
    268319           END IF
    269            IF ( k .EQ. 8 ) THEN
     320           IF ( kref .EQ. 8 ) THEN
    270321              IF (latgeo.GE.50.0) srcbe(i,k)=0.175
    271322              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.16
    … …  
    275326              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.1
    276327           END IF
    277            IF ( k .EQ. 9 ) THEN
     328           IF ( kref .EQ. 9 ) THEN
    278329              IF (latgeo.GE.50.0) srcbe(i,k)=0.245
    279330              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.21
    … …  
    283334              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.12
    284335           END IF
    285            IF ( k .EQ. 10 ) THEN
     336           IF ( kref .EQ. 10 ) THEN
    286337              IF (latgeo.GE.50.0) srcbe(i,k)=0.31
    287338              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.28
    … …  
    291342              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.14
    292343           END IF
    293            IF ( k .EQ. 11 ) THEN
     344           IF ( kref .EQ. 11 ) THEN
    294345              IF (latgeo.GE.50.0) srcbe(i,k)=0.35
    295346              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.3
    … …  
    299350              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.16
    300351           END IF
    301            IF ( k .EQ. 12 ) THEN
     352           IF ( kref .EQ. 12 ) THEN
    302353              IF (latgeo.GE.40.0) srcbe(i,k)=0.5
    303354              IF (latgeo.GE.30.0 .and. latgeo.LT.40.0) srcbe(i,k)=0.4
    … …  
    306357              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.25
    307358           END IF
    308            IF ( k .EQ. 13 ) THEN
     359           IF ( kref .EQ. 13 ) THEN
    309360              IF (latgeo.GE.50.0) srcbe(i,k)=0.8
    310361              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=0.7
    … …  
    314365              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.35
    315366           END IF
    316            IF ( k .EQ. 14 ) THEN
     367           IF ( kref .EQ. 14 ) THEN
    317368              IF (latgeo.GE.50.0) srcbe(i,k)=1.2
    318369              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=1.
    … …  
    322373              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.4
    323374           END IF
    324            IF ( k .EQ. 15 ) THEN
     375           IF ( kref .EQ. 15 ) THEN
    325376              IF (latgeo.GE.60.0) srcbe(i,k)=1.75
    326377              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=1.8
    … …  
    331382              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.65
    332383           END IF
    333            IF ( k .EQ. 16 ) THEN
     384           IF ( kref .EQ. 16 ) THEN
    334385              IF (latgeo.GE.50.0) srcbe(i,k)=3.
    335386              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=2.5
    … …  
    339390              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=0.9
    340391           END IF
    341            IF ( k .EQ. 17 ) THEN
     392           IF ( kref .EQ. 17 ) THEN
    342393              IF (latgeo.GE.50.0) srcbe(i,k)=4.
    343394              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=3.
    … …  
    347398              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=1.4
    348399           END IF
    349            IF ( k .EQ. 18 ) THEN
     400           IF ( kref .EQ. 18 ) THEN
    350401              IF (latgeo.GE.50.0) srcbe(i,k)=7.
    351402              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=6.
    … …  
    355406              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=2.
    356407           END IF
    357            IF ( k .EQ. 19 ) THEN
     408           IF ( kref .EQ. 19 ) THEN
    358409              IF (latgeo.GE.50.0) srcbe(i,k)=8.5
    359410              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=8.
    … …  
    363414              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.
    364415           END IF
    365            IF ( k .EQ. 20 ) THEN
     416           IF ( kref .EQ. 20 ) THEN
    366417              IF (latgeo.GE.50.0) srcbe(i,k)=12.5
    367418              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=12.
    … …  
    371422              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=3.5
    372423           END IF
    373            IF ( k .EQ. 21 ) THEN
     424           IF ( kref .EQ. 21 ) THEN
    374425              IF (latgeo.GE.50.0) srcbe(i,k)=16.
    375426              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=13.
    … …  
    379430              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.
    380431           END IF
    381            IF ( k .EQ. 22 ) THEN
     432           IF ( kref .EQ. 22 ) THEN
    382433              IF (latgeo.GE.50.0) srcbe(i,k)=20.
    383434              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=17.5
    … …  
    387438              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=4.5
    388439           END IF
    389            IF ( k .EQ. 23 ) THEN
     440           IF ( kref .EQ. 23 ) THEN
    390441              IF (latgeo.GE.50.0) srcbe(i,k)=25.
    391442              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=22.
    … …  
    395446              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
    396447           END IF
    397            IF ( k .EQ. 24 ) THEN
     448           IF ( kref .EQ. 24 ) THEN
    398449              IF (latgeo.GE.50.0) srcbe(i,k)=28.
    399450              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=26.
    … …  
    403454              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.
    404455           END IF
    405            IF ( k .EQ. 25 ) THEN
     456           IF ( kref .EQ. 25 ) THEN
    406457              IF (latgeo.GE.50.0) srcbe(i,k)=33.
    407458              IF (latgeo.GE.40.0 .and. latgeo.LT.50.0) srcbe(i,k)=28.
    … …  
    411462              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=8.5
    412463           END IF
    413            IF ( k .EQ. 26 ) THEN
     464           IF ( kref .EQ. 26 ) THEN
    414465              IF (latgeo.GE.60.0) srcbe(i,k)=38.
    415466              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=36.
    … …  
    419470              IF (latgeo.GE.10.0 .and. latgeo.LT.20.0) srcbe(i,k)=11.5
    420471              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
    421               IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
    422            END IF
    423            IF ( k .EQ. 27 ) THEN
     472!!jyg
     473!!              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=6.
     474!!jyg end
     475           END IF
     476           IF ( kref .EQ. 27 ) THEN
    424477              IF (latgeo.GE.60.0) srcbe(i,k)=46.
    425478              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=44.
    … …  
    430483              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
    431484           END IF
    432            IF ( k .EQ. 28 ) THEN
     485           IF ( kref .EQ. 28 ) THEN
    433486              IF (latgeo.GE.60.0) srcbe(i,k)=53.
    434487              IF (latgeo.GE.50.0 .and. latgeo.LT.60.0) srcbe(i,k)=48.
    … …  
    439492              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
    440493           END IF
    441            IF ( k .EQ. 29 ) THEN
     494           IF ( kref .EQ. 29 ) THEN
    442495              IF (latgeo.GE.70.0) srcbe(i,k)=58.
    443496              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=56.
    … …  
    449502              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=10.
    450503           END IF
    451            IF ( k .EQ. 30 ) THEN
     504           IF ( kref .EQ. 30 ) THEN
    452505              IF (latgeo.GE.70.0) srcbe(i,k)=65.
    453506              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=60.
    … …  
    459512              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=9.
    460513           END IF
    461            IF ( k .EQ. 31 ) THEN
     514           IF ( kref .EQ. 31 ) THEN
    462515              IF (latgeo.GE.70.0) srcbe(i,k)=70.
    463516              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=62.
    … …  
    469522              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.6
    470523           END IF
    471            IF ( k .EQ. 32 ) THEN
     524           IF ( kref .EQ. 32 ) THEN
    472525              IF (latgeo.GE.70.0) srcbe(i,k)=80.
    473526              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=60.
    … …  
    479532              IF (latgeo.GE.0.0 .and. latgeo.LT.10.0) srcbe(i,k)=7.4
    480533           END IF
    481            IF ( k .GE. 33 ) THEN
     534           IF ( kref .GE. 33 ) THEN
    482535              IF (latgeo.GE.70.0) srcbe(i,k)=80.
    483536              IF (latgeo.GE.60.0 .and. latgeo.LT.70.0) srcbe(i,k)=70.
    … …  
    491544        END DO
    492545     END DO
    493   END IF ! fin de 39 niveaux verticaux
     546!!!!!!  END IF ! fin de 39 niveaux verticaux
    494547
    495548
    … …  
    503556       ! 1/(60*1.295) = 0.01287
    504557       srcbe(i,k)=srcbe(i,k)*0.01287
     558!!       print *,'  k, srcbe(i,k) ',   &
     559!!                  k, srcbe(i,k)
    505560       ! La source est  at/min/m3 -> at/s/m3
    506561       ! srcbe(i,k)=srcbe(i,k)*0.0166667

  • LMDZ5/trunk/libf/phylmd/orografi_strato.F

    r1492 r1742  
    20042004
    20052005      DO 110 JK=1,NLEV
    2006       ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2+1,1)
     2006      ZPM1R=pplay_glo(klon_glo/2+1,jk)/paprs_glo(klon_glo/2+1,1)
    20072007      IF(ZPM1R.GE.ZSIGT)THEN
    20082008         nktopg=JK
    20092009      ENDIF
    2010       ZPM1R=pplay_glo(klon_glo/2,jk)/paprs_glo(klon_glo/2+1,1)
     2010      ZPM1R=pplay_glo(klon_glo/2+1,jk)/paprs_glo(klon_glo/2+1,1)
    20112011      IF(ZPM1R.GE.ZTOP)THEN
    20122012         nstra=JK

  • LMDZ5/trunk/libf/phylmd/phys_local_var_mod.F90

    r1539 r1742  
    2424      REAL, SAVE, ALLOCATABLE :: d_u_dyn(:,:), d_v_dyn(:,:)
    2525      !$OMP THREADPRIVATE(d_u_dyn, d_v_dyn)
     26!!!!
     27      REAL, SAVE, ALLOCATABLE :: d_tr_dyn(:,:,:)
     28      !$OMP THREADPRIVATE(d_tr_dyn)
     29!!!!
    2630      REAL, SAVE, ALLOCATABLE :: d_t_con(:,:),d_q_con(:,:)
    2731      !$OMP THREADPRIVATE(d_t_con,d_q_con)
    … …  
    199203      allocate(d_t_dyn(klon,klev),d_q_dyn(klon,klev))
    200204      allocate(d_u_dyn(klon,klev),d_v_dyn(klon,klev))
     205      allocate(d_tr_dyn(klon,klev,nbtr))                   !RomP
    201206      allocate(d_t_con(klon,klev),d_q_con(klon,klev))
    202207      allocate(d_u_con(klon,klev),d_v_con(klon,klev))
    … …  
    287292      deallocate(d_t_dyn,d_q_dyn)
    288293      deallocate(d_u_dyn,d_v_dyn)
     294      deallocate(d_tr_dyn)                      !RomP
    289295      deallocate(d_t_con,d_q_con)
    290296      deallocate(d_u_con,d_v_con)

  • LMDZ5/trunk/libf/phylmd/phys_output_mod.F90

    r1724 r1742  
    190190       ctrl_out((/ 1, 6, 10, 10, 10, 10 /),'tsol_sic') /)
    191191
    192   type(ctrl_out),save,dimension(4) :: o_evappot_srf     = (/ ctrl_out((/ 1, 6, 10, 10, 10, 10 /),'evappot_ter'), &
     192  type(ctrl_out),save,dimension(4) :: o_evappot_srf  = (/ ctrl_out((/ 1, 6, 10, 10, 10, 10 /),'evappot_ter'), &
    193193       ctrl_out((/ 4, 6, 10, 10, 10, 10 /),'evappot_lic'), &
    194194       ctrl_out((/ 4, 6, 10, 10, 10, 10 /),'evappot_oce'), &
    … …  
    480480  type(ctrl_out),save :: o_rneb         = ctrl_out((/ 2, 5, 10, 10, 10, 10 /),'rneb')
    481481  type(ctrl_out),save :: o_rnebcon      = ctrl_out((/ 2, 5, 10, 10, 10, 10 /),'rnebcon')
     482  type(ctrl_out),save :: o_rnebls       = ctrl_out((/ 2, 5, 10, 10, 10, 10 /),'rnebls')
    482483  type(ctrl_out),save :: o_rhum         = ctrl_out((/ 2, 5, 10, 10, 10, 10 /),'rhum')
    483484  type(ctrl_out),save :: o_ozone        = ctrl_out((/ 2, 10, 10, 10, 10, 10 /),'ozone')
    … …  
    539540  type(ctrl_out),save :: o_wake_deltaq  = ctrl_out((/ 4, 5, 10, 10, 10, 10 /),'wake_deltaq')
    540541  type(ctrl_out),save :: o_wake_omg     = ctrl_out((/ 4, 5, 10, 10, 10, 10 /),'wake_omg')
     542  type(ctrl_out),save :: o_wdtrainA     = ctrl_out((/ 4, 1, 10,  4,  1, 10 /),'wdtrainA') !<<RomP
     543  type(ctrl_out),save :: o_wdtrainM     = ctrl_out((/ 4, 1, 10,  4,  1, 10 /),'wdtrainM') !<<RomP
    541544  type(ctrl_out),save :: o_Vprecip      = ctrl_out((/ 10, 10, 10, 10, 10, 10 /),'Vprecip')
    542545  type(ctrl_out),save :: o_ftd          = ctrl_out((/ 4, 5, 10, 10, 10, 10 /),'ftd')
    … …  
    545548  type(ctrl_out),save :: o_dtlschr      = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtlschr')
    546549  type(ctrl_out),save :: o_dqlsc        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqlsc')
     550  type(ctrl_out),save :: o_beta_prec    = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'beta_prec')
    547551  type(ctrl_out),save :: o_dtvdf        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dtvdf')
    548552  type(ctrl_out),save :: o_dqvdf        = ctrl_out((/ 4, 10, 10, 10, 10, 10 /),'dqvdf')
    … …  
    658662    USE infotrac
    659663    USE ioipsl
    660 !    USE phys_cal_mod, only : hour
     664    USE phys_cal_mod, only : hour
    661665    USE mod_phys_lmdz_para
    662666    USE aero_mod, only : naero_spc,name_aero
    … …  
    845849
    846850          idayref = day_ref
    847           CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)       
     851!          CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)       
     852! correction pour l heure initiale                               !jyg
     853!                                                                !jyg
     854          CALL ymds2ju(annee_ref, 1, idayref, hour, zjulian)         !jyg
     855! correction pour l heure initiale                               !jyg
     856!                                                                !jyg
     857!!!      CALL ymds2ju(annee_ref, 1, idayref, 0.0, zjulian)       !jyg
    848858! correction pour l heure initiale                               !jyg
    849859!                                                                !jyg
    … …  
    10471057                  o_tsol_srf(nsrf)%flag,o_tsol_srf(nsrf)%name,"Temperature "//clnsurf(nsrf),"K")
    10481058             CALL histdef2d(iff,clef_stations(iff), &
    1049                   o_evappot_srf(nsrf)%flag,o_evappot_srf(nsrf)%name,"Temperature "//clnsurf(nsrf),"K")
     1059                  o_evappot_srf(nsrf)%flag,o_evappot_srf(nsrf)%name,"Temperature"//clnsurf(nsrf),"K")
    10501060             CALL histdef2d(iff,clef_stations(iff), &
    10511061                  o_ustar_srf(nsrf)%flag,o_ustar_srf(nsrf)%name,"Friction velocity "//clnsurf(nsrf),"m/s")
    … …  
    14621472               o_rnebcon%flag,o_rnebcon%name, "Convective Cloud Fraction", "-")
    14631473          CALL histdef3d(iff,clef_stations(iff), &
     1474               o_rnebls%flag,o_rnebls%name, "LS Cloud fraction", "-")
     1475          CALL histdef3d(iff,clef_stations(iff), &
    14641476               o_rhum%flag,o_rhum%name, "Relative humidity", "-")
    14651477          CALL histdef3d(iff,clef_stations(iff), &
    … …  
    15671579                CALL histdef3d(iff,clef_stations(iff),o_wake_omg%flag,o_wake_omg%name, "wake_omg", " ")
    15681580             ENDIF
    1569              CALL histdef3d(iff,clef_stations(iff),o_Vprecip%flag,o_Vprecip%name, "precipitation vertical profile", "-")
     1581!!! RomP             CALL histdef3d(iff,clef_stations(iff),o_Vprecip%flag,o_Vprecip%name, "precipitation vertical profile", "-")
    15701582             CALL histdef3d(iff,clef_stations(iff),o_ftd%flag,o_ftd%name, "tend temp due aux descentes precip", "-")
    15711583             CALL histdef3d(iff,clef_stations(iff),o_fqd%flag,o_fqd%name,"tend vap eau due aux descentes precip", "-")
    15721584          ENDIF !(iflag_con.EQ.3)
     1585
     1586          IF(iflag_con.GE.3) THEN   !  RomP >>>
     1587            CALL histdef3d(iff,clef_stations(iff),o_wdtrainA%flag,o_wdtrainA%name, "precipitation from AA", "-")
     1588            CALL histdef3d(iff,clef_stations(iff),o_wdtrainM%flag,o_wdtrainM%name, "precipitation from mixture", "-")
     1589            CALL histdef3d(iff,clef_stations(iff),o_Vprecip%flag,o_Vprecip%name, "precipitation vertical profile", "-")
     1590          ENDIF !(iflag_con.GE.3)   ! <<< RomP
    15731591
    15741592!!! nrlmd le 10/04/2012
    … …  
    15971615          CALL histdef3d(iff,clef_stations(iff),o_dtlschr%flag,o_dtlschr%name,"Large-scale condensational heating rate","K/s")
    15981616          CALL histdef3d(iff,clef_stations(iff),o_dqlsc%flag,o_dqlsc%name, "Condensation dQ", "(kg/kg)/s")
     1617          CALL histdef3d(iff,clef_stations(iff),o_beta_prec%flag,o_beta_prec%name, "LS Conversion rate to prec", "(kg/kg)/s")
    15991618          CALL histdef3d(iff,clef_stations(iff),o_dtvdf%flag,o_dtvdf%name, "Boundary-layer dT", "K/s")
    16001619          CALL histdef3d(iff,clef_stations(iff),o_dqvdf%flag,o_dqvdf%name, "Boundary-layer dQ", "(kg/kg)/s")

  • LMDZ5/trunk/libf/phylmd/phys_output_write.h

    r1724 r1742  
    916916     $o_fqd%name,itau_w,fqd)
    917917        ENDIF
    918       ENDIF !(iflag_con.EQ.3)
     918
     919      ELSEIF (iflag_con.EQ.30) THEN
     920! sortie RomP convection descente insaturee iflag_con=30
     921       IF (o_Vprecip%flag(iff)<=lev_files(iff)) THEN
     922      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     923     $o_Vprecip%name,itau_w,Vprecip)
     924       ENDIF
     925      IF (o_wdtrainA%flag(iff)<=lev_files(iff)) THEN
     926      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     927     $o_wdtrainA%name,itau_w,wdtrainA)
     928      ENDIF
     929      IF (o_wdtrainM%flag(iff)<=lev_files(iff)) THEN
     930      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     931     $o_wdtrainM%name,itau_w,wdtrainM)
     932      ENDIF
     933
     934      ENDIF !(iflag_con.EQ.3.or.iflag_con.EQ.30)
    919935 
    920936!!! nrlmd le 10/04/2012
    … …  
    15891605       ENDIF
    15901606
     1607       IF (o_rnebls%flag(iff)<=lev_files(iff)) THEN
     1608      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     1609     $o_rnebls%name,itau_w,rneb)
     1610       ENDIF
     1611
    15911612       IF (o_rhum%flag(iff)<=lev_files(iff)) THEN
    15921613      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    … …  
    17831804      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
    17841805     $o_dqlsc%name,itau_w,zx_tmp_fi3d)
     1806       ENDIF
     1807
     1808       IF (o_beta_prec%flag(iff)<=lev_files(iff)) THEN
     1809      zx_tmp_fi3d(1:klon,1:klev)=beta_prec(1:klon,1:klev)
     1810      CALL histwrite_phy(nid_files(iff),clef_stations(iff),
     1811     $o_beta_prec%name,itau_w,zx_tmp_fi3d)
    17851812       ENDIF
    17861813

  • LMDZ5/trunk/libf/phylmd/phys_state_var_mod.F90

    r1670 r1742  
    5252      REAL, ALLOCATABLE, SAVE :: u_ancien(:,:), v_ancien(:,:)
    5353!$OMP THREADPRIVATE(u_ancien, v_ancien)
     54!!! RomP >>>
     55      REAL, ALLOCATABLE, SAVE :: tr_ancien(:,:,:)
     56!$OMP THREADPRIVATE(tr_ancien)
     57!!! RomP <<<
    5458      LOGICAL, SAVE :: ancien_ok
    5559!$OMP THREADPRIVATE(ancien_ok)
    … …  
    358362USE control_mod
    359363use aero_mod
     364use infotrac, ONLY : nbtr
    360365IMPLICIT NONE
    361366
    … …  
    384389      ALLOCATE(t_ancien(klon,klev), q_ancien(klon,klev))
    385390      ALLOCATE(u_ancien(klon,klev), v_ancien(klon,klev))
     391!!! Rom P >>>
     392      ALLOCATE(tr_ancien(klon,klev,nbtr))
     393!!! Rom P <<<
    386394      ALLOCATE(clwcon(klon,klev),rnebcon(klon,klev))
    387395      ALLOCATE(ratqs(klon,klev))
    … …  
    521529      deallocate(rugoro, t_ancien, q_ancien, clwcon, rnebcon)
    522530      deallocate(        u_ancien, v_ancien                 )
     531      deallocate(        tr_ancien)                           !RomP
    523532      deallocate(ratqs, pbl_tke)
    524533      deallocate(zmax0, f0)

  • LMDZ5/trunk/libf/phylmd/physiq.F

    r1740 r1742  
    227227      REAL d_qx(klon,klev,nqtot)
    228228      REAL d_ps(klon)
     229! Variables pour le transport convectif
    229230      real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
     231! Variables pour le lessivage convectif
     232! RomP >>>
     233      real phi2(klon,klev,klev)
     234      real d1a(klon,klev),dam(klon,klev)
     235      real ev(klon,klev),ep(klon,klev)
     236      real clw(klon,klev),elij(klon,klev,klev)
     237      real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
     238      real wdtrainA(klon,klev),wdtrainM(klon,klev)
     239! RomP <<<
    230240!IM definition dynamique o_trac dans phys_output_open
    231241!      type(ctrl_out) :: o_trac(nqtot)
    … …  
    697707      REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
    698708      REAL frac_nucl(klon,klev) ! idem (nucleation)
     709! RomP >>>
     710      REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
     711      REAL beta_prec(klon,klev)       ! taux de conv de l'eau cond (utilise)
     712! RomP <<<
    699713      INTEGER       :: iii
    700714      REAL          :: calday
    … …  
    17471761      mp(:,:)=0.
    17481762      phi(:,:,:)=0.
     1763! RomP >>>
     1764      phi2(:,:,:)=0.
     1765      beta_prec_fisrt(:,:)=0.
     1766      beta_prec(:,:)=0.
     1767      epmlmMm(:,:,:)=0.
     1768      eplaMm(:,:)=0.
     1769      d1a(:,:)=0.
     1770      dam(:,:)=0.
     1771! RomP <<<
     1772
    17491773c
    17501774c Ne pas affecter les valeurs entrees de u, v, h, et q
    … …  
    18121836         ENDDO
    18131837         ENDDO
     1838!!! RomP >>>   td dyn traceur
     1839       IF (nqtot.GE.3) THEN
     1840          DO iq = 3, nqtot
     1841          DO k = 1, klev
     1842          DO i = 1, klon
     1843            d_tr_dyn(i,k,iq-2)=
     1844     $       (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
     1845!         iiq=niadv(iq)
     1846!         print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
     1847          ENDDO
     1848          ENDDO
     1849          ENDDO
     1850       ENDIF
     1851!!! RomP <<<
    18141852      ELSE
    18151853         DO k = 1, klev
    … …  
    18211859         ENDDO
    18221860         ENDDO
     1861!!! RomP >>>   td dyn traceur
     1862        IF (nqtot.GE.3) THEN
     1863          DO iq = 3, nqtot
     1864          DO k = 1, klev
     1865          DO i = 1, klon
     1866            d_tr_dyn(i,k,iq-2)= 0.0
     1867          ENDDO
     1868          ENDDO
     1869          ENDDO
     1870       ENDIF
     1871!!! RomP <<<
    18231872         ancien_ok = .TRUE.
    18241873      ENDIF
    … …  
    23002349     .        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
    23012350     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
    2302      .        pmflxr,pmflxs,da,phi,mp,
    2303      .        ftd,fqd,lalim_conv,wght_th)
     2351! RomP >>>
     2352!!     .        pmflxr,pmflxs,da,phi,mp,
     2353!!     .        ftd,fqd,lalim_conv,wght_th)
     2354     .        pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij,
     2355     .        ftd,fqd,lalim_conv,wght_th,
     2356     .        ev, ep,epmlmMm,eplaMm,
     2357     .        wdtrainA,wdtrainM)
     2358! RomP <<<
    23042359
    23052360cIM begin
    … …  
    27832838     .           rain_lsc, snow_lsc,
    27842839     .           pfrac_impa, pfrac_nucl, pfrac_1nucl,
    2785      .           frac_impa, frac_nucl,
     2840     .           frac_impa, frac_nucl, beta_prec_fisrt,
    27862841     .           prfl, psfl, rhcl,
    27872842     .           zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
    … …  
    36383693     I     itap,     days_elapsed+1,    jH_cur,   debut,
    36393694     I     lafin,    dtime,     u, v,     t,
    3640      I     paprs,    pplay,     pmfu,     pmfd, 
     3695     I     paprs,    pplay,     pmfu,     pmfd,
    36413696     I     pen_u,    pde_u,     pen_d,    pde_d,
    36423697     I     cdragh,   coefh,     fm_therm, entr_therm,
    36433698     I     u1,       v1,        ftsol,    pctsrf,
    36443699     I     ustar,     u10m,      v10m,
    3645      I     rlat,     frac_impa, frac_nucl,rlon,
     3700     I     rlat,     rlon,
     3701     I     frac_impa,frac_nucl, beta_prec_fisrt,beta_prec,
    36463702     I     presnivs, pphis,     pphi,     albsol1,
    3647      I     sh_in,    rhcl,      cldfra,   rneb, 
     3703     I     sh_in,    rhcl,      cldfra,   rneb,
    36483704     I     diafra,   cldliq,    itop_con, ibas_con,
    36493705     I     pmflxr,   pmflxs,    prfl,     psfl,
    3650      I     da,       phi,       mp,       upwd,     
     3706     I     da,       phi,       mp,       upwd,
     3707     I     phi2,     d1a,       dam,      sij,         !<<RomP
     3708     I     wdtrainA, wdtrainM,  sigd,     clw,elij,    !<<RomP
     3709     I     ev,       ep,        epmlmMm,  eplaMm,      !<<RomP
    36513710     I     dnwd,     aerosol_couple,      flxmass_w,
    36523711     I     tau_aero, piz_aero,  cg_aero,  ccm,
    36533712     I     rfname,
     3713     I     d_tr_dyn,                                   !<<RomP
    36543714     O     tr_seri)
    36553715
    … …  
    38233883      ENDDO
    38243884      ENDDO
    3825 c
     3885
     3886!!! RomP >>>
     3887      IF (nqtot.GE.3) THEN
     3888        DO iq = 3, nqtot
     3889        DO k = 1, klev
     3890        DO i = 1, klon
     3891           tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
     3892        ENDDO
     3893        ENDDO
     3894        ENDDO
     3895      ENDIF
     3896!!! RomP <<<
    38263897!==========================================================================
    38273898! Sorties des tendances pour un point particulier

  • LMDZ5/trunk/libf/phylmd/phytrac.F90

    r1670 r1742  
    33SUBROUTINE phytrac(                            &
    44     nstep,     julien,   gmtime,   debutphy,  &
    5      lafin,     pdtphys,  u, v,     t_seri,     &
     5     lafin,     pdtphys,  u, v,     t_seri,    &
    66     paprs,     pplay,    pmfu,     pmfd,      &
    77     pen_u,     pde_u,    pen_d,    pde_d,     &
    … …  
    99     yu1,       yv1,      ftsol,    pctsrf,    &
    1010     ustar,     u10m,      v10m,               &
    11      xlat,      frac_impa,frac_nucl,xlon,      &
     11     xlat,      xlon,                          &
     12     frac_impa,frac_nucl,beta_fisrt,beta_v1,   &
    1213     presnivs,  pphis,    pphi,     albsol,    &
    1314     sh,        rh,       cldfra,   rneb,      &
    … …  
    1516     pmflxr,    pmflxs,   prfl,     psfl,      &
    1617     da,        phi,      mp,       upwd,      &
     18     phi2,      d1a,      dam,      sij,       &   ! RomP
     19     wdtrainA,  wdtrainM, sigd,     clw,elij,  &   ! RomP
     20     evap,      ep,       epmlmMm,  eplaMm,    &   ! RomP
    1721     dnwd,      aerosol_couple,     flxmass_w, &
    1822     tau_aero,  piz_aero,  cg_aero, ccm,       &
    1923     rfname,                                   &
     24     d_tr_dyn,                                 &   ! RomP
    2025     tr_seri)         
    2126!
    … …  
    2328! Auteur(s) FH
    2429! Objet: Moniteur general des tendances traceurs
     30! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
     31! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
    2532!======================================================================
    2633
    2734  USE ioipsl
     35  USE phys_cal_mod, only : hour
     36  USE phys_output_mod, only : convers_timesteps
    2837  USE dimphy
    2938  USE infotrac
    … …  
    3645  USE tracreprobus_mod
    3746  USE control_mod
    38 
    3947
    4048  IMPLICIT NONE
    … …  
    6876!--------
    6977  REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri  ! Temperature
    70   REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used
    71   REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used
     78  REAL,DIMENSION(klon,klev),INTENT(IN)   :: u       ! variable not used 
     79  REAL,DIMENSION(klon,klev),INTENT(IN)   :: v       ! variable not used 
    7280  REAL,DIMENSION(klon,klev),INTENT(IN)   :: sh      ! humidite specifique
    7381  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rh      ! humidite relative
    … …  
    8189  REAL,DIMENSION(klon,klev),INTENT(IN)   :: diafra  ! fraction nuageuse (convection ou stratus artificiels)
    8290  REAL,DIMENSION(klon,klev),INTENT(IN)   :: rneb    ! fraction nuageuse (grande echelle)
     91!
     92  REAL                                   :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
     93  REAL,DIMENSION(klon,klev),INTENT(IN)   :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
     94  REAL,DIMENSION(klon,klev),INTENT(out)  :: beta_v1    ! -- (originale version)
     95
     96!
    8397  INTEGER,DIMENSION(klon),INTENT(IN)     :: itop_con
    8498  INTEGER,DIMENSION(klon),INTENT(IN)     :: ibas_con
    8599  REAL,DIMENSION(klon),INTENT(IN)        :: albsol  ! albedo surface
     100!
     101!Dynamique
     102!--------
     103  REAL,DIMENSION(klon,klev,nbtr),INTENT(IN)    :: d_tr_dyn
    86104!
    87105!Convection:
    … …  
    108126  REAL,DIMENSION(klon,klev),INTENT(IN)     :: da
    109127  REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
     128! RomP >>>
     129  REAL,DIMENSION(klon,klev),INTENT(IN)      :: d1a,dam
     130  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
     131!
     132  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainA
     133  REAL,DIMENSION(klon,klev),INTENT(IN)      :: wdtrainM
     134  REAL,DIMENSION(klon),INTENT(IN)           :: sigd
     135! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
     136  REAL,DIMENSION(klon,klev),INTENT(IN)      :: evap
     137  REAL,DIMENSION(klon,klev),INTENT(IN)      :: ep
     138  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
     139  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
     140  REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm
     141  REAL,DIMENSION(klon,klev),INTENT(IN)      :: eplaMm
     142  REAL,DIMENSION(klon,klev),INTENT(IN)      :: clw
     143! RomP <<<
     144
     145!
    110146  REAL,DIMENSION(klon,klev),INTENT(IN)     :: mp
    111147  REAL,DIMENSION(klon,klev),INTENT(IN)     :: upwd      ! saturated updraft mass flux
    … …  
    120156!--------------
    121157!
    122   REAL,DIMENSION(klon),INTENT(IN)     :: cdragh ! coeff drag pour T et Q
    123   REAL,DIMENSION(klon,klev),INTENT(IN):: coefh  ! coeff melange CL (m**2/s)
    124   REAL,DIMENSION(klon),INTENT(IN)     :: ustar,u10m,v10m ! u* & vent a 10m (m/s)
    125   REAL,DIMENSION(klon),INTENT(IN)     :: yu1    ! vents au premier niveau
    126   REAL,DIMENSION(klon),INTENT(IN)     :: yv1    ! vents au premier niveau
     158  REAL,DIMENSION(klon),INTENT(IN)      :: cdragh ! coeff drag pour T et Q
     159  REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh  ! coeff melange CL (m**2/s)
     160  REAL,DIMENSION(klon),INTENT(IN)      :: ustar,u10m,v10m ! u* & vent a 10m (m/s)
     161  REAL,DIMENSION(klon),INTENT(IN)      :: yu1    ! vents au premier niveau
     162  REAL,DIMENSION(klon),INTENT(IN)      :: yv1    ! vents au premier niveau
    127163!
    128164!Lessivage:
    … …  
    141177! Output argument
    142178!----------------
    143   REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] 
    144 
     179  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
     180  REAL,DIMENSION(klon,klev)                    :: sourceBE
    145181!=======================================================================================
    146182!                        -- LOCAL VARIABLES --
    … …  
    153189!--------------------------------------------
    154190!
    155   REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit 
     191  REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source  ! a voir lorsque le flux de surface est prescrit
    156192!$OMP THREADPRIVATE(source)
    157193
    … …  
    166202  REAL,DIMENSION(klon)      :: zx_tmp_fi2d ! variable temporaire grille physique
    167203  INTEGER                   :: itau_w      ! pas de temps ecriture = nstep + itau_phy
    168   LOGICAL,PARAMETER :: ok_sync=.TRUE.
     204  LOGICAL,PARAMETER         :: ok_sync=.TRUE.
     205  CHARACTER(len=20)         :: chtratimestep
    169206
    170207!
    … …  
    175212  REAL,DIMENSION(klon,klev)             :: delp     ! epaisseur de couche (Pa)
    176213!
    177 ! Tendances de traceurs (Td):
     214! Tendances de traceurs (Td) et flux de traceurs:
    178215!------------------------
    179 !
    180216  REAL,DIMENSION(klon,klev)      :: d_tr     ! Td dans l'atmosphere
    181217  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cl  ! Td couche limite/traceur
     218  REAL,DIMENSION(klon,nbtr)      :: d_tr_dry ! Td depot sec/traceur (1st layer)  jyg
     219  REAL,DIMENSION(klon,nbtr)      :: flux_tr_dry ! depot sec/traceur (surface)    jyg
     220  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_dec                            !RomP
    182221  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_cv  ! Td convection/traceur
     222! RomP >>>
     223  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_insc
     224  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_bcscav
     225  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_evapls
     226  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_ls
     227  REAL,DIMENSION(klon,nbtr)      :: qPrls      !jyg: concentration tra dans pluie LS a la surf.
     228  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_trsp
     229  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_sscav
     230  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_sat
     231  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_uscav
     232  REAL,DIMENSION(klon,klev,nbtr) :: qPr,qDi ! concentration tra dans pluie,air descente insaturee
     233  REAL,DIMENSION(klon,klev,nbtr) :: qPa,qMel
     234  REAL,DIMENSION(klon,klev,nbtr) :: qTrdi,dtrcvMA ! conc traceur descente air insaturee et td convective MA
     235  REAL,DIMENSION(klon,klev)      :: Mint
     236  REAL,DIMENSION(klon,klev,nbtr) :: zmfd1a
     237  REAL,DIMENSION(klon,klev,nbtr) :: zmfdam
     238  REAL,DIMENSION(klon,klev,nbtr) :: zmfphi2
     239! RomP <<<
    183240  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_th  ! Td thermique
    184241  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_impa ! Td du lessivage par impaction
    185   REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation 
     242  REAL,DIMENSION(klon,klev,nbtr) :: d_tr_lessi_nucl ! Td du lessivage par nucleation
    186243!
    187244! Physique
    188 !----------  
     245!----------
    189246  REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
    190247  REAL,DIMENSION(klon,klev)      :: zmasse    ! densité atmosphérique Kg/m2
    191248  REAL,DIMENSION(klon,klev)      :: ztra_th
     249!PhH
     250  REAL,DIMENSION(klon,klev)      :: zrho
     251  REAL,DIMENSION(klon,klev)      :: zdz
     252  REAL                           :: evaplsc,dx,beta ! variable pour lessivage Genthon
     253  REAL,DIMENSION(klon)           :: his_dh          ! ---
     254! in-cloud scav variables
     255  REAL           :: ql_incloud_ref     ! ref value of in-cloud condensed water content
    192256 
    193257!Controles:
    … …  
    199263
    200264  CHARACTER(len=8),DIMENSION(nbtr) :: solsym
    201 
    202 
     265!RomP >>>
     266  INTEGER,SAVE  :: iflag_lscav
     267  LOGICAL,SAVE  :: convscav
     268!$OMP THREADPRIVATE(iflag_lscav,convscav)
     269!RomP <<<
    203270!######################################################################
    204271!                    -- INITIALIZATION --
    205272!######################################################################
     273  DO k=1,klev
     274     DO i=1,klon
     275      sourceBE(i,k)=0.
     276      Mint(i,k)=0.
     277      zrho(i,k)=0.
     278      zdz(i,k)=0.
     279     END DO
     280  END DO
     281
     282  DO it=1, nbtr
     283   DO k=1,klev
     284    DO i=1,klon
     285    d_tr_insc(i,k,it)=0.
     286    d_tr_bcscav(i,k,it)=0.
     287    d_tr_evapls(i,k,it)=0.
     288    d_tr_ls(i,k,it)=0.
     289    d_tr_cv(i,k,it)=0.
     290    d_tr_cl(i,k,it)=0.
     291    d_tr_trsp(i,k,it)=0.
     292    d_tr_sscav(i,k,it)=0.
     293    d_tr_sat(i,k,it)=0.
     294    d_tr_uscav(i,k,it)=0.
     295    d_tr_lessi_impa(i,k,it)=0.
     296    d_tr_lessi_nucl(i,k,it)=0.
     297    qDi(i,k,it)=0.
     298    qPr(i,k,it)=0.
     299    qPa(i,k,it)=0.
     300    qMel(i,k,it)=0.
     301    qTrdi(i,k,it)=0.
     302    dtrcvMA(i,k,it)=0.
     303    zmfd1a(i,k,it)=0.
     304    zmfdam(i,k,it)=0.
     305    zmfphi2(i,k,it)=0.
     306    END DO
     307   END DO
     308  END DO
    206309  IF (debutphy) THEN
    207      IF (prt_level >9) WRITE(lunout,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
     310!!jyg
     311   chtratimestep='DefFreq'
     312   CALL getin('tra_time_step',chtratimestep)
     313   IF (chtratimestep .NE. 'DefFreq') THEN
     314     call convers_timesteps(chtratimestep,pdtphys,ecrit_tra)
     315   ENDIF
     316!RomP >>>
     317!
     318!Config Key  = convscav
     319!Config Desc = Convective scavenging switch: 0=off, 1=on.
     320!Config Def  = .false.
     321!Config Help =
     322!
     323  convscav=.false.
     324  call getin('convscav', convscav)
     325  print*,'phytrac passage dans routine conv avec lessivage', convscav
     326!
     327!Config Key  = iflag_lscav
     328!Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
     329!              2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
     330!Config Def  = 1
     331!Config Help =
     332!
     333  iflag_lscav=1
     334  call getin('iflag_lscav', iflag_lscav)
     335!
     336  SELECT CASE(iflag_lscav)
     337  CASE(0)
     338   PRINT*, 'Large scale scavenging: none'
     339  CASE(1)
     340   PRINT*, 'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
     341  CASE(2)
     342   PRINT*, 'Large scale scavenging: C. Genthon, modified P. Heinrich'
     343  CASE(3)
     344   PRINT*, 'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
     345  CASE(4)
     346   PRINT*, 'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
     347  END SELECT
     348!RomP <<<
     349     WRITE(*,*) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
    208350     ALLOCATE( source(klon,nbtr), stat=ierr)
    209351     IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
    … …  
    237379     END DO
    238380  END DO
     381!
     382  IF (id_be .GT. 0) THEN
     383  DO k=1,klev
     384     DO i=1,klon
     385        sourceBE(i,k)=srcbe(i,k)       !RomP  -> pour sortie histrac
     386     END DO
     387  END DO
     388  ENDIF
    239389
    240390!===============================================================================
    … …  
    246396     !    -- Traitement des traceurs avec traclmdz
    247397     CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
    248           cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,couchelimite,sh,&
     398          cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,couchelimite,sh, &
    249399          rh, pphi, ustar, u10m, v10m, &
    250           tr_seri, source, solsym, d_tr_cl, zmasse)
     400!!          tr_seri, source, solsym, d_tr_cl, zmasse)                      !RomP
     401          tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse)               !RomP
    251402  CASE('inca')
    252403     !    -- CHIMIE INCA  config_inca = aero or chem --
    … …  
    273424     
    274425  END SELECT
    275 
    276426!======================================================================
    277427!       -- Calcul de l'effet de la convection --
    278428!======================================================================
     429
    279430  IF (convection) THEN
    280431     DO it=1, nbtr
    281432        IF ( conv_flg(it) == 0 ) CYCLE
    282        
    283433        IF (iflag_con.LT.2) THEN
    284            d_tr_cv(:,:,:)=0.
     434           d_tr_cv(:,:,it)=0.
    285435        ELSE IF (iflag_con.EQ.2) THEN
    286436!..Tiedke
    287437           CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
    288438                pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
    289         ELSE
    290 !..K.Emanuel
    291            CALL cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(:,:,it),&
    292                 upwd,dnwd,d_tr_cv(:,:,it))
     439! RomP >>>               
     440        ELSE   
     441!..K.Emanuel                  !RomP modif arg
     442        if (convscav.and.aerosol(it)) then    ! lessivage convectif pour aerosol
     443!
     444          CALL cvltr(pdtphys, da, phi,phi2,d1a,dam, mp,ep,              &
     445               sigd,sij,clw,elij,epmlmMm,eplaMm,                        &
     446               pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM,             &   
     447               paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con,            &
     448               d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,&
     449               qPa,qMel,qTrdi,dtrcvMA,Mint,                             &
     450               zmfd1a,zmfphi2,zmfdam)
     451        else  !pas de lessivage convectif ou n'est pas un aerosol
     452           CALL cvltrorig(it,pdtphys, da, phi,mp,paprs,pplay,tr_seri,&
     453                    upwd,dnwd,d_tr_cv)
     454        endif
    293455        END IF
     456! RomP <<<
    294457
    295458        DO k = 1, klev
    … …  
    357520                tr_seri(:,:,it), source(:,it),      &
    358521                paprs, pplay, delp,                 &
    359                 d_tr_cl(:,:,it))
     522                d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
    360523           
    361524           DO k = 1, klev
    … …  
    372535
    373536!======================================================================
    374 !   Calcul de l'effet de la precipitation
    375 !======================================================================
    376 
     537!   Calcul de l'effet de la precipitation grande echelle
     538!======================================================================
    377539  IF (lessivage) THEN
    378      
     540
     541   ql_incloud_ref = 10.e-4
     542   ql_incloud_ref =  5.e-4
     543
     544
     545! calcul du contenu en eau liquide au sein du nuage
     546     ql_incl = ql_incloud_ref
     547! choix du lessivage
     548!
     549  IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
     550! ********  Olivier Boucher version (3) possibly with modified ql_incl (4)
     551!
     552   DO it = 1, nbtr
     553!  incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
     554! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
     555! Liu (2001) proposed to use 1.5e-3 kg/kg
     556
     557    CALL lsc_scav(pdtphys,it,iflag_lscav,ql_incl,prfl,psfl,rneb,beta_fisrt,  &
     558                  beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc,   &
     559                  d_tr_bcscav,d_tr_evapls,qPrls)
     560
     561!large scale scavenging tendency
     562   DO k = 1, klev
     563    DO i = 1, klon
     564    d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
     565    tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
     566    ENDDO
     567   ENDDO
     568     CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//solsym(it))
     569   END DO  !tr
     570
     571 ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
     572! *********   modified  old version
     573
     574     d_tr_lessi_nucl(:,:,:) = 0.
     575     d_tr_lessi_impa(:,:,:) = 0.
     576     flestottr(:,:,:) = 0.
     577! Tendance des aerosols nuclees et impactes
     578     DO it = 1, nbtr
     579        IF (aerosol(it)) THEN
     580        his_dh(:)=0.
     581           DO k = 1, klev
     582              DO i = 1, klon
     583!PhH
     584              zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
     585              zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
     586!
     587              END DO
     588           END DO
     589
     590          DO k=klev-1, 1, -1
     591            DO i=1, klon
     592!             d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.)
     593             dx=d_tr_ls(i,k,it)
     594             his_dh(i)=his_dh(i)-dx*zrho(i,k)*zdz(i,k)/pdtphys !  kg/m2/s
     595             evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
     596! Evaporation Partielle -> Liberation Partielle 0.5*evap
     597            IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
     598                evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
     599! evaplsc est donc positif, his_dh(i) est positif
     600!-------------- 
     601             d_tr_evapls(i,k,it)=0.5*evaplsc*(d_tr_lessi_nucl(i,k+1,it) &
     602                                  +d_tr_lessi_impa(i,k+1,it))
     603!-------------   d_tr_evapls(i,k,it)=-0.5*evaplsc*(d_tr_lsc(i,k+1,it))
     604             beta=0.5*evaplsc
     605              if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
     606               beta=1.0*evaplsc
     607              endif
     608            dx=beta*his_dh(i)/zrho(i,k)/zdz(i,k)*pdtphys
     609            his_dh(i)=(1.-beta)*his_dh(i)   ! tracer from
     610            d_tr_evapls(i,k,it)=dx
     611            ENDIF
     612            d_tr_ls(i,k,it)=tr_seri(i,k,it)*(frac_impa(i,k)*frac_nucl(i,k)-1.) &
     613                            +d_tr_evapls(i,k,it)
     614
     615!--------------
     616                 d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +    &
     617                      ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
     618                 d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +    &
     619                      ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
     620!
     621! Flux lessivage total
     622                 flestottr(i,k,it) = flestottr(i,k,it) -   &
     623                      ( d_tr_lessi_nucl(i,k,it)   +        &
     624                      d_tr_lessi_impa(i,k,it) ) *          &
     625                      ( paprs(i,k)-paprs(i,k+1) ) /        &
     626                      (RG * pdtphys)
     627!! Mise a jour des traceurs due a l'impaction,nucleation
     628!                 tr_seri(i,k,it)=tr_seri(i,k,it)*frac_impa(i,k)*frac_nucl(i,k)
     629!!  calcul de la tendance liee au lessivage stratiforme
     630!                 d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
     631!                                (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
     632!--------------
     633              END DO
     634           END DO
     635        END IF
     636     END DO
     637! *********   end modified old version
     638
     639 ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
     640! *********    old version
     641
    379642     d_tr_lessi_nucl(:,:,:) = 0.
    380643     d_tr_lessi_impa(:,:,:) = 0.
    … …  
    412675     END DO
    413676     
    414   END IF ! lessivage
     677! *********   end old version
     678  ENDIF  !  iflag_lscav . EQ. 1, 2, 3 or 4
     679!
     680  END IF !  lessivage
    415681
    416682!=============================================================

  • LMDZ5/trunk/libf/phylmd/surf_land_bucket_mod.F90

    r1724 r1742  
    101101       CALL calbeta_clim(knon,jour,rlatd(knindex(:)),beta)
    102102    endif
    103 
    104103       
    105104! calculate temperature, heat capacity and conduction flux in soil
    … …  
    112111    ELSE
    113112       cal(:) = RCPD * capsol(:)
     113       IF (klon_glo .EQ. 1) THEN
     114         cal(:) = 0.
     115       ENDIF
    114116    ENDIF
    115117   

  • LMDZ5/trunk/libf/phylmd/traclmdz_mod.F90

    r1670 r1742  
    22!
    33MODULE traclmdz_mod
     4
    45!
    56! In this module all tracers specific to LMDZ are treated. This module is used
    … …  
    117118    REAL, DIMENSION(klev)          :: mintmp, maxtmp
    118119    LOGICAL                        :: zero
    119 
     120! RomP >>> profil initial Be7
     121      integer ilesfil
     122      parameter (ilesfil=1)
     123      integer  irr,kradio
     124      real     beryllium(klon,klev)
     125! profil initial Pb210
     126      integer ilesfil2
     127      parameter (ilesfil2=1)
     128      integer  irr2,kradio2
     129      real     plomb(klon,klev)
     130!! RomP <<<
    120131! --------------------------------------------
    121132! Allocation
    … …  
    148159
    149160    lessivage  = .TRUE.
     161!!jyg(20130206) : le choix d activation du lessivage est fait dans phytrac avec iflag_lscav
     162!!    call getin('lessivage',lessivage)
     163!!    if(lessivage) then
     164!!     print*,'lessivage lsc ON'
     165!!    else
     166!!     print*,'lessivage lsc OFF'
     167!!    endif
    150168    aerosol(:) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
    151169   
    … …  
    161179       ELSE IF ( tname(iiq) == "PB") THEN
    162180          id_pb=it ! plomb
     181! RomP >>> profil initial de PB210
     182     open (ilesfil2,file='prof.pb210',status='old',iostat=irr2)
     183     IF (irr2 == 0) THEN
     184      read(ilesfil2,*) kradio2
     185      print*,'number of levels for pb210 profile ',kradio2
     186      do k=kradio2,1,-1
     187       read (ilesfil2,*) plomb(:,k)
     188      enddo
     189      close(ilesfil2)
     190      do k=1,klev
     191       do i=1,klon
     192         tr_seri(i,k,id_pb)=plomb(i,k)
     193!!        print*, 'tr_seri',i,k,tr_seri(i,k,id_pb)
     194        enddo
     195      enddo
     196     ELSE
     197       print *, 'Prof.pb210 does not exist: use restart values'
     198     ENDIF
     199! RomP <<<
    163200       ELSE IF ( tname(iiq) == "Aga" .OR. tname(iiq) == "AGA" ) THEN
    164201          ! Age of stratospheric air
    … …  
    183220          radio(id_be) = .TRUE.
    184221          aerosol(id_be) = .TRUE. ! le Be est un aerosol
    185           CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
     222!jyg le 13/03/2013 ; ajout de pplay en argument de init_be
     223!!!          CALL init_be(pctsrf,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
     224          CALL init_be(pctsrf,pplay,masktr(:,id_be),tautr(id_be),vdeptr(id_be),scavtr(id_be),srcbe)
    186225          WRITE(lunout,*) 'Initialisation srcBe: OK'
     226! RomP >>> profil initial de Be7
     227      open (ilesfil,file='prof.be7',status='old',iostat=irr)
     228      IF (irr == 0) THEN
     229       read(ilesfil,*) kradio
     230       print*,'number of levels for Be7 profile ',kradio
     231       do k=kradio,1,-1
     232        read (ilesfil,*) beryllium(:,k)
     233       enddo
     234       close(ilesfil)
     235       do k=1,klev
     236         do i=1,klon
     237          tr_seri(i,k,id_be)=beryllium(i,k)
     238!!        print*, 'tr_seri',i,k,tr_seri(i,k,id_be)
     239         enddo
     240       enddo
     241     ELSE
     242       print *, 'Prof.Be7 does not exist: use restart values'
     243     ENDIF
     244! RomP <<<
    187245       ELSE IF (tname(iiq)=="O3" .OR. tname(iiq)=="o3") THEN
    188246          ! Recherche de l'ozone : parametrization de la chimie par Cariolle
    … …  
    280338       cdragh,  coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon, couchelimite, sh, &
    281339       rh, pphi, ustar, zu10m, zv10m, &
    282        tr_seri, source, solsym, d_tr_cl, zmasse)
     340!!          tr_seri, source, solsym, d_tr_cl, zmasse)                      !RomP
     341          tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse)               !RomP
    283342   
    284343    USE dimphy
    … …  
    316375!--------------
    317376!
    318     REAL,DIMENSION(klon),INTENT(IN)      :: cdragh  ! coeff drag pour T et Q
    319     REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh   ! diffusivite turb (m**2/s)
    320     REAL,DIMENSION(klon),INTENT(IN)      :: yu1     ! vents au premier niveau
    321     REAL,DIMENSION(klon),INTENT(IN)      :: yv1     ! vents au premier niveau
     377    REAL,DIMENSION(klon),INTENT(IN)      :: cdragh     ! coeff drag pour T et Q
     378    REAL,DIMENSION(klon,klev),INTENT(IN) :: coefh      ! coeff melange CL (m**2/s)
     379    REAL,DIMENSION(klon),INTENT(IN)      :: yu1        ! vents au premier niveau
     380    REAL,DIMENSION(klon),INTENT(IN)      :: yv1        ! vents au premier niveau
    322381    LOGICAL,INTENT(IN)                   :: couchelimite
    323     REAL,DIMENSION(klon,klev),INTENT(IN) :: sh      ! humidite specifique
     382    REAL,DIMENSION(klon,klev),INTENT(IN) :: sh         ! humidite specifique
    324383    REAL,DIMENSION(klon,klev),INTENT(IN) :: rh      ! Humidite relative
    325384    REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi    ! geopotentie

  • LMDZ5/trunk/libf/phylmd/write_histrac.h

    r1577 r1742  
    1010     CALL histwrite_phy(nid_tra,.FALSE.,"aire",itau_w,airephy)
    1111     CALL histwrite_phy(nid_tra,.FALSE.,"zmasse",itau_w,zmasse)
     12! RomP >>>
     13     CALL histwrite_phy(nid_tra,.FALSE.,"sourceBE",itau_w,sourceBE)
     14! RomP <<<
    1215
    1316!TRACEURS
    … …  
    2023
    2124! TD LESSIVAGE       
    22         IF (lessivage .AND. aerosol(it)) THEN
    23            CALL histwrite_phy(nid_tra,.FALSE.,"fl"//tname(iiq),itau_w,flestottr(:,:,it))
     25      IF (lessivage .AND. aerosol(it)) THEN
     26        CALL histwrite_phy(nid_tra,.FALSE.,"fl"//tname(iiq),itau_w,flestottr(:,:,it))
     27        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_ls_"//tname(iiq),itau_w,d_tr_ls(:,:,it))
     28        IF(iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) then
     29          CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_insc_"//tname(iiq),itau_w,d_tr_insc(:,:,it))
     30          CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_bcscav_"//tname(iiq),itau_w,d_tr_bcscav(:,:,it))
     31          CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_evls_"//tname(iiq),itau_w,d_tr_evapls(:,:,it))
     32          CALL histwrite_phy(nid_tra,.FALSE.,"qpr_ls_"//tname(iiq),itau_w,qPrls(:,it))
    2433        ENDIF
     34      ENDIF
    2535
    2636! TD THERMIQUES
    … …  
    3545
    3646! TD COUCHE-LIMITE
     47      IF (couchelimite) THEN
    3748        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_cl_"//tname(iiq),itau_w,d_tr_cl(:,:,it))
     49        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_dry_"//tname(iiq),itau_w,d_tr_dry(:,it))
     50        CALL histwrite_phy(nid_tra,.FALSE.,"flux_tr_dry_"//tname(iiq),itau_w,flux_tr_dry(:,it))
     51      ENDIF
     52
     53! TD radio-decroissance
     54        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_dec_"//tname(iiq),itau_w,d_tr_dec(:,:,it))
     55
     56! RomP >>>
     57        IF (iflag_con.EQ.30) THEN
     58        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_cvMA_"//tname(iiq),itau_w,dtrcvMA(:,:,it))
     59        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_trsp_"//tname(iiq),itau_w,d_tr_trsp(:,:,it))
     60        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_sscav_"//tname(iiq),itau_w,d_tr_sscav(:,:,it))
     61        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_sat_"//tname(iiq),itau_w,d_tr_sat(:,:,it))
     62        CALL histwrite_phy(nid_tra,.FALSE.,"d_tr_uscav_"//tname(iiq),itau_w,d_tr_uscav(:,:,it))
     63        CALL histwrite_phy(nid_tra,.FALSE.,"tr_pr_"//tname(iiq),itau_w,qPr(:,:,it))
     64        CALL histwrite_phy(nid_tra,.FALSE.,"tr_aa_"//tname(iiq),itau_w,qPa(:,:,it))
     65        CALL histwrite_phy(nid_tra,.FALSE.,"tr_mel_"//tname(iiq),itau_w,qMel(:,:,it))
     66        CALL histwrite_phy(nid_tra,.FALSE.,"tr_di_"//tname(iiq),itau_w,qDi(:,:,it))
     67        CALL histwrite_phy(nid_tra,.FALSE.,"tr_trspdi_"//tname(iiq),itau_w,qTrdi(:,:,it))
     68        CALL histwrite_phy(nid_tra,.FALSE.,"zmfd1a_"//tname(iiq),itau_w,zmfd1a(:,:,it))
     69        CALL histwrite_phy(nid_tra,.FALSE.,"zmfphi2_"//tname(iiq),itau_w,zmfphi2(:,:,it))
     70        CALL histwrite_phy(nid_tra,.FALSE.,"zmfdam_"//tname(iiq),itau_w,zmfdam(:,:,it))
     71        ENDIF
     72        CALL histwrite_phy(nid_tra,.FALSE.,"dtrdyn_"//tname(iiq),itau_w,d_tr_dyn(:,:,it))
     73! RomP <<<
    3874     ENDDO
    3975!---------------
    … …  
    65101 
    66102! DIVERS   
     103     CALL histwrite_phy(nid_tra,.FALSE.,"Mint",itau_w,Mint(:,:))
     104     CALL histwrite_phy(nid_tra,.FALSE.,"frac_impa",itau_w,frac_impa(:,:))   
     105     CALL histwrite_phy(nid_tra,.FALSE.,"frac_nucl",itau_w,frac_nucl(:,:))
     106
     107
    67108     CALL histwrite_phy(nid_tra,.FALSE.,"pplay",itau_w,pplay)     
    68109     CALL histwrite_phy(nid_tra,.FALSE.,"T",itau_w,t_seri)     
Note: See TracChangeset for help on using the changeset viewer.