Changeset 879


Ignore:
Timestamp:
Jan 17, 2008, 10:20:32 PM (17 years ago)
Author:
Laurent Fairhead
Message:

Suite de la bascule vers une physique avec thermiques, nouvelle convection, poche froide ...
LF

Location:
LMDZ4/trunk/libf/phylmd
Files:
19 added
3 deleted
21 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/phylmd/aaam_bud.F

    r644 r879  
     1!
     2! $Header$
     3!
    14      subroutine aaam_bud (iam,nlon,nlev,rjour,rsec,
    25     i                   rea,rg,ome,     
    … …  
    210213C South Pole
    211214
     215      if (jjm.GT.1) then
    212216      l=l+1
    213217      ub(1,jjm+1)=0.
    … …  
    229233      vb(i,jjm+1)=vb(1,jjm+1)                               
    230234      enddo
     235      endif
    231236
    232237C

  • LMDZ4/trunk/libf/phylmd/calltherm.F90

    r878 r879  
    77     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs  &
    88     &      ,fm_therm,entr_therm,zqasc,clwcon0,lmax,ratqscth,  &
    9      &       ratqsdiff,zqsatth)
     9     &       ratqsdiff,zqsatth,Ale_bl,Alp_bl,lalim_conv,wght_th)
    1010
    1111      implicit none
    … …  
    2727      REAL pphi(klon,klev)
    2828      real zlev(klon,klev+1)
     29!test: on sort lentr et a* pour alimenter KE
     30      REAL wght_th(klon,klev)
     31      INTEGER lalim_conv(klon)
    2932
    3033!FH Update Thermiques
    … …  
    5053      real ratqsdiff(klon,klev)
    5154      real zqsatth(klon,klev) 
     55!nouvelles variables pour la convection
     56      real Ale_bl(klon)
     57      real Alp_bl(klon)
     58      real Ale(klon)
     59      real Alp(klon)
     60!RC
    5261!********************************************************
    5362
    … …  
    7180         fm_therm(:,:)=0.
    7281         entr_therm(:,:)=0.
     82         Ale_bl(:)=0.
     83         Alp_bl(:)=0.
    7384       print*,'thermV4 nsplit: ',nsplit_thermals,' weak_inversion'
    7485
    … …  
    160171     &      ,ratqscth,ratqsdiff,zqsatth  &
    161172     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
    162      &      ,tho_thermals)
     173     &      ,tho_thermals,Ale,Alp,lalim_conv,wght_th)
    163174         endif
    164175
    … …  
    195206           ENDIF
    196207       ENDDO
     208       ENDDO
     209
     210       DO i=1,klon
     211            fm_therm(i,klev+1)=0.
     212            Ale_bl(i)=Ale_bl(i)+Ale(i)/float(nsplit_thermals)
     213!            write(22,*)'ALE CALLTHERM',Ale_bl(i),Ale(i)
     214            Alp_bl(i)=Alp_bl(i)+Alp(i)/float(nsplit_thermals)
     215!            write(23,*)'ALP CALLTHERM',Alp_bl(i),Alp(i)
    197216       ENDDO
    198217

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

    r766 r879  
    22! $Header$
    33!
    4       SUBROUTINE concvl (iflag_con,dtime,paprs,pplay,t,q,u,v,tra,ntra,
    5      .             work1,work2,d_t,d_q,d_u,d_v,d_tra,
    6      .             rain, snow, kbas, ktop,
    7      .             upwd,dnwd,dnwdbis,Ma,cape,tvp,iflag,
     4      SUBROUTINE concvl (iflag_con,iflag_clos,
     5     .             dtime,paprs,pplay,
     6     .             t,q,t_wake,q_wake,u,v,tra,ntra,
     7     .             ALE,ALP,work1,work2,
     8     .             d_t,d_q,d_u,d_v,d_tra,
     9     .             rain, snow, kbas, ktop, sigd,
     10     .             upwd,dnwd,dnwdbis,Ma,mip,Vprecip,
     11     .             cape,cin,tvp,Tconv,iflag,
    812     .             pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,
    9      .             qcondc,wd,
    10      .             pmflxr,pmflxs,
    11      .             da,phi,mp)
    12  
    13 c
    14       USE dimphy
     13     .             qcondc,wd,pmflxr,pmflxs,
     14     .             da,phi,mp,dd_t,dd_q,lalim_conv,wght_th)
     15***************************************************************
     16*                                                             *
     17* CONCVL                                                      *
     18*                                                             *
     19*                                                             *
     20* written by   : Sandrine Bony-Lena, 17/05/2003, 11.16.04    *
     21* modified by :                                               *
     22***************************************************************
     23*
     24c
     25c      USE dimphy
    1526      IMPLICIT none
    1627c======================================================================
    17 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930818
     28c Auteur(s): S. Bony-Lena (LMD/CNRS) date: ???
    1829c Objet: schema de convection de Emanuel (1991) interface
    1930c======================================================================
    … …  
    3142c                            on peut les mettre a 0 au debut
    3243c ALE-----input-R-energie disponible pour soulevement
     44c ALP-----input-R-puissance disponible pour soulevement
    3345c
    3446C d_h-----output-R-increment de l'enthalpie potentielle (h)
    … …  
    3951c dnwd----output-R-saturated downdraft mass flux (kg/m**2/s)
    4052c dnwd0---output-R-unsaturated downdraft mass flux (kg/m**2/s)
     53c Ma------output-R-adiabatic ascent mass flux (kg/m2/s)
     54c mip-----output-R-mass flux shed by adiabatic ascent (kg/m2/s)
     55c Vprecip-output-R-vertical profile of precipitations (kg/m2/s)
     56c Tconv---output-R-environment temperature seen by convective scheme (K)
    4157c Cape----output-R-CAPE (J/kg)
     58c Cin ----output-R-CIN  (J/kg)
    4259c Tvp-----output-R-Temperature virtuelle d'une parcelle soulevee
    4360c                  adiabatiquement a partir du niveau 1 (K)
    4461c deltapb-output-R-distance entre LCL et base de la colonne (<0 ; Pa)
    4562c Ice_flag-input-L-TRUE->prise en compte de la thermodynamique de la glace
     63c dd_t-----output-R-increment de la temperature du aux descentes precipitantes
     64c dd_q-----output-R-increment de la vapeur d'eau du aux desc precip
    4665c======================================================================
    4766c
    4867#include "dimensions.h"
    49 cym#include "dimphy.h"
     68#include "dimphy.h"
    5069c
    5170      integer NTRAC
    5271      PARAMETER (NTRAC=nqmx-2)
    5372c
    54        INTEGER iflag_con
     73       INTEGER iflag_con,iflag_clos
    5574c
    5675       REAL dtime, paprs(klon,klev+1),pplay(klon,klev)
    5776       REAL t(klon,klev),q(klon,klev),u(klon,klev),v(klon,klev)
     77       REAL t_wake(klon,klev),q_wake(klon,klev)
    5878       REAL tra(klon,klev,ntrac)
    5979       INTEGER ntra
    60        REAL work1(klon,klev),work2(klon,klev)
     80       REAL work1(klon,klev),work2(klon,klev),ptop2(klon)
    6181       REAL pmflxr(klon,klev+1),pmflxs(klon,klev+1)
     82       REAL ALE(klon),ALP(klon)
    6283c
    6384       REAL d_t(klon,klev),d_q(klon,klev),d_u(klon,klev),d_v(klon,klev)
     85       REAL dd_t(klon,klev),dd_q(klon,klev)
    6486       REAL d_tra(klon,klev,ntrac)
    6587       REAL rain(klon),snow(klon)
    … …  
    6890       REAL em_ph(klon,klev+1),em_p(klon,klev)
    6991       REAL upwd(klon,klev),dnwd(klon,klev),dnwdbis(klon,klev)
    70        REAL Ma(klon,klev),cape(klon),tvp(klon,klev)
     92       REAL Ma(klon,klev), mip(klon,klev),Vprecip(klon,klev)
    7193       real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
     94       REAL cape(klon),cin(klon),tvp(klon,klev)
     95       REAL Tconv(klon,klev)
     96c
     97cCR:test: on passe lentr et alim_star des thermiques
     98       INTEGER lalim_conv(klon)
     99       REAL wght_th(klon,klev)
     100       REAL em_sig1feed ! sigma at lower bound of feeding layer
     101       REAL em_sig2feed ! sigma at upper bound of feeding layer
     102       REAL em_wght(klev) ! weight density determining the feeding mixture
     103con enleve le save
     104c       SAVE em_sig1feed,em_sig2feed,em_wght
     105c
    72106       INTEGER iflag(klon)
    73107       REAL rflag(klon)
    … …  
    77111       REAL qcondc(klon,klev)
    78112       REAL wd(klon)
    79 c
     113       REAL Plim1(klon),Plim2(klon),asupmax(klon,klev)
     114       REAL supmax0(klon),asupmaxmin(klon)
     115c
     116       REAL sigd(klon)
    80117       REAL zx_t,zdelta,zx_qs,zcor
    81118c
     119       INTEGER iflag_mix
     120       SAVE iflag_mix
    82121       INTEGER noff, minorig
    83122       INTEGER i,k,itra
    84        REAL qs(klon,klev)
    85 cym       REAL cbmf(klon)
    86 cym       SAVE cbmf
    87        REAL,ALLOCATABLE,SAVE :: cbmf(:)
    88 c$OMP THREADPRIVATE(cbmf)
     123       REAL qs(klon,klev),qs_wake(klon,klev)
     124       REAL cbmf(klon)
     125       SAVE cbmf
     126!       REAL cbmflast(klon)
    89127       INTEGER ifrst
    90128       SAVE ifrst
    … …  
    92130c$OMP THREADPRIVATE(ifrst)
    93131
     132c
     133C     Variables supplementaires liees au bilan d'energie
     134c      Real paire(klon)
     135      Real ql(klon,klev)
     136c      Save paire
     137      Save ql
     138      Real t1(klon,klev),q1(klon,klev)
     139      Save t1,q1
     140c      Data paire /1./
     141c
     142C     Variables liees au bilan d'energie et d'enthalpi
     143      REAL ztsol(klon)
     144      REAL      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
     145     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
     146      SAVE      h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
     147     $        , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
     148      REAL      d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
     149      REAL      d_h_vcol_phy
     150      REAL      fs_bound, fq_bound
     151      SAVE      d_h_vcol_phy
     152      REAL      zero_v(klon)
     153      CHARACTER*15 ztit
     154      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
     155      SAVE      ip_ebil
     156      DATA      ip_ebil/2/
     157      INTEGER   if_ebil ! level for energy conserv. dignostics
     158      SAVE      if_ebil
     159      DATA      if_ebil/2/
     160c+jld ec_conser
     161      REAL d_t_ec(klon,klev)    ! tendance du a la conersion Ec -> E thermique
     162      REAL ZRCPD
     163c-jld ec_conser
     164c
    94165#include "YOMCST.h"
     166#include "YOMCST2.h"
    95167#include "YOETHF.h"
    96168#include "FCTTRE.h"
    97169c
     170
     171c    Copy T into Tconv
     172      DO k = 1,klev
     173        DO i = 1,klon
     174          Tconv(i,k) = T(i,k)
     175        ENDDO
     176      ENDDO
     177c
     178      IF (if_ebil.ge.1) THEN
     179        DO i=1,klon
     180          ztsol(i) = t(i,1)
     181          zero_v(i)=0.
     182          Do k = 1,klev
     183            ql(i,k) = 0.
     184          ENDDO
     185        END DO
     186      END IF
    98187c
    99188cym
    … …  
    102191      IF (ifrst .EQ. 0) THEN
    103192         ifrst = 1
    104          allocate(cbmf(klon))
     193c
     194C===========================================================================
     195C    READ IN PARAMETERS FOR THE CLOSURE AND THE MIXING DISTRIBUTION
     196C===========================================================================
     197C
     198      if (iflag_con.eq.3) then
     199c     CALL cv3_inicp(iflag_clos,iflag_mix)
     200      CALL cv3_inip(iflag_mix)
     201      endif
     202c
     203C===========================================================================
     204C    READ IN PARAMETERS FOR CONVECTIVE INHIBITION BY TROPOS. DRYNESS
     205C===========================================================================
     206C
     207         open (56,file='supcrit.data')
     208         read (56,*) Supcrit1, Supcrit2
     209         close (56)
     210c
     211         print*, 'supcrit1, supcrit2' ,supcrit1, supcrit2
     212C
     213C===========================================================================
     214C      Initialisation pour les bilans d'eau et d'energie
     215C===========================================================================
     216         IF (if_ebil.ge.1) d_h_vcol_phy=0.
     217c
    105218         DO i = 1, klon
    106219          cbmf(i) = 0.
    107220         ENDDO
    108       ENDIF
     221      ENDIF   !(ifrst .EQ. 0)
    109222
    110223      DO k = 1, klev+1
    … …  
    120233      ENDDO
    121234      ENDDO
    122 
    123 c
     235c
     236!
     237!  Feeding layer
     238!
     239      em_sig1feed = 1.
     240      em_sig2feed = 0.97
     241c      em_sig2feed = 0.8
     242! Relative Weight densities
     243       do k=1,klev
     244         em_wght(k)=1.
     245       end do
     246cCRtest: couche alim des tehrmiques ponderee par a*
     247c       DO i = 1, klon
     248c         do k=1,lalim_conv(i)
     249c         em_wght(k)=wght_th(i,k)
     250c         print*,'em_wght=',em_wght(k),wght_th(i,k)
     251c       end do
     252c      END DO
     253
    124254      if (iflag_con .eq. 4) then
    125255      DO k = 1, klev
    … …  
    130260         zcor=1./(1.-retv*zx_qs)
    131261         qs(i,k)=zx_qs*zcor
     262        ENDDO
     263        DO i = 1, klon
     264         zx_t = t_wake(i,k)
     265         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
     266         zx_qs= MIN(0.5 , r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0)
     267         zcor=1./(1.-retv*zx_qs)
     268         qs_wake(i,k)=zx_qs*zcor
    132269        ENDDO
    133270      ENDDO
    … …  
    143280         qs(i,k)=zx_qs
    144281        ENDDO
     282        DO i = 1, klon
     283         zx_t = t_wake(i,k)
     284         zdelta=MAX(0.,SIGN(1.,rtt-zx_t))
     285         zx_qs= r2es * FOEEW(zx_t,zdelta)/em_p(i,k)/100.0
     286         zx_qs= MIN(0.5,zx_qs)
     287         zcor=1./(1.-retv*zx_qs)
     288         zx_qs=zx_qs*zcor
     289         qs_wake(i,k)=zx_qs
     290        ENDDO
    145291      ENDDO
    146292      endif ! iflag_con
    … …  
    149295
    150296C Main driver for convection:
    151 C               iflag_con = 3  -> equivalent to convect3
     297C               iflag_con=3 -> nvlle version de KE (JYG)
     298C               iflag_con = 30  -> equivalent to convect3
    152299C               iflag_con = 4  -> equivalent to convect1/2
     300c
     301c
     302      if (iflag_con.eq.30) then
    153303
    154304      CALL cv_driver(klon,klev,klev+1,ntra,iflag_con,
    … …  
    160310     $              dtime,Ma,upwd,dnwd,dnwdbis,qcondc,wd,cape,
    161311     $              da,phi,mp)
    162 
     312     
     313      else
     314
     315      CALL cva_driver(klon,klev,klev+1,ntra,
     316     $              iflag_con,iflag_mix,iflag_clos,dtime,
     317     :              t,q,qs,t_wake,q_wake,qs_wake,u,v,tra,
     318     $              em_p,em_ph,
     319     .              ALE,ALP,
     320     .              em_sig1feed,em_sig2feed,em_wght,
     321     .              iflag,d_t,d_q,d_u,d_v,d_tra,rain,kbas,ktop,
     322     $              cbmf,work1,work2,ptop2,sigd,
     323     $              Ma,mip,Vprecip,upwd,dnwd,dnwdbis,qcondc,wd,
     324     $              cape,cin,tvp,
     325     $              dd_t,dd_q,Plim1,Plim2,asupmax,supmax0,
     326     $              asupmaxmin,lalim_conv)
     327      endif 
    163328C------------------------------------------------------------------
    164329
    … …  
    176341        ENDDO
    177342      ENDDO
     343c
     344       if (iflag_con.eq.30) then
    178345       DO itra = 1,ntra
    179346        DO k = 1, klev
    … …  
    182349         ENDDO
    183350        ENDDO
    184        ENDDO
     351       ENDDO
     352       endif
     353
     354      DO k = 1, klev
     355        DO i = 1, klon
     356          t1(i,k) = t(i,k)+ d_t(i,k)
     357          q1(i,k) = q(i,k)+ d_q(i,k)
     358        ENDDO
     359      ENDDO
     360c
     361cc      IF (if_ebil.ge.2) THEN
     362cc        ztit='after convect'
     363cc        CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
     364cc     e      , t1,q1,ql,qs,u,v,paprs,pplay
     365cc     s      , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
     366cc         call diagphy(paire,ztit,ip_ebil
     367cc     e      , zero_v, zero_v, zero_v, zero_v, zero_v
     368cc     e      , zero_v, rain, zero_v, ztsol
     369cc     e      , d_h_vcol, d_qt, d_ec
     370cc     s      , fs_bound, fq_bound )
     371cc      END IF
     372C
     373c
    185374c les traceurs ne sont pas mis dans cette version de convect4:
    186375      if (iflag_con.eq.4) then
    … …  
    193382       ENDDO
    194383      endif
    195  
     384      print*, 'concvl->: dd_t,dd_q ',dd_t(1,1),dd_q(1,1)
     385
    196386      RETURN
    197387      END

  • LMDZ4/trunk/libf/phylmd/conema3.h

    r766 r879  
    11!
    22! $Header$
     3!-- Modified by : Filiberti M-A 06/2005
    34!
    45      real epmax             ! 0.993
    56      logical ok_adj_ema      ! F
    67      integer iflag_clw      ! 0
     8          integer iflag_cvl_sigd
     9      real sig1feed      ! 1.
     10      real sig2feed      ! 0.95
    711
    8       common/comconema/epmax,ok_adj_ema,iflag_clw
    9 !$OMP THREADPRIVATE(/comconema/)
     12      common/comconema1/epmax,ok_adj_ema,iflag_clw,sig1feed,sig2feed
     13      common/comconema2/iflag_cvl_sigd
     14
     15!      common/comconema/epmax,ok_adj_ema,iflag_clw
     16!$OMP THREADPRIVATE(/comconema1/)
     17!$OMP THREADPRIVATE(/comconema2/)

  • LMDZ4/trunk/libf/phylmd/conf_phys.F90

    r878 r879  
    77  subroutine conf_phys(ocean, ok_veget, ok_journe, ok_mensuel, ok_instan, ok_hf, &
    88 &                     seuil_inversion, &
    9  &                     fact_cldcon, facttemps,ok_newmicro,iflag_cldcon, &
     9 &                     fact_cldcon, facttemps,ok_newmicro,iflag_radia,&
     10 &                     iflag_cldcon, &
    1011 &                     iflag_ratqs,ratqsbas,ratqshaut, &
    1112 &                     ok_ade, ok_aie, &
    1213 &                     bl95_b0, bl95_b1,&
    13  &                     iflag_thermals,nsplit_thermals)
     14 &                     iflag_thermals,nsplit_thermals, &
     15 &                     iflag_coupl,iflag_clos,iflag_wake )
    1416
    1517   use IOIPSL
    … …  
    4648  character (len = 6)  :: ocean
    4749  logical              :: ok_veget, ok_newmicro
     50  integer              :: iflag_radia
    4851  logical              :: ok_journe, ok_mensuel, ok_instan, ok_hf
    4952  LOGICAL              :: ok_ade, ok_aie
    … …  
    6164  real,SAVE           :: fact_cldcon_omp, facttemps_omp,ratqsbas_omp
    6265  real,SAVE           :: ratqshaut_omp
     66  integer,SAVE        :: iflag_radia_omp
    6367  integer,SAVE        :: iflag_cldcon_omp, ip_ebil_phy_omp
    6468  integer,SAVE        :: iflag_ratqs_omp
    … …  
    7377  integer :: iflag_thermals,nsplit_thermals
    7478  integer,SAVE :: iflag_thermals_omp,nsplit_thermals_omp
     79  integer :: iflag_coupl
     80  integer :: iflag_clos
     81  integer :: iflag_wake
     82  integer,SAVE :: iflag_coupl_omp,iflag_clos_omp,iflag_wake_omp
     83  integer,SAVE :: iflag_cvl_sigd_omp
    7584
    7685  REAL,SAVE :: R_ecc_omp,R_peri_omp,R_incl_omp,solaire_omp,co2_ppm_omp
    … …  
    461470
    462471!
     472!Config Key  = iflag_radia
     473!Config Desc = 
     474!Config Def  = 1
     475!Config Help =
     476!
     477  iflag_radia_omp = 1
     478  call getin('iflag_radia',iflag_radia_omp)
     479
    463480!Config Key  = iflag_cldcon
    464481!Config Desc = 
    … …  
    658675  call getin('nsplit_thermals',nsplit_thermals_omp)
    659676
    660 
     677!
     678!Config Key  = iflag_coupl
     679!Config Desc =
     680!Config Def  = 0
     681!Config Help =
     682!
     683  iflag_coupl = 0
     684  call getin('iflag_coupl',iflag_coupl_omp)
     685
     686!
     687!Config Key  = iflag_clos
     688!Config Desc = 
     689!Config Def  = 0
     690!Config Help =
     691!
     692  iflag_clos = 1
     693  call getin('iflag_clos',iflag_clos_omp)
     694!
     695!Config Key  = iflag_cvl_sigd
     696!Config Desc = 
     697!Config Def  = 0
     698!Config Help =
     699!
     700  iflag_cvl_sigd = 0
     701  call getin('iflag_cvl_sigd',iflag_cvl_sigd_omp)
     702
     703!Config Key  = iflag_wake
     704!Config Desc = 
     705!Config Def  = 0
     706!Config Help =
     707!
     708  iflag_wake = 0
     709  call getin('iflag_wake',iflag_wake_omp)
    661710
    662711!
    … …  
    856905    ratqsbas = ratqsbas_omp
    857906    ratqshaut = ratqshaut_omp
     907    iflag_radia = iflag_radia_omp
    858908    iflag_cldcon = iflag_cldcon_omp
    859909    iflag_ratqs = iflag_ratqs_omp
    … …  
    861911    iflag_thermals = iflag_thermals_omp
    862912    nsplit_thermals = nsplit_thermals_omp
     913    iflag_coupl = iflag_coupl_omp
     914    iflag_clos = iflag_clos_omp
     915    iflag_wake = iflag_wake_omp
     916    iflag_cvl_sigd = iflag_cvl_sigd_omp
    863917    type_run = type_run_omp
    864918    ok_isccp = ok_isccp_omp
    … …  
    915969  write(numout,*)' iflag_pdf = ', iflag_pdf
    916970  write(numout,*)' iflag_cldcon = ', iflag_cldcon
     971  write(numout,*)' iflag_radia = ', iflag_radia
    917972  write(numout,*)' iflag_ratqs = ', iflag_ratqs
    918973  write(numout,*)' seuil_inversion = ', seuil_inversion

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

    r829 r879  
    3030C   ***                     IT MUST BE LESS THAN 0              ***
    3131
    32 #include "cvparam3.h"
     32#include "cv3param.h"
    3333#include "conema3.h"
    3434
    … …  
    4848
    4949c -- "microphysical" parameters:
    50 
    51       sigd   = 0.01
     50       sigdz=0.01
     51c      sigd=0.003
     52c     sigd   = 0.01
     53cCR:test sur la fraction des descentes precipitantes
    5254      spfac  = 0.15
    5355      pbcrit = 150.0
    5456      ptcrit = 500.0
    55 cIM cf. FH     epmax  = 0.993
     57      epmax  = 0.993
    5658
    5759      omtrain = 45.0 ! used also for snow (no disctinction rain/snow)
    … …  
    6163      dtovsh = -0.2 ! dT for overshoot
    6264      dpbase = -40. ! definition cloud base (400m above LCL)
    63       dttrig = 5.   ! (loose) condition for triggering
     65ccc      dttrig = 5.   ! (loose) condition for triggering
     66      dttrig = 10.   ! (loose) condition for triggering
    6467
    6568c -- rate of approach to quasi-equilibrium:
    6669
    6770      dtcrit = -2.0
    68       tau    = 8000.
     71c      dtcrit = -5.0
     72c      tau    = 3000.
     73cc      tau = 1800.
     74      tau=8000.
    6975      beta   = 1.0 - delt/tau
    70       alpha  = 1.5E-3 * delt/tau
     76      alpha1 = 1.5e-3
     77      alpha  = alpha1 * delt/tau
    7178c increase alpha to compensate W decrease:
    7279      alpha  = alpha*1.5
    … …  
    109116
    110117#include "cvthermo.h"
    111 #include "cvparam3.h"
     118#include "cv3param.h"
    112119
    113120
    … …  
    138145        gz(i,k)=gz(i,k-1)+0.5*rrd*(tvx+tvy)     !convect3
    139146     &          *(p(i,k-1)-p(i,k))/ph(i,k)      !convect3
    140 
     147c
     148cc        print *,' gz(',k,')',gz(i,k),' tvx',tvx,' tvy ',tvy
     149c
    141150c ori         gz(i,k)=gz(i,k-1)+hrd*(tv(i,k-1)+tv(i,k))
    142151c ori    &         *(p(i,k-1)-p(i,k))/ph(i,k)
    … …  
    158167      end
    159168
    160       SUBROUTINE cv3_feed(len,nd,t,q,qs,p,ph,hm,gz
    161      :                  ,nk,icb,icbmax,iflag,tnk,qnk,gznk,plcl)
     169      SUBROUTINE cv3_feed(len,nd,t,q,u,v,p,ph,hm,gz
     170     :                  ,p1feed,p2feed,wght
     171     :                  ,wghti,tnk,thnk,qnk,qsnk,unk,vnk
     172     :                  ,cpnk,hnk,nk,icb,icbmax,iflag,gznk,plcl)
    162173      implicit none
    163174
    … …  
    176187C================================================================
    177188
    178 #include "cvparam3.h"
     189#include "cv3param.h"
     190#include "cvthermo.h"
    179191
    180192c inputs:
    181193          integer len, nd
    182       real t(len,nd), q(len,nd), qs(len,nd), p(len,nd)
     194      real t(len,nd), q(len,nd), p(len,nd)
     195      real u(len,nd), v(len,nd)
    183196      real hm(len,nd), gz(len,nd)
    184197      real ph(len,nd+1)
    185 
     198      real p1feed(len)
     199c,  wght(len)
     200      real wght(nd)
     201c input-output
     202      real p2feed(len)
    186203c outputs:
    187204          integer iflag(len), nk(len), icb(len), icbmax
    188       real tnk(len), qnk(len), gznk(len), plcl(len)
     205c      real   wghti(len)
     206      real wghti(len,nd)
     207      real   tnk(len), thnk(len), qnk(len), qsnk(len)
     208      real   unk(len), vnk(len)
     209      real   cpnk(len), hnk(len), gznk(len)
     210      real   plcl(len)
    189211
    190212c local variables:
    191       integer i, k
     213      integer i, k, iter, niter
    192214      integer ihmin(len)
    193215      real work(len)
    194       real pnk(len), qsnk(len), rh(len), chi(len)
    195       real A, B ! convect3
    196 cym
    197       plcl=0.0
    198 c@ !-------------------------------------------------------------------
    199 c@ ! --- Find level of minimum moist static energy
    200 c@ ! --- If level of minimum moist static energy coincides with
    201 c@ ! --- or is lower than minimum allowable parcel origin level,
    202 c@ ! --- set iflag to 6.
    203 c@ !-------------------------------------------------------------------
    204 c@
    205 c@       do 180 i=1,len
    206 c@        work(i)=1.0e12
    207 c@        ihmin(i)=nl
    208 c@  180  continue
    209 c@       do 200 k=2,nlp
    210 c@         do 190 i=1,len
    211 c@          if((hm(i,k).lt.work(i)).and.
    212 c@      &      (hm(i,k).lt.hm(i,k-1)))then
    213 c@            work(i)=hm(i,k)
    214 c@            ihmin(i)=k
    215 c@          endif
    216 c@  190    continue
    217 c@  200  continue
    218 c@       do 210 i=1,len
    219 c@         ihmin(i)=min(ihmin(i),nlm)
    220 c@         if(ihmin(i).le.minorig)then
    221 c@           iflag(i)=6
    222 c@         endif
    223 c@  210  continue
    224 c@ c
    225 c@ !-------------------------------------------------------------------
    226 c@ ! --- Find that model level below the level of minimum moist static
    227 c@ ! --- energy that has the maximum value of moist static energy
    228 c@ !-------------------------------------------------------------------
    229 c@ 
    230 c@       do 220 i=1,len
    231 c@        work(i)=hm(i,minorig)
    232 c@        nk(i)=minorig
    233 c@  220  continue
    234 c@       do 240 k=minorig+1,nl
    235 c@         do 230 i=1,len
    236 c@          if((hm(i,k).gt.work(i)).and.(k.le.ihmin(i)))then
    237 c@            work(i)=hm(i,k)
    238 c@            nk(i)=k
    239 c@          endif
    240 c@  230     continue
    241 c@  240  continue
    242 
     216      real pup(len),plo(len),pfeed(len)
     217      real plclup(len),plcllo(len),plclfeed(len)
     218      real posit(len)
     219      logical nocond(len)
     220!
    243221!-------------------------------------------------------------------
    244222! --- Origin level of ascending parcels for convect3:
    … …  
    247225         do 220 i=1,len
    248226          nk(i)=minorig
     227          gznk(i)=gz(i,nk(i))
    249228  220    continue
    250 
     229!
     230!-------------------------------------------------------------------
     231! --- Adjust feeding layer thickness so that lifting up to the top of
     232! --- the feeding layer does not induce condensation (i.e. so that
     233! --- plcl < p2feed).
     234! --- Method : iterative secant method.
     235!-------------------------------------------------------------------
     236!
     237c 1- First bracketing of the solution : ph(nk+1), p2feed
     238c
     239c 1.a- LCL associated to p2feed
     240      do i = 1,len
     241        pup(i) = p2feed(i)
     242      enddo
     243         call cv3_vertmix(len,nd,iflag,p1feed,pup,p,ph
     244     i              ,t,q,u,v,wght
     245     o              ,wghti,nk,tnk,thnk,qnk,qsnk,unk,vnk,plclup)
     246c 1.b- LCL associated to ph(nk+1)
     247      do i = 1,len
     248        plo(i) = ph(i,nk(i)+1)
     249      enddo
     250         call cv3_vertmix(len,nd,iflag,p1feed,plo,p,ph
     251     i              ,t,q,u,v,wght
     252     o              ,wghti,nk,tnk,thnk,qnk,qsnk,unk,vnk,plcllo)
     253c 2- Iterations
     254      niter = 5
     255      do iter = 1,niter
     256        do i = 1,len
     257          plcllo(i) = min(plo(i),plcllo(i))
     258          plclup(i) = max(pup(i),plclup(i))
     259          nocond(i) = plclup(i).le.pup(i)
     260        enddo
     261        do i = 1,len
     262          if(nocond(i)) then
     263             pfeed(i)=pup(i)
     264          else
     265             pfeed(i) = (pup(i)*(plo(i)-plcllo(i))+
     266     :                plo(i)*(plclup(i)-pup(i)))/
     267     :            (plo(i)-plcllo(i)+plclup(i)-pup(i))
     268          endif
     269        enddo
     270         call cv3_vertmix(len,nd,iflag,p1feed,pfeed,p,ph
     271     i              ,t,q,u,v,wght
     272     o              ,wghti,nk,tnk,thnk,qnk,qsnk,unk,vnk,plclfeed)
     273        do i = 1,len
     274          posit(i) = (sign(1.,plclfeed(i)-pfeed(i))+1.)*0.5
     275          if (plclfeed(i) .eq. pfeed(i)) posit(i) = 1.
     276c- posit = 1 when lcl is below top of feeding layer (plclfeed>pfeed)
     277c-               => pup=pfeed
     278c- posit = 0 when lcl is above top of feeding layer (plclfeed<pfeed)
     279c-               => plo=pfeed
     280          pup(i) = posit(i)*pfeed(i) + (1.-posit(i))*pup(i)
     281          plo(i) = (1.-posit(i))*pfeed(i) + posit(i)*plo(i)
     282          plclup(i) = posit(i)*plclfeed(i) + (1.-posit(i))*plclup(i)
     283          plcllo(i) = (1.-posit(i))*plclfeed(i) + posit(i)*plcllo(i)
     284        enddo
     285      enddo       !  iter
     286      do i = 1,len
     287        p2feed(i) = pfeed(i)
     288        plcl(i) = plclfeed(i)
     289      enddo
     290!
     291      do 175 i=1,len
     292         cpnk(i)=cpd*(1.0-qnk(i))+cpv*qnk(i)
     293         hnk(i)=gz(i,1)+cpnk(i)*tnk(i)
     294 175  continue
     295!
    251296!-------------------------------------------------------------------
    252297! --- Check whether parcel level temperature and specific humidity
    … …  
    254299!-------------------------------------------------------------------
    255300       do 250 i=1,len
    256        if( (     ( t(i,nk(i)).lt.250.0    )
    257      &       .or.( q(i,nk(i)).le.0.0      )     )
    258 c@      &       .or.( p(i,ihmin(i)).lt.400.0 )  )
     301       if( (     ( tnk(i).lt.250.0    )
     302     &       .or.( qnk(i).le.0.0      ) )
    259303     &   .and.
    260304     &       ( iflag(i).eq.0) ) iflag(i)=7
    261305 250   continue
    262 !-------------------------------------------------------------------
    263 ! --- Calculate lifted condensation level of air at parcel origin level
    264 ! --- (Within 0.2% of formula of Bolton, MON. WEA. REV.,1980)
    265 !-------------------------------------------------------------------
    266 
    267        A = 1669.0 ! convect3
    268        B = 122.0  ! convect3
    269 
    270        do 260 i=1,len
    271 
    272         if (iflag(i).ne.7) then ! modif sb Jun7th 2002
    273 
    274         tnk(i)=t(i,nk(i))
    275         qnk(i)=q(i,nk(i))
    276         gznk(i)=gz(i,nk(i))
    277         pnk(i)=p(i,nk(i))
    278         qsnk(i)=qs(i,nk(i))
    279 c
    280         rh(i)=qnk(i)/qsnk(i)
    281 c ori        rh(i)=min(1.0,rh(i)) ! removed for convect3
    282 c ori        chi(i)=tnk(i)/(1669.0-122.0*rh(i)-tnk(i))
    283         chi(i)=tnk(i)/(A-B*rh(i)-tnk(i)) ! convect3
    284         plcl(i)=pnk(i)*(rh(i)**chi(i))
    285         if(((plcl(i).lt.200.0).or.(plcl(i).ge.2000.0))
    286      &   .and.(iflag(i).eq.0))iflag(i)=8
    287  
    288         endif ! iflag=7 
    289 
    290  260   continue
    291 
     306c
    292307!-------------------------------------------------------------------
    293308! --- Calculate first level above lcl (=icb)
    … …  
    322337 290  continue
    323338c
     339
     340c     print*,'icb dans cv3_feed '
     341c     write(*,'(64i2)') icb(2:len-1)
     342c     call dump2d(64,43,'plcl dans cv3_feed ',plcl(2:len-1))
     343
    324344      do 300 i=1,len
    325345c@        if((icb(i).ge.nlm).and.(iflag(i).eq.0))iflag(i)=9
    … …  
    342362      end
    343363
    344       SUBROUTINE cv3_undilute1(len,nd,t,q,qs,gz,plcl,p,nk,icb
     364      SUBROUTINE cv3_undilute1(len,nd,t,qs,gz,plcl,p,icb,tnk,qnk,gznk
    345365     :                       ,tp,tvp,clw,icbs)
    346366      implicit none
    … …  
    360380
    361381#include "cvthermo.h"
    362 #include "cvparam3.h"
     382#include "cv3param.h"
    363383
    364384c inputs:
    365385      integer len, nd
    366       integer nk(len), icb(len)
    367       real t(len,nd), q(len,nd), qs(len,nd), gz(len,nd)
    368       real p(len,nd)
     386      integer icb(len)
     387      real t(len,nd), qs(len,nd), gz(len,nd)
     388      real tnk(len), qnk(len), gznk(len)
     389      real p(len,nd)
    369390      real plcl(len) ! convect3
    370391
    … …  
    377398      real tg, qg, alv, s, ahg, tc, denom, es, rg
    378399      real ah0(len), cpp(len)
    379       real tnk(len), qnk(len), gznk(len), ticb(len), gzicb(len)
     400      real ticb(len), gzicb(len)
    380401      real qsicb(len) ! convect3
    381402      real cpinv(len) ! convect3
    … …  
    388409!-------------------------------------------------------------------
    389410
    390       do 320 i=1,len
    391         tnk(i)=t(i,nk(i))
    392         qnk(i)=q(i,nk(i))
    393         gznk(i)=gz(i,nk(i))
    394 c ori        ticb(i)=t(i,icb(i))
    395 c ori        gzicb(i)=gz(i,icb(i))
    396  320  continue
    397411c
    398412c   ***  Calculate certain parcel quantities, including static energy   ***
    … …  
    547561c
    548562
    549         do i=1,len             
    550          ticb(i)=t(i,icb(i)+1)   
    551          gzicb(i)=gz(i,icb(i)+1) 
    552          qsicb(i)=qs(i,icb(i)+1) 
    553         enddo                   
     563        do i=1,len
     564         ticb(i)=t(i,icb(i)+1)
     565         gzicb(i)=gz(i,icb(i)+1)
     566         qsicb(i)=qs(i,icb(i)+1)
     567        enddo
    554568
    555569        do 460 i=1,len
    … …  
    627641      end
    628642
    629       SUBROUTINE cv3_trigger(len,nd,icb,plcl,p,th,tv,tvp
    630      o                ,pbase,buoybase,iflag,sig,w0)
     643      SUBROUTINE cv3_trigger(len,nd,icb,plcl,p,th,tv,tvp,thnk,
     644     o                pbase,buoybase,iflag,sig,w0)
    631645      implicit none
    632646
    … …  
    646660!-------------------------------------------------------------------
    647661
    648 #include "cvparam3.h"
     662#include "cv3param.h"
    649663
    650664c input:
    … …  
    653667      real plcl(len), p(len,nd)
    654668      real th(len,nd), tv(len,nd), tvp(len,nd)
     669      real thnk(len)
    655670
    656671c output:
    … …  
    714729
    715730       tdif = buoybase(i)
    716        ath1 = th(i,1)
     731       ath1 = thnk(i)
    717732       ath  = th(i,icb(i)-1) - dttrig
    718733
    … …  
    738753     :    ,t1,q1,qs1,u1,v1,gz1,th1
    739754     :    ,tra1
    740      :    ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1 
     755     :    ,h1,lv1,cpn1,p1,ph1,tv1,tp1,tvp1,clw1
    741756     :    ,sig1,w01
    742757     o    ,iflag,nk,icb,icbs
    … …  
    748763      implicit none
    749764
    750 #include "cvparam3.h"
     765#include "cv3param.h"
    751766
    752767c inputs:
    … …  
    773788      real tvp(nloc,nd),clw(nloc,nd)
    774789      real th(nloc,nd)
    775       real sig(nloc,nd), w0(nloc,nd) 
     790      real sig(nloc,nd), w0(nloc,nd)
    776791      real tra(nloc,nd,ntra)
    777792
    … …  
    807822 110  continue
    808823
    809 c      do 121 j=1,ntra
    810 c      do 111 k=1,nd
    811 c       nn=0
    812 c      do 101 i=1,len
    813 c      if(iflag1(i).eq.0)then
    814 c       nn=nn+1
    815 c       tra(nn,k,j)=tra1(i,k,j)
    816 c      endif
    817 c 101  continue
    818 c 111  continue
    819 c 121  continue
     824      do 121 j=1,ntra
     825ccccc      do 111 k=1,nl+1
     826      do 111 k=1,nd
     827       nn=0
     828      do 101 i=1,len
     829      if(iflag1(i).eq.0)then
     830       nn=nn+1
     831       tra(nn,k,j)=tra1(i,k,j)
     832      endif
     833 101  continue
     834 111  continue
     835 121  continue
    820836
    821837      if (nn.ne.ncum) then
    … …  
    845861
    846862      SUBROUTINE cv3_undilute2(nloc,ncum,nd,icb,icbs,nk
    847      :                       ,tnk,qnk,gznk,t,q,qs,gz
     863     :                       ,tnk,qnk,gznk,hnk,t,q,qs,gz
    848864     :                       ,p,h,tv,lv,pbase,buoybase,plcl
    849865     o                       ,inb,tp,tvp,clw,hp,ep,sigp,buoy)
    … …  
    854870C     FIND THE REST OF THE LIFTED PARCEL TEMPERATURES
    855871C     &
    856 C     COMPUTE THE PRECIPITATION EFFICIENCIES AND THE 
     872C     COMPUTE THE PRECIPITATION EFFICIENCIES AND THE
    857873C     FRACTION OF PRECIPITATION FALLING OUTSIDE OF CLOUD
    858874C     &
    … …  
    869885
    870886#include "cvthermo.h"
    871 #include "cvparam3.h"
     887#include "cv3param.h"
    872888#include "conema3.h"
    873889
    … …  
    878894      real p(nloc,nd)
    879895      real tnk(nloc), qnk(nloc), gznk(nloc)
     896      real hnk(nloc)
    880897      real lv(nloc,nd), tv(nloc,nd), h(nloc,nd)
    881898      real pbase(nloc), buoybase(nloc), plcl(nloc)
    … …  
    893910      real ah0(nloc), cape(nloc), capem(nloc), byp(nloc)
    894911      logical lcape(nloc)
     912      integer iposit(nloc)
    895913
    896914!=====================================================================
    … …  
    10541072c=====================================================================
    10551073c  --- EFFECTIVE VERTICAL PROFILE OF BUOYANCY (convect3 only):
    1056 c===================================================================== 
     1074c=====================================================================
    10571075
    10581076c-- this is for convect3 only:
    … …  
    10621080      do 500 i=1,ncum
    10631081       do 501 k=1,nl
    1064         buoy(i,k)=tvp(i,k)-tv(i,k) 
     1082        buoy(i,k)=tvp(i,k)-tv(i,k)
    10651083 501   continue
    10661084 500  continue
    … …  
    10751093        endif
    10761094 506   continue
    1077 cIM cf. CRio/JYG 270807   buoy(icb(i),k)=buoybase(i)
    1078        buoy(i,icb(i))=buoybase(i)
     1095c       buoy(icb(i),k)=buoybase(i)
     1096      buoy(i,icb(i))=buoybase(i)
    10791097 505  continue
    10801098
    … …  
    10901108      do 510 i=1,ncum
    10911109       inb(i)=nl-1
     1110       iposit(i) = nl
    10921111 510  continue
     1112
     1113c
     1114c--    iposit(i) = first level, above icb, with positive buoyancy
     1115      do k = 1,nl-1
     1116       do i = 1,ncum
     1117        if (k .ge. icb(i) .and. buoy(i,k) .gt. 0.) then
     1118          iposit(i) = min(iposit(i),k)
     1119        endif
     1120       enddo
     1121      enddo
     1122
     1123      do i = 1,ncum
     1124       if (iposit(i) .eq. nl) then
     1125         iposit(i) = icb(i)
     1126       endif
     1127      enddo
    10931128
    10941129      do 530 i=1,ncum
    10951130       do 535 k=1,nl-1
    1096         if ((k.ge.icb(i)).and.(buoy(i,k).lt.dtovsh)) then
     1131        if ((k.ge.iposit(i)).and.(buoy(i,k).lt.dtovsh)) then
    10971132         inb(i)=MIN(inb(i),k)
    10981133        endif
    … …  
    11961231c=====================================================================
    11971232c
    1198 cym      do i=1,ncum*nlp
    1199 cym       hp(i,1)=h(i,1)
    1200 cym      enddo
    1201 
    1202       do k=1,nlp
    1203         do i=1,ncum
    1204           hp(i,k)=h(i,k)
    1205         enddo
     1233      do k = 1,nd
     1234      do i=1,ncum
     1235         hp(i,k)=h(i,k)
     1236      enddo
    12061237      enddo
    12071238
    … …  
    12091240        do 590 i=1,ncum
    12101241        if((k.ge.icb(i)).and.(k.le.inb(i)))then
    1211           hp(i,k)=h(i,nk(i))+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
     1242          hp(i,k)=hnk(i)+(lv(i,k)+(cpd-cpv)*t(i,k))*ep(i,k)*clw(i,k)
    12121243        endif
    12131244 590    continue
    … …  
    12191250      SUBROUTINE cv3_closure(nloc,ncum,nd,icb,inb
    12201251     :                      ,pbase,p,ph,tv,buoy
    1221      o                      ,sig,w0,cape,m)
     1252     o                      ,sig,w0,cape,m,iflag)
    12221253      implicit none
    12231254
    … …  
    12291260
    12301261#include "cvthermo.h"
    1231 #include "cvparam3.h"
     1262#include "cv3param.h"
    12321263
    12331264c input:
    … …  
    12401271c input/output:
    12411272      real sig(nloc,nd), w0(nloc,nd)
     1273      integer iflag(nloc)
    12421274
    12431275c output:
    … …  
    12491281      real deltap, fac, w, amu
    12501282      real dtmin(nloc,nd), sigold(nloc,nd)
     1283      real cbmflast(nloc)
    12511284
    12521285
    … …  
    12621295
    12631296c -------------------------------------------------------
    1264 c -- Reset sig(i) and w0(i) for i>inb and i<icb   
     1297c -- Reset sig(i) and w0(i) for i>inb and i<icb
    12651298c -------------------------------------------------------
    1266      
     1299
    12671300c update sig and w0 above LNB:
    12681301
    … …  
    13161349c -- and after 10 time steps of no convection
    13171350c -------------------------------------------------------------
    1318      
     1351
    13191352      do 400 k=1,nl-1
    13201353       do 410 i=1,ncum
    … …  
    13271360
    13281361c -------------------------------------------------------------
    1329 c -- Calculate convective available potential energy (cape), 
    1330 c -- vertical velocity (w), fractional area covered by   
    1331 c -- undilute updraft (sig), and updraft mass flux (m) 
     1362c -- Calculate convective available potential energy (cape),
     1363c -- vertical velocity (w), fractional area covered by
     1364c -- undilute updraft (sig), and updraft mass flux (m)
    13321365c -------------------------------------------------------------
    13331366
    … …  
    13401373      do i=1,ncum
    13411374       do k=1,nl
    1342          dtmin(i,k)=100.0 
     1375         dtmin(i,k)=100.0
    13431376       enddo
    13441377      enddo
    … …  
    13931426       sig(i,icb(i)-1)=sig(i,icb(i))
    13941427 700  continue
    1395 
    1396 
     1428c
     1429cccc 3. Compute final cloud base mass flux and set iflag to 3 if
     1430cccc    cloud base mass flux is exceedingly small and is decreasing (i.e. if
     1431cccc    the final mass flux (cbmflast) is greater than the target mass flux
     1432cccc    (cbmf) ??).
     1433ccc
     1434cc      do i = 1,ncum
     1435cc       cbmflast(i) = 0.
     1436cc      enddo
     1437ccc
     1438cc      do k= 1,nl
     1439cc       do i = 1,ncum
     1440cc        IF (k .ge. icb(i) .and. k .le. inb(i)) THEN
     1441cc         cbmflast(i) = cbmflast(i)+M(i,k)
     1442cc        ENDIF
     1443cc       enddo
     1444cc      enddo
     1445ccc
     1446cc      do i = 1,ncum
     1447cc       IF (cbmflast(i) .lt. 1.e-6) THEN
     1448cc         iflag(i) = 3
     1449cc       ENDIF
     1450cc      enddo
     1451ccc
     1452cc      do k= 1,nl
     1453cc       do i = 1,ncum
     1454cc        IF (iflag(i) .ge. 3) THEN
     1455cc         M(i,k) = 0.
     1456cc         sig(i,k) = 0.
     1457cc         w0(i,k) = 0.
     1458cc        ENDIF
     1459cc       enddo
     1460cc      enddo
     1461ccc
    13971462c!      cape=0.0
    13981463c!      do 98 i=icb+1,inb
    … …  
    14281493      SUBROUTINE cv3_mixing(nloc,ncum,nd,na,ntra,icb,nk,inb
    14291494     :                    ,ph,t,rr,rs,u,v,tra,h,lv,qnk
    1430      :                    ,hp,tv,tvp,ep,clw,m,sig
     1495     :                    ,unk,vnk,hp,tv,tvp,ep,clw,m,sig
    14311496     :   ,ment,qent,uent,vent,sij,elij,ments,qents,traent)
    14321497      implicit none
    … …  
    14341499!---------------------------------------------------------------------
    14351500! a faire:
    1436 !       - changer rr(il,1) -> qnk(il)
    14371501!   - vectorisation de la partie normalisation des flux (do 789...)
    14381502!---------------------------------------------------------------------
    14391503
    14401504#include "cvthermo.h"
    1441 #include "cvparam3.h"
     1505#include "cv3param.h"
    14421506
    14431507c inputs:
    … …  
    14451509      integer icb(nloc), inb(nloc), nk(nloc)
    14461510      real sig(nloc,nd)
    1447       real qnk(nloc)
     1511      real qnk(nloc),unk(nloc),vnk(nloc)
    14481512      real ph(nloc,nd+1)
    14491513      real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
    … …  
    14581522      real uent(nloc,na,na), vent(nloc,na,na)
    14591523      real sij(nloc,na,na), elij(nloc,na,na)
    1460       real traent(nloc,nd,nd,ntra) 
     1524      real traent(nloc,nd,nd,ntra)
    14611525      real ments(nloc,nd,nd), qents(nloc,nd,nd)
    14621526      real sigij(nloc,nd,nd)
    … …  
    15061570      sij(1:ncum,1:nd,1:nd)=0.0
    15071571     
    1508 c      do k=1,ntra
    1509 c       do j=1,nd  ! instead nlp
    1510 c        do i=1,nd ! instead nlp
    1511 c         do il=1,ncum
    1512 c            traent(il,i,j,k)=tra(il,j,k)
    1513 c         enddo
    1514 c        enddo
    1515 c       enddo
    1516 c      enddo
     1572      do k=1,ntra
     1573       do j=1,nd  ! instead nlp
     1574        do i=1,nd ! instead nlp
     1575         do il=1,ncum
     1576            traent(il,i,j,k)=tra(il,j,k)
     1577         enddo
     1578        enddo
     1579       enddo
     1580      enddo
    15171581      zm(:,:)=0.
    15181582
    … …  
    15301594     :      (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
    15311595
    1532           rti=rr(il,1)-ep(il,i)*clw(il,i)
     1596          rti=qnk(il)-ep(il,i)*clw(il,i)
    15331597          bf2=1.+lv(il,j)*lv(il,j)*rs(il,j)/(rrv*t(il,j)*t(il,j)*cpd)
    15341598          anum=h(il,j)-hp(il,i)+(cpv-cpd)*t(il,j)*(rti-rr(il,j))
    … …  
    15531617         if(sij(il,i,j).gt.0.0.and.sij(il,i,j).lt.0.95)then
    15541618          qent(il,i,j)=sij(il,i,j)*rr(il,i)+(1.-sij(il,i,j))*rti
    1555           uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*u(il,nk(il))
    1556           vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*v(il,nk(il))
     1619          uent(il,i,j)=sij(il,i,j)*u(il,i)+(1.-sij(il,i,j))*unk(il)
     1620          vent(il,i,j)=sij(il,i,j)*v(il,i)+(1.-sij(il,i,j))*vnk(il)
    15571621c!!!      do k=1,ntra
    15581622c!!!      traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
    … …  
    15701634 710  continue
    15711635
    1572 c       do k=1,ntra
    1573 c        do j=minorig,nl
    1574 c         do il=1,ncum
    1575 c          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
    1576 c     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
    1577 c            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
    1578 c     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
    1579 c          endif
    1580 c         enddo
    1581 c        enddo
    1582 c       enddo
     1636       do k=1,ntra
     1637        do j=minorig,nl
     1638         do il=1,ncum
     1639          if( (i.ge.icb(il)).and.(i.le.inb(il)).and.
     1640     :       (j.ge.(icb(il)-1)).and.(j.le.inb(il)))then
     1641            traent(il,i,j,k)=sij(il,i,j)*tra(il,i,k)
     1642     :            +(1.-sij(il,i,j))*tra(il,nk(il),k)
     1643          endif
     1644         enddo
     1645        enddo
     1646       enddo
    15831647
    15841648c
    … …  
    15901654
    15911655      do 740 il=1,ncum
    1592       if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then 
     1656      if ((i.ge.icb(il)).and.(i.le.inb(il)).and.(nent(il,i).eq.0)) then
    15931657c@      if(nent(il,i).eq.0)then
    15941658      ment(il,i,i)=m(il,i)
    1595       qent(il,i,i)=rr(il,nk(il))-ep(il,i)*clw(il,i)
    1596       uent(il,i,i)=u(il,nk(il))
    1597       vent(il,i,i)=v(il,nk(il))
     1659      qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
     1660      uent(il,i,i)=unk(il)
     1661      vent(il,i,i)=vnk(il)
    15981662      elij(il,i,i)=clw(il,i)
    15991663cMAF      sij(il,i,i)=1.0
    … …  
    16021666 740  continue
    16031667 750  continue
    1604  
    1605 c      do j=1,ntra
    1606 c       do i=minorig+1,nl
    1607 c        do il=1,ncum
    1608 c         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
    1609 c          traent(il,i,i,j)=tra(il,nk(il),j)
    1610 c         endif
    1611 c        enddo
    1612 c       enddo
    1613 c      enddo
     1668
     1669      do j=1,ntra
     1670       do i=minorig+1,nl
     1671        do il=1,ncum
     1672         if (i.ge.icb(il) .and. i.le.inb(il) .and. nent(il,i).eq.0) then
     1673          traent(il,i,i,j)=tra(il,nk(il),j)
     1674         endif
     1675        enddo
     1676       enddo
     1677      enddo
    16141678
    16151679      do 100 j=minorig,nl
    … …  
    16321696c=====================================================================
    16331697
    1634 cym      call zilch(asum,ncum*nd)
    1635 cym      call zilch(bsum,ncum*nd)
    1636 cym      call zilch(csum,ncum*nd)
    16371698      call zilch(asum,nloc*nd)
    16381699      call zilch(csum,nloc*nd)
    … …  
    16431704      enddo
    16441705
    1645       DO 789 i=minorig+1,nl 
     1706      DO 789 i=minorig+1,nl
    16461707
    16471708      num1=0
    … …  
    16551716       if ( i.ge.icb(il) .and. i.le.inb(il) ) then
    16561717        lwork(il)=(nent(il,i).ne.0)
    1657         qp=rr(il,1)-ep(il,i)*clw(il,i)
     1718        qp=qnk(il)-ep(il,i)*clw(il,i)
    16581719        anum=h(il,i)-hp(il,i)-lv(il,i)*(qp-rs(il,i))
    16591720     :           +(cpv-cpd)*t(il,i)*(qp-rr(il,i))
    … …  
    16741735      do il=1,ncum
    16751736       if ( i.ge.icb(il) .and. i.le.inb(il) .and.
    1676      :      j.ge.(icb(il)-1) .and. j.le.inb(il) 
     1737     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
    16771738     :      .and. lwork(il) ) num2=num2+1
    16781739      enddo
    … …  
    16811742      do 782 il=1,ncum
    16821743      if ( i.ge.icb(il) .and. i.le.inb(il) .and.
    1683      :      j.ge.(icb(il)-1) .and. j.le.inb(il) 
     1744     :      j.ge.(icb(il)-1) .and. j.le.inb(il)
    16841745     :      .and. lwork(il) ) then
    16851746
    … …  
    17711832        nent(il,i)=0
    17721833        ment(il,i,i)=m(il,i)
    1773         qent(il,i,i)=rr(il,1)-ep(il,i)*clw(il,i)
    1774         uent(il,i,i)=u(il,nk(il))
    1775         vent(il,i,i)=v(il,nk(il))
     1834        qent(il,i,i)=qnk(il)-ep(il,i)*clw(il,i)
     1835        uent(il,i,i)=unk(il)
     1836        vent(il,i,i)=vnk(il)
    17761837        elij(il,i,i)=clw(il,i)
    17771838cMAF        sij(il,i,i)=1.0
    … …  
    17801841      enddo ! il
    17811842
    1782 c      do j=1,ntra
    1783 c       do il=1,ncum
    1784 c        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
    1785 c     :     .and. csum(il,i).lt.m(il,i) ) then
    1786 c         traent(il,i,i,j)=tra(il,nk(il),j)
    1787 c        endif
    1788 c       enddo
    1789 c      enddo
     1843      do j=1,ntra
     1844       do il=1,ncum
     1845        if ( i.ge.icb(il) .and. i.le.inb(il) .and. lwork(il)
     1846     :     .and. csum(il,i).lt.m(il,i) ) then
     1847         traent(il,i,i,j)=tra(il,nk(il),j)
     1848        endif
     1849       enddo
     1850      enddo
    17901851789   continue
    17911852c     
    17921853c MAF: renormalisation de MENT
     1854      call zilch(zm,nloc*na)
    17931855      do jm=1,nd
    17941856        do im=1,nd
    … …  
    18211883      end
    18221884
    1823 
    1824       SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb
     1885      SUBROUTINE cv3_unsat(nloc,ncum,nd,na,ntra,icb,inb,iflag
    18251886     :              ,t,rr,rs,gz,u,v,tra,p,ph
    18261887     :              ,th,tv,lv,cpn,ep,sigp,clw
    1827      :              ,m,ment,elij,delt,plcl
    1828      :              ,mp,rp,up,vp,trap,wt,water,evap,b)
     1888     :              ,m,ment,elij,delt,plcl,coef_clos
     1889     o              ,mp,rp,up,vp,trap,wt,water,evap,b,sigd)
    18291890      implicit none
    18301891
    18311892
    18321893#include "cvthermo.h"
    1833 #include "cvparam3.h"
     1894#include "cv3param.h"
    18341895#include "cvflag.h"
    18351896
    … …  
    18381899      integer icb(nloc), inb(nloc)
    18391900      real delt, plcl(nloc)
    1840       real t(nloc,nd), rr(nloc,nd), rs(nloc,nd)
     1901      real t(nloc,nd), rr(nloc,nd), rs(nloc,nd),gz(nloc,na)
    18411902      real u(nloc,nd), v(nloc,nd)
    18421903      real tra(nloc,nd,ntra)
    18431904      real p(nloc,nd), ph(nloc,nd+1)
    1844       real th(nloc,na), gz(nloc,na)
    1845       real lv(nloc,na), ep(nloc,na), sigp(nloc,na), clw(nloc,na)
    1846       real cpn(nloc,na), tv(nloc,na)
     1905      real ep(nloc,na), sigp(nloc,na), clw(nloc,na)
     1906      real th(nloc,na),tv(nloc,na),lv(nloc,na),cpn(nloc,na)
    18471907      real m(nloc,na), ment(nloc,na,na), elij(nloc,na,na)
    1848 
     1908      real coef_clos(nloc)
     1909c
     1910c input/output
     1911      integer iflag(nloc)
     1912c
    18491913c outputs:
    18501914      real mp(nloc,na), rp(nloc,na), up(nloc,na), vp(nloc,na)
    18511915      real water(nloc,na), evap(nloc,na), wt(nloc,na)
    18521916      real trap(nloc,na,ntra)
    1853       real b(nloc,na)
     1917      real b(nloc,na), sigd(nloc)
    18541918
    18551919c local variables
    1856       integer i,j,k,il,num1
     1920      integer i,j,k,il,num1,ndp1
    18571921      real tinv, delti
    18581922      real awat, afac, afac1, afac2, bfac
    … …  
    18611925      real ampmax
    18621926      real lvcp(nloc,na)
     1927      real h(nloc,na),hm(nloc,na)
    18631928      real wdtrain(nloc)
    18641929      logical lwork(nloc)
    … …  
    18851950         enddo
    18861951        enddo
    1887 
    1888 c        do k=1,ntra
    1889 c         do i=1,nd
    1890 c          do il=1,ncum
    1891 c           trap(il,i,k)=tra(il,i,k)
    1892 c          enddo
    1893 c         enddo
    1894 c        enddo
    1895 
     1952        do k=1,ntra
     1953         do i=1,nd
     1954          do il=1,ncum
     1955           trap(il,i,k)=tra(il,i,k)
     1956          enddo
     1957         enddo
     1958        enddo
    18961959c
    18971960c   ***  check whether ep(inb)=0, if so, skip precipitating    ***
    … …  
    19051968
    19061969        call zilch(wdtrain,ncum)
    1907  
     1970
     1971c   ***  Set the fractionnal area sigd of precipitating downdraughts
     1972        do il = 1,ncum
     1973          sigd(il) = sigdz*coef_clos(il)
     1974        enddo
     1975
    19081976        DO 400 i=nl+1,1,-1
    19091977
    … …  
    19872055      if(i.eq.inb(il))afac=0.0
    19882056      afac=amax1(afac,0.0)
    1989       bfac=1./(sigd*wt(il,i))
     2057      bfac=1./(sigd(il)*wt(il,i))
    19902058c
    19912059cjyg1
    … …  
    20062074cjyg2
    20072075c
    2008       b6=bfac*50.*sigd*(ph(il,i)-ph(il,i+1))*sigt*afac
     2076cjyg----
     2077        b6 = bfac*100.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
     2078        c6 = water(il,i+1) + wdtrain(il)*bfac
     2079         revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
     2080         evap(il,i)=sigt*afac*revap
     2081         water(il,i)=revap*revap
     2082cc        print *,' i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il) ',
     2083cc     $            i,b6,c6,revap,evap(il,i),water(il,i),wdtrain(il)
     2084cc---end jyg---
     2085c
     2086c--------retour à la formulation originale d'Emanuel.
     2087      if (1.eq.0) then
     2088      b6=bfac*50.*sigd(il)*(ph(il,i)-ph(il,i+1))*sigt*afac
    20092089      c6=water(il,i+1)+bfac*wdtrain(il)
    2010      :                -50.*sigd*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
     2090     :    -50.*sigd(il)*bfac*(ph(il,i)-ph(il,i+1))*evap(il,i+1)
    20112091      if(c6.gt.0.0)then
    20122092       revap=0.5*(-b6+sqrt(b6*b6+4.*c6))
    … …  
    20152095      else
    20162096       evap(il,i)=-evap(il,i+1)
    2017      :            +0.02*(wdtrain(il)+sigd*wt(il,i)*water(il,i+1))
    2018      :                 /(sigd*(ph(il,i)-ph(il,i+1)))
     2097     :            +(wdtrain(il)+sigd(il)*wt(il,i)*water(il,i+1))
     2098     :                 /(sigd(il)*(ph(il,i)-ph(il,i+1))*50.)
    20192099      end if
    2020 c
     2100      endif
     2101ccc
    20212102c    ***  calculate precipitating downdraft mass flux under     ***
    20222103c    ***              hydrostatic approximation                 ***
    … …  
    20272108      delth=amax1(0.001,(th(il,i)-th(il,i-1)))
    20282109      if (cvflag_grav) then
    2029        mp(il,i)=100.*ginv*lvcp(il,i)*sigd*tevap
     2110       mp(il,i)=100.*ginv*lvcp(il,i)*sigd(il)*tevap
    20302111     :              *(p(il,i-1)-p(il,i))/delth
    20312112      else
    2032        mp(il,i)=10.*lvcp(il,i)*sigd*tevap*(p(il,i-1)-p(il,i))/delth
     2113       mp(il,i)=10.*lvcp(il,i)*sigd(il)*tevap
     2114     :         *(p(il,i-1)-p(il,i))/delth
    20332115      endif
    20342116c
    … …  
    20372119c    ***  and mass flux from two simultaneous differential eqns ***
    20382120c
    2039       amfac=sigd*sigd*70.0*ph(il,i)*(p(il,i-1)-p(il,i))
     2121      amfac=sigd(il)*sigd(il)*70.0*ph(il,i)*(p(il,i-1)-p(il,i))
    20402122     :          *(th(il,i)-th(il,i-1))/(tv(il,i)*th(il,i))
    20412123      amp2=abs(mp(il,i+1)*mp(il,i+1)-mp(il,i)*mp(il,i))
    20422124      if(amp2.gt.(0.1*amfac))then
    2043        xf=100.0*sigd*sigd*sigd*(ph(il,i)-ph(il,i+1))
     2125       xf=100.0*sigd(il)*sigd(il)*sigd(il)*(ph(il,i)-ph(il,i+1))
    20442126       tf=b(il,i)-5.0*(th(il,i)-th(il,i-1))*t(il,i)
    2045      :               /(lvcp(il,i)*sigd*th(il,i))
     2127     :               /(lvcp(il,i)*sigd(il)*th(il,i))
    20462128       af=xf*tf+mp(il,i+1)*mp(il,i+1)*tinv
    20472129       bf=2.*(tinv*mp(il,i+1))**3+tinv*mp(il,i+1)*xf*tf
    … …  
    20652147
    20662148       if (cvflag_grav) then
    2067 Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante: 
    2068 C il faut diviser par (mp(il,i)*sigd*grav) et non par (mp(il,i)+sigd*0.1).
     2149Cjyg : il y a vraisemblablement une erreur dans la ligne 2 suivante:
     2150C il faut diviser par (mp(il,i)*sigd(il)*grav) et non par (mp(il,i)+sigd(il)*0.1).
    20692151C Et il faut bien revoir les facteurs 100.
    20702152        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
    2071      2   /(mp(il,i)+sigd*0.1)
    2072      3   -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i))
     2153     2   /(mp(il,i)+sigd(il)*0.1)
     2154     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
     2155     : *sigd(il)*th(il,i))
    20732156       else
    20742157        b(il,i-1)=b(il,i)+100.0*(p(il,i-1)-p(il,i))*tevap
    2075      2   /(mp(il,i)+sigd*0.1)
    2076      3   -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)*sigd*th(il,i))
     2158     2   /(mp(il,i)+sigd(il)*0.1)
     2159     3 -10.0*(th(il,i)-th(il,i-1))*t(il,i)/(lvcp(il,i)
     2160     : *sigd(il)*th(il,i))
    20772161       endif
    20782162       b(il,i-1)=amax1(b(il,i-1),0.0)
    … …  
    20892173c    ***       between cloud base and the surface                   ***
    20902174c
    2091       if(p(il,i).gt.p(il,icb(il)))then
    2092        mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
     2175c
     2176cc      if(p(il,i).gt.p(il,icb(il)))then
     2177cc       mp(il,i)=mp(il,icb(il))*(p(il,1)-p(il,i))/(p(il,1)-p(il,icb(il)))
     2178cc      endif
     2179      if(ph(il,i) .gt. 0.9*plcl(il)) then
     2180       mp(il,i) = mp(il,i)*(ph(il,1)-ph(il,i))/
     2181     $                     (ph(il,1)-0.9*plcl(il))
    20932182      endif
    20942183
    … …  
    21072196       if (cvflag_grav) then
    21082197        rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1))
    2109      :   +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1))
     2198     :   +100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))
    21102199     :                     *(evap(il,i+1)+evap(il,i))
    21112200       else
    21122201        rp(il,i)=rp(il,i+1)*mp(il,i+1)+rr(il,i)*(mp(il,i)-mp(il,i+1))
    2113      :   +5.*sigd*(ph(il,i)-ph(il,i+1))
     2202     :   +5.*sigd(il)*(ph(il,i)-ph(il,i+1))
    21142203     :                      *(evap(il,i+1)+evap(il,i))
    21152204       endif
    … …  
    21202209      vp(il,i)=vp(il,i)/mp(il,i)
    21212210
    2122 c      do j=1,ntra
    2123 c      trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
    2124 ctestmaf     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
    2125 c     :            +tra(il,i,j)*(mp(il,i)-mp(il,i+1))
    2126 c      trap(il,i,j)=trap(il,i,j)/mp(il,i)
    2127 c      end do
     2211      do j=1,ntra
     2212      trap(il,i,j)=trap(il,i+1,j)*mp(il,i+1)
     2213     :            +trap(il,i,j)*(mp(il,i)-mp(il,i+1))
     2214      trap(il,i,j)=trap(il,i,j)/mp(il,i)
     2215      end do
    21282216
    21292217      else
    … …  
    21322220        if (cvflag_grav) then
    21332221         rp(il,i)=rp(il,i+1)
    2134      :            +100.*ginv*0.5*sigd*(ph(il,i)-ph(il,i+1))
     2222     :            +100.*ginv*0.5*sigd(il)*(ph(il,i)-ph(il,i+1))
    21352223     :            *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
    21362224        else
    21372225         rp(il,i)=rp(il,i+1)
    2138      :           +5.*sigd*(ph(il,i)-ph(il,i+1))
     2226     :           +5.*sigd(il)*(ph(il,i)-ph(il,i+1))
    21392227     :           *(evap(il,i+1)+evap(il,i))/mp(il,i+1)
    21402228        endif
    … …  
    21422230       vp(il,i)=vp(il,i+1)
    21432231
    2144 c       do j=1,ntra
    2145 c       trap(il,i,j)=trap(il,i+1,j)
    2146 c       end do
     2232       do j=1,ntra
     2233       trap(il,i,j)=trap(il,i+1,j)
     2234       end do
    21472235
    21482236       endif
    … …  
    21602248       end
    21612249
    2162       SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra 
     2250      SUBROUTINE cv3_yield(nloc,ncum,nd,na,ntra
    21632251     :                    ,icb,inb,delt
    2164      :                    ,t,rr,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
     2252     :                    ,t,rr,t_wake,rr_wake,u,v,tra
     2253     :                    ,gz,p,ph,h,hp,lv,cpn,th
    21652254     :                    ,ep,clw,m,tp,mp,rp,up,vp,trap
    2166      :                    ,wt,water,evap,b
    2167      :                    ,ment,qent,uent,vent,nent,elij,traent,sig
    2168      :                    ,tv,tvp
    2169      :                    ,iflag,precip,VPrecip,ft,fr,fu,fv,ftra
    2170      :                    ,upwd,dnwd,dnwd0,ma,mike,tls,tps,qcondc,wd)
     2255     :                    ,wt,water,evap,b,sigd
     2256     :                    ,ment,qent,hent,iflag_mix,uent,vent
     2257     :                    ,nent,elij,traent,sig
     2258     :                    ,tv,tvp,wghti
     2259     :                    ,iflag,precip,Vprecip,ft,fr,fu,fv,ftra
     2260     :                    ,cbmf,upwd,dnwd,dnwd0,ma,mip
     2261     :                    ,tls,tps,qcondc,wd
     2262     :                    ,ftd,fqd)
     2263     
    21712264      implicit none
    21722265
    21732266#include "cvthermo.h"
    2174 #include "cvparam3.h"
     2267#include "cv3param.h"
    21752268#include "cvflag.h"
    21762269#include "conema3.h"
    21772270
    21782271c inputs:
     2272c      print*,'cv3_yield apres include'
     2273      integer iflag_mix
    21792274      integer ncum,nd,na,ntra,nloc
    21802275      integer icb(nloc), inb(nloc)
    21812276      real delt
    21822277      real t(nloc,nd), rr(nloc,nd), u(nloc,nd), v(nloc,nd)
     2278      real t_wake(nloc,nd), rr_wake(nloc,nd)
    21832279      real tra(nloc,nd,ntra), sig(nloc,nd)
    21842280      real gz(nloc,na), ph(nloc,nd+1), h(nloc,na), hp(nloc,na)
    … …  
    21872283      real m(nloc,na), mp(nloc,na), rp(nloc,na), up(nloc,na)
    21882284      real vp(nloc,na), wt(nloc,nd), trap(nloc,nd,ntra)
    2189       real water(nloc,na), evap(nloc,na), b(nloc,na)
     2285      real water(nloc,na), evap(nloc,na), b(nloc,na), sigd(nloc)
    21902286      real ment(nloc,na,na), qent(nloc,na,na), uent(nloc,na,na)
    2191 cym      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
    2192       real vent(nloc,na,na), elij(nloc,na,na)
    2193       integer nent(nloc,na)
     2287      real hent(nloc,na,na)
     2288      real vent(nloc,na,na), nent(nloc,na), elij(nloc,na,na)
    21942289      real traent(nloc,na,na,ntra)
    2195       real tv(nloc,nd), tvp(nloc,nd)
    2196 
     2290      real tv(nloc,nd), tvp(nloc,nd), wghti(nloc,nd)
     2291c      print*,'cv3_yield declarations 1'
    21972292c input/output:
    21982293      integer iflag(nloc)
    … …  
    22002295c outputs:
    22012296      real precip(nloc)
    2202       real VPrecip(nloc,nd+1)
    22032297      real ft(nloc,nd), fr(nloc,nd), fu(nloc,nd), fv(nloc,nd)
     2298      real ftd(nloc,nd), fqd(nloc,nd)
    22042299      real ftra(nloc,nd,ntra)
    22052300      real upwd(nloc,nd), dnwd(nloc,nd), ma(nloc,nd)
    2206       real dnwd0(nloc,nd), mike(nloc,nd)
     2301      real dnwd0(nloc,nd), mip(nloc,nd)
     2302      real Vprecip(nloc,nd)
    22072303      real tls(nloc,nd), tps(nloc,nd)
    22082304      real qcondc(nloc,nd)                               ! cld
    22092305      real wd(nloc)                                      ! gust
    2210 
     2306      real cbmf(nloc)
     2307c      print*,'cv3_yield declarations 2'
    22112308c local variables:
    22122309      integer i,k,il,n,j,num1
    2213       real rat, awat, delti
     2310      real rat, delti
    22142311      real ax, bx, cx, dx, ex
    22152312      real cpinv, rdcp, dpinv
     2313      real awat(nloc)
    22162314      real lvcp(nloc,na), mke(nloc,na)
    22172315      real am(nloc), work(nloc), ad(nloc), amp1(nloc)
    … …  
    22222320      real siga(nloc,nd), sax(nloc,nd), mac(nloc,nd)      ! cld
    22232321
    2224 
     2322c      print*,'cv3_yield declarations 3'
    22252323c-------------------------------------------------------------
    22262324
    … …  
    22282326
    22292327      delti = 1.0/delt
    2230 
     2328c      print*,'cv3_yield initialisation delt', delt
     2329cprecip,Vprecip,ft,fr,fu,fv,ftra
     2330c     :                    ,cbmf,upwd,dnwd,dnwd0,ma,mip
     2331c     :                    ,tls,tps,qcondc,wd
     2332c     :                    ,ftd,fqd  )
    22312333      do il=1,ncum
    22322334       precip(il)=0.0
     2335c       Vprecip(il,nd+1)=0.0
    22332336       wd(il)=0.0     ! gust
    2234        VPrecip(il,nd+1)=0.
    22352337      enddo
    22362338
    22372339      do i=1,nd
    22382340       do il=1,ncum
    2239          VPrecip(il,i)=0.0
     2341         Vprecip(il,i)=0.0
    22402342         ft(il,i)=0.0
    22412343         fr(il,i)=0.0
    22422344         fu(il,i)=0.0
    22432345         fv(il,i)=0.0
     2346         upwd(il,i)=0.0
     2347         dnwd(il,i)=0.0
     2348         dnwd0(il,i)=0.0
     2349         mip(il,i)=0.0
     2350         ftd(il,i)=0.0
     2351         fqd(il,i)=0.0
    22442352         qcondc(il,i)=0.0                                ! cld
    22452353         qcond(il,i)=0.0                                 ! cld
    … …  
    22472355       enddo
    22482356      enddo
    2249 
    2250 c      do j=1,ntra
    2251 c       do i=1,nd
    2252 c        do il=1,ncum
    2253 c          ftra(il,i,j)=0.0
    2254 c        enddo
    2255 c       enddo
    2256 c      enddo
    2257 
     2357c       print*,'cv3_yield initialisation 2'
     2358      do j=1,ntra
     2359       do i=1,nd
     2360        do il=1,ncum
     2361          ftra(il,i,j)=0.0
     2362        enddo
     2363       enddo
     2364      enddo
     2365c       print*,'cv3_yield initialisation 3'
    22582366      do i=1,nl
    22592367       do il=1,ncum
    … …  
    22662374c   ***  calculate surface precipitation in mm/day     ***
    22672375c
    2268       do il=1,ncum 
    2269        if(ep(il,inb(il)).ge.0.0001)then
     2376      do il=1,ncum
     2377       if(ep(il,inb(il)).ge.0.0001 .and. iflag(il) .le. 1)then
    22702378        if (cvflag_grav) then
    2271          precip(il)=wt(il,1)*sigd*water(il,1)*86400.*1000./(rowl*grav)
     2379           precip(il)=wt(il,1)*sigd(il)*water(il,1)*86400.*1000.
     2380     :               /(rowl*grav)
    22722381        else
    2273          precip(il)=wt(il,1)*sigd*water(il,1)*8640.
     2382         precip(il)=wt(il,1)*sigd(il)*water(il,1)*8640.
    22742383        endif
    2275        endif 
    2276       enddo 
    2277 
    2278 C   ***  CALCULATE VERTICAL PROFILE OF  PRECIPITATIONs IN kg/m2/s  ===
     2384       endif
     2385      enddo
     2386c      print*,'cv3_yield apres calcul precip'
     2387
    22792388C
    2280 c MAF rajout pour lessivage
    2281        do k=1,nl
    2282          do il=1,ncum
    2283           if (k.le.inb(il)) then
    2284             if (cvflag_grav) then
    2285              VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/grav
    2286             else
    2287              VPrecip(il,k) = wt(il,k)*sigd*water(il,k)/10.
    2288             endif
    2289           endif
    2290          end do
    2291        end do
     2389C   ===  calculate vertical profile of  precipitation in kg/m2/s  ===
     2390C
     2391      do i = 1,nl
     2392      do il=1,ncum
     2393       if(ep(il,inb(il)).ge.0.0001 .and. i.le.inb(il)
     2394     :    .and. iflag(il) .le. 1)then
     2395        if (cvflag_grav) then
     2396           VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/grav
     2397        else
     2398           VPrecip(il,i) = wt(il,i)*sigd(il)*water(il,i)/10.
     2399        endif
     2400       endif
     2401      enddo
     2402      enddo
    22922403C
    22932404c
    … …  
    22972408c!      do il=1,ncum
    22982409c!        wd(il)=betad*abs(mp(il,icb(il)))*0.01*rrd*t(il,icb(il))
    2299 c!     :                                  /(sigd*p(il,icb(il)))
     2410c!     :                                  /(sigd(il)*p(il,icb(il)))
    23002411c!      enddo
    23012412
    … …  
    23062417      do il=1,ncum
    23072418       work(il)=1.0/(ph(il,1)-ph(il,2))
    2308        am(il)=0.0
     2419       cbmf(il)=0.0
    23092420      enddo
    23102421
    23112422      do k=2,nl
    23122423       do il=1,ncum
    2313         if (k.le.inb(il)) then
    2314          am(il)=am(il)+m(il,k)
     2424        if (k.ge.icb(il)) then
     2425         cbmf(il)=cbmf(il)+m(il,k)
    23152426        endif
    23162427       enddo
    23172428      enddo
    23182429
     2430c      print*,'cv3_yield avant ft'
     2431c AM is the part of cbmf taken from the first level
    23192432      do il=1,ncum
    2320 
     2433        am(il)=cbmf(il)*wghti(il,1)
     2434      enddo
     2435c
     2436      do il=1,ncum
     2437        if (iflag(il) .le. 1) then
    23212438c convect3      if((0.1*dpinv*am).ge.delti)iflag(il)=4
     2439cjyg  Correction pour conserver l'eau
     2440ccc       ft(il,1)=-0.5*lvcp(il,1)*sigd(il)*(evap(il,1)+evap(il,2)) !precip
     2441       ft(il,1)=-lvcp(il,1)*sigd(il)*evap(il,1)                  !precip
     2442
    23222443      if (cvflag_grav) then
     2444        ft(il,1)=ft(il,1)-0.009*grav*sigd(il)*mp(il,2)
     2445     :                              *t_wake(il,1)*b(il,1)*work(il)
     2446      else
     2447        ft(il,1)=ft(il,1)-0.09*sigd(il)*mp(il,2)
     2448     :                              *t_wake(il,1)*b(il,1)*work(il)
     2449      endif
     2450
     2451      ft(il,1)=ft(il,1)+0.01*sigd(il)*wt(il,1)*(cl-cpd)*water(il,2)
     2452     :     *(t_wake(il,2)-t_wake(il,1))*work(il)/cpn(il,1)
     2453
     2454      ftd(il,1) = ft(il,1)                        ! fin precip
     2455
     2456      if (cvflag_grav) then                  !sature
    23232457      if((0.01*grav*work(il)*am(il)).ge.delti)iflag(il)=1!consist vect
    2324        ft(il,1)=0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)
     2458       ft(il,1)=ft(il,1)+0.01*grav*work(il)*am(il)*(t(il,2)-t(il,1)
    23252459     :            +(gz(il,2)-gz(il,1))/cpn(il,1))
    23262460      else
    23272461       if((0.1*work(il)*am(il)).ge.delti)iflag(il)=1 !consistency vect
    2328        ft(il,1)=0.1*work(il)*am(il)*(t(il,2)-t(il,1)
     2462       ft(il,1)=ft(il,1)+0.1*work(il)*am(il)*(t(il,2)-t(il,1)
    23292463     :            +(gz(il,2)-gz(il,1))/cpn(il,1))
    23302464      endif
    2331 
    2332       ft(il,1)=ft(il,1)-0.5*lvcp(il,1)*sigd*(evap(il,1)+evap(il,2))
    2333 
    2334       if (cvflag_grav) then
    2335        ft(il,1)=ft(il,1)-0.009*grav*sigd*mp(il,2)
    2336      :                             *t(il,1)*b(il,1)*work(il)
    2337       else
    2338        ft(il,1)=ft(il,1)-0.09*sigd*mp(il,2)*t(il,1)*b(il,1)*work(il)
    2339       endif
    2340 
    2341       ft(il,1)=ft(il,1)+0.01*sigd*wt(il,1)*(cl-cpd)*water(il,2)*(t(il,2)
    2342      :-t(il,1))*work(il)/cpn(il,1)
    2343 
    2344       if (cvflag_grav) then
     2465      endif  ! iflag
     2466      enddo
     2467
     2468       do j=2,nl
     2469      IF (iflag_mix .gt. 0) then
     2470        do il=1,ncum
     2471         if (j.le.inb(il) .and. iflag(il) .le. 1) then
     2472         if (cvflag_grav) then
     2473          ft(il,1)=ft(il,1)
     2474     :       +0.01*grav*work(il)*ment(il,j,1)*(hent(il,j,1)-h(il,1)
     2475     :       +t(il,1)*(cpv-cpd)*(rr(il,1)-Qent(il,j,1)))*cpinv
     2476         else
     2477          ft(il,1)=ft(il,1)
     2478     :       +0.1*work(il)*ment(il,j,1)*(hent(il,j,1)-h(il,1)
     2479     :       +t(il,1)*(cpv-cpd)*(rr(il,1)-Qent(il,j,1)))*cpinv
     2480         endif  ! cvflag_grav
     2481        endif ! j
     2482       enddo
     2483       ENDIF
     2484        enddo
     2485         ! fin sature
     2486
     2487
     2488      do il=1,ncum
     2489        if (iflag(il) .le. 1) then
     2490          if (cvflag_grav) then
    23452491Cjyg1  Correction pour mieux conserver l'eau (conformite avec CONVECT4.3)
    2346 c (sb: pour l'instant, on ne fait que le chgt concernant grav, pas evap)
    2347        fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr(il,1))*work(il)
    2348      :          +sigd*0.5*(evap(il,1)+evap(il,2))
    2349 c+tard     :          +sigd*evap(il,1)
    2350 
    2351        fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)
     2492       fr(il,1)=0.01*grav*mp(il,2)*(rp(il,2)-rr_wake(il,1))*work(il)
     2493     :          +sigd(il)*evap(il,1)
     2494ccc     :          +sigd(il)*0.5*(evap(il,1)+evap(il,2))
     2495
     2496       fqd(il,1)=fr(il,1)     !precip
     2497
     2498       fr(il,1)=fr(il,1)+0.01*grav*am(il)*(rr(il,2)-rr(il,1))*work(il)  !sature
    23522499
    23532500       fu(il,1)=fu(il,1)+0.01*grav*work(il)*(mp(il,2)*(up(il,2)-u(il,1))
    … …  
    23562503     :         +am(il)*(v(il,2)-v(il,1)))
    23572504      else  ! cvflag_grav
    2358        fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr(il,1))*work(il)
    2359      :          +sigd*0.5*(evap(il,1)+evap(il,2))
     2505       fr(il,1)=0.1*mp(il,2)*(rp(il,2)-rr_wake(il,1))*work(il)
     2506     :          +sigd(il)*evap(il,1)
     2507ccc     :          +sigd(il)*0.5*(evap(il,1)+evap(il,2))
     2508       fqd(il,1)=fr(il,1)  !precip
    23602509       fr(il,1)=fr(il,1)+0.1*am(il)*(rr(il,2)-rr(il,1))*work(il)
    23612510       fu(il,1)=fu(il,1)+0.1*work(il)*(mp(il,2)*(up(il,2)-u(il,1))
    … …  
    23632512       fv(il,1)=fv(il,1)+0.1*work(il)*(mp(il,2)*(vp(il,2)-v(il,1))
    23642513     :         +am(il)*(v(il,2)-v(il,1)))
    2365       endif ! cvflag_grav
    2366 
     2514         endif ! cvflag_grav
     2515       endif  ! iflag
    23672516      enddo ! il
    23682517
    2369 c      do j=1,ntra
    2370 c       do il=1,ncum
    2371 c        if (cvflag_grav) then
    2372 c         ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
    2373 c     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2374 c     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2375 c        else
    2376 c         ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
    2377 c     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
    2378 c     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
    2379 c        endif
    2380 c       enddo
    2381 c      enddo
    2382 
    2383       do j=2,nl
     2518
     2519      do j=1,ntra
    23842520       do il=1,ncum
    2385         if (j.le.inb(il)) then
     2521        if (iflag(il) .le. 1) then
     2522        if (cvflag_grav) then
     2523         ftra(il,1,j)=ftra(il,1,j)+0.01*grav*work(il)
     2524     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     2525     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     2526        else
     2527         ftra(il,1,j)=ftra(il,1,j)+0.1*work(il)
     2528     :                     *(mp(il,2)*(trap(il,2,j)-tra(il,1,j))
     2529     :             +am(il)*(tra(il,2,j)-tra(il,1,j)))
     2530        endif
     2531        endif  ! iflag
     2532       enddo
     2533      enddo
     2534
     2535       do j=2,nl
     2536       do il=1,ncum
     2537        if (j.le.inb(il) .and. iflag(il) .le. 1) then
    23862538         if (cvflag_grav) then
    23872539          fr(il,1)=fr(il,1)
    … …  
    23972549     :         +0.1*work(il)*ment(il,j,1)*(uent(il,j,1)-u(il,1))
    23982550          fv(il,1)=fv(il,1)
    2399      :         +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))
     2551     :         +0.1*work(il)*ment(il,j,1)*(vent(il,j,1)-v(il,1))  ! fin sature
    24002552         endif  ! cvflag_grav
    24012553        endif ! j
    … …  
    24032555      enddo
    24042556
    2405 c      do k=1,ntra
    2406 c       do j=2,nl
    2407 c        do il=1,ncum
    2408 c         if (j.le.inb(il)) then
    2409 
    2410 c          if (cvflag_grav) then
    2411 c           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
    2412 c     :                *(traent(il,j,1,k)-tra(il,1,k))
    2413 c          else
    2414 c           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
    2415 c     :                *(traent(il,j,1,k)-tra(il,1,k))
    2416 c          endif
    2417 
    2418 c         endif
    2419 c        enddo
    2420 c       enddo
    2421 c      enddo
    2422 
     2557      do k=1,ntra
     2558       do j=2,nl
     2559        do il=1,ncum
     2560         if (j.le.inb(il) .and. iflag(il) .le. 1) then
     2561
     2562          if (cvflag_grav) then
     2563           ftra(il,1,k)=ftra(il,1,k)+0.01*grav*work(il)*ment(il,j,1)
     2564     :                *(traent(il,j,1,k)-tra(il,1,k))
     2565          else
     2566           ftra(il,1,k)=ftra(il,1,k)+0.1*work(il)*ment(il,j,1)
     2567     :                *(traent(il,j,1,k)-tra(il,1,k))
     2568          endif
     2569
     2570         endif
     2571        enddo
     2572       enddo
     2573      enddo
     2574c      print*,'cv3_yield apres ft'
    24232575c
    24242576c   ***  calculate tendencies of potential temperature and mixing ratio  ***
    … …  
    24332585       num1=0
    24342586       do il=1,ncum
    2435         if(i.le.inb(il))num1=num1+1
     2587        if(i.le.inb(il) .and. iflag(il) .le. 1)num1=num1+1
    24362588       enddo
    24372589       if(num1.le.0)go to 500
    … …  
    24402592       call zilch(ad,ncum)
    24412593
    2442       do 440 k=i+1,nl+1
     2594      do 440 k=1,nl+1
    24432595       do 441 il=1,ncum
    2444         if (i.le.inb(il) .and. k.le.(inb(il)+1)) then
    2445          amp1(il)=amp1(il)+m(il,k)
     2596        if(i.ge.icb(il)) then
     2597          if(k.ge.i+1.and. k.le.(inb(il)+1)) then
     2598            amp1(il)=amp1(il)+m(il,k)
     2599          endif
     2600         else
     2601c AMP1 is the part of cbmf taken from layers I and lower
     2602          if(k.le.i) then
     2603           amp1(il)=amp1(il)+cbmf(il)*wghti(il,k)
     2604          endif
    24462605        endif
    24472606 441   continue
    … …  
    24692628 
    24702629      do 1350 il=1,ncum
    2471       if (i.le.inb(il)) then
     2630      if (i.le.inb(il) .and. iflag(il) .le. 1) then
    24722631       dpinv=1.0/(ph(il,i)-ph(il,i+1))
    24732632       cpinv=1.0/cpn(il,i)
    … …  
    24802639      endif
    24812640
     2641       ! precip
     2642ccc       ft(il,i)= -0.5*sigd(il)*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
     2643       ft(il,i)= -sigd(il)*lvcp(il,i)*evap(il,i)
     2644        rat=cpn(il,i-1)*cpinv
     2645c
    24822646      if (cvflag_grav) then
    2483        ft(il,i)=0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i)
     2647       ft(il,i)=ft(il,i)-0.009*grav*sigd(il)
     2648     :   *(mp(il,i+1)*t_wake(il,i)*b(il,i)
     2649     :   -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
     2650       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
     2651     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     2652       ftd(il,i)=ft(il,i)
     2653        ! fin precip
     2654c
     2655           ! sature
     2656       ft(il,i)=ft(il,i)+0.01*grav*dpinv*(amp1(il)*(t(il,i+1)-t(il,i)
    24842657     :    +(gz(il,i+1)-gz(il,i))*cpinv)
    24852658     :    -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
    2486      :    -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
    2487        rat=cpn(il,i-1)*cpinv
    2488        ft(il,i)=ft(il,i)-0.009*grav*sigd*(mp(il,i+1)*t(il,i)*b(il,i)
    2489      :   -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
     2659
     2660c
     2661      IF (iflag_mix .eq. 0) then
    24902662       ft(il,i)=ft(il,i)+0.01*grav*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i)
    24912663     :    +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
     2664      ENDIF
     2665c
    24922666      else  ! cvflag_grav
    2493        ft(il,i)=0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i)
     2667       ft(il,i)=ft(il,i)-0.09*sigd(il)
     2668     :   *(mp(il,i+1)*t_wake(il,i)*b(il,i)
     2669     :   -mp(il,i)*t_wake(il,i-1)*rat*b(il,i-1))*dpinv
     2670       ft(il,i)=ft(il,i)+0.01*sigd(il)*wt(il,i)*(cl-cpd)*water(il,i+1)
     2671     :           *(t_wake(il,i+1)-t_wake(il,i))*dpinv*cpinv
     2672       ftd(il,i)=ft(il,i)
     2673        ! fin precip
     2674c
     2675           ! sature
     2676       ft(il,i)=ft(il,i)+0.1*dpinv*(amp1(il)*(t(il,i+1)-t(il,i)
    24942677     :    +(gz(il,i+1)-gz(il,i))*cpinv)
    24952678     :    -ad(il)*(t(il,i)-t(il,i-1)+(gz(il,i)-gz(il,i-1))*cpinv))
    2496      :    -0.5*sigd*lvcp(il,i)*(evap(il,i)+evap(il,i+1))
    2497        rat=cpn(il,i-1)*cpinv
    2498        ft(il,i)=ft(il,i)-0.09*sigd*(mp(il,i+1)*t(il,i)*b(il,i)
    2499      :   -mp(il,i)*t(il,i-1)*rat*b(il,i-1))*dpinv
     2679
     2680c
     2681      IF (iflag_mix .eq. 0) then
    25002682       ft(il,i)=ft(il,i)+0.1*dpinv*ment(il,i,i)*(hp(il,i)-h(il,i)
    25012683     :    +t(il,i)*(cpv-cpd)*(rr(il,i)-qent(il,i,i)))*cpinv
     2684      ENDIF
    25022685      endif ! cvflag_grav
    25032686
    25042687
    2505       ft(il,i)=ft(il,i)+0.01*sigd*wt(il,i)*(cl-cpd)*water(il,i+1)
    2506      :           *(t(il,i+1)-t(il,i))*dpinv*cpinv
     2688        if (cvflag_grav) then
     2689c sb: on ne fait pas encore la correction permettant de mieux
     2690c conserver l'eau:
     2691c jyg: correction permettant de mieux conserver l'eau:
     2692ccc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
     2693         fr(il,i)=sigd(il)*evap(il,i)
     2694     :        +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i))
     2695     :        -mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
     2696         fqd(il,i)=fr(il,i)    ! precip
     2697
     2698         fu(il,i)=0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i))
     2699     :             -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
     2700         fv(il,i)=0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i))
     2701     :             -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
     2702        else  ! cvflag_grav
     2703ccc         fr(il,i)=0.5*sigd(il)*(evap(il,i)+evap(il,i+1))
     2704         fr(il,i)=sigd(il)*evap(il,i)
     2705     :        +0.1*(mp(il,i+1)*(rp(il,i+1)-rr_wake(il,i))
     2706     :             -mp(il,i)*(rp(il,i)-rr_wake(il,i-1)))*dpinv
     2707         fqd(il,i)=fr(il,i)    ! precip
     2708
     2709         fu(il,i)=0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i))
     2710     :             -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
     2711         fv(il,i)=0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i))
     2712     :             -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
     2713        endif ! cvflag_grav
     2714
    25072715
    25082716      if (cvflag_grav) then
    2509        fr(il,i)=0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i))
     2717       fr(il,i)=fr(il,i)+0.01*grav*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i))
    25102718     :           -ad(il)*(rr(il,i)-rr(il,i-1)))
    25112719       fu(il,i)=fu(il,i)+0.01*grav*dpinv*(amp1(il)*(u(il,i+1)-u(il,i))
    … …  
    25142722     :             -ad(il)*(v(il,i)-v(il,i-1)))
    25152723      else  ! cvflag_grav
    2516        fr(il,i)=0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i))
     2724       fr(il,i)=fr(il,i)+0.1*dpinv*(amp1(il)*(rr(il,i+1)-rr(il,i))
    25172725     :           -ad(il)*(rr(il,i)-rr(il,i-1)))
    25182726       fu(il,i)=fu(il,i)+0.1*dpinv*(amp1(il)*(u(il,i+1)-u(il,i))
    … …  
    252527331350  continue
    25262734
    2527 c      do k=1,ntra
    2528 c       do il=1,ncum
    2529 c        if (i.le.inb(il)) then
    2530 c         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2531 c         cpinv=1.0/cpn(il,i)
    2532 c         if (cvflag_grav) then
    2533 c           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
    2534 c     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    2535 c     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    2536 c         else
    2537 c           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
    2538 c     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
    2539 c     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
    2540 c         endif
    2541 c        endif
    2542 c       enddo
    2543 c      enddo
    2544 
    2545       do 480 k=1,i-1
    2546        do 1370 il=1,ncum
    2547         if (i.le.inb(il)) then
     2735      do k=1,ntra
     2736       do il=1,ncum
     2737        if (i.le.inb(il) .and. iflag(il) .le. 1) then
    25482738         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    25492739         cpinv=1.0/cpn(il,i)
    2550 
    2551       awat=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
    2552       awat=amax1(awat,0.0)
    2553 
    2554       if (cvflag_grav) then
     2740         if (cvflag_grav) then
     2741           ftra(il,i,k)=ftra(il,i,k)+0.01*grav*dpinv
     2742     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     2743     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     2744         else
     2745           ftra(il,i,k)=ftra(il,i,k)+0.1*dpinv
     2746     :         *(amp1(il)*(tra(il,i+1,k)-tra(il,i,k))
     2747     :           -ad(il)*(tra(il,i,k)-tra(il,i-1,k)))
     2748         endif
     2749        endif
     2750       enddo
     2751      enddo
     2752
     2753      do 480 k=1,i-1
     2754c
     2755       do il = 1,ncum
     2756        awat(il)=elij(il,k,i)-(1.-ep(il,i))*clw(il,i)
     2757        awat(il)=amax1(awat(il),0.0)
     2758       enddo
     2759c
     2760      IF (iflag_mix .ne. 0) then
     2761       do il=1,ncum
     2762        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     2763         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2764         cpinv=1.0/cpn(il,i)
     2765       if (cvflag_grav) then
     2766      ft(il,i)=ft(il,i)
     2767     :       +0.01*grav*dpinv*ment(il,k,i)*(hent(il,k,i)-h(il,i)
     2768     :       +t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
     2769
     2770c
     2771c
     2772       else
     2773      ft(il,i)=ft(il,i)
     2774     :       +0.1*dpinv*ment(il,k,i)*(hent(il,k,i)-h(il,i)
     2775     :       +t(il,i)*(cpv-cpd)*(rr(il,i)+awat(il)-Qent(il,k,i)))*cpinv
     2776       endif  !cvflag_grav
     2777       endif  ! i
     2778       enddo
     2779      ENDIF
     2780c
     2781       do 1370 il=1,ncum
     2782        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     2783         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2784         cpinv=1.0/cpn(il,i)
     2785       if (cvflag_grav) then
    25552786      fr(il,i)=fr(il,i)
    2556      :   +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
     2787     :   +0.01*grav*dpinv*ment(il,k,i)*(qent(il,k,i)-awat(il)-rr(il,i))
    25572788      fu(il,i)=fu(il,i)
    25582789     :         +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
    … …  
    25612792      else  ! cvflag_grav
    25622793      fr(il,i)=fr(il,i)
    2563      :   +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat-rr(il,i))
     2794     :   +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-awat(il)-rr(il,i))
    25642795      fu(il,i)=fu(il,i)
    25652796     :   +0.01*grav*dpinv*ment(il,k,i)*(uent(il,k,i)-u(il,i))
    … …  
    25692800
    25702801c (saturated updrafts resulting from mixing)        ! cld
    2571         qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat) ! cld
     2802        qcond(il,i)=qcond(il,i)+(elij(il,k,i)-awat(il)) ! cld
    25722803        nqcond(il,i)=nqcond(il,i)+1.                ! cld
    25732804      endif ! i
    … …  
    25752806480   continue
    25762807
    2577 c      do j=1,ntra
    2578 c       do k=1,i-1
    2579 c        do il=1,ncum
    2580 c         if (i.le.inb(il)) then
    2581 c          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2582 c          cpinv=1.0/cpn(il,i)
    2583 c          if (cvflag_grav) then
    2584 c           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    2585 c     :        *(traent(il,k,i,j)-tra(il,i,j))
    2586 c          else
    2587 c           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    2588 c     :        *(traent(il,k,i,j)-tra(il,i,j))
    2589 c          endif
    2590 c         endif
    2591 c        enddo
    2592 c       enddo
    2593 c      enddo
     2808      do j=1,ntra
     2809       do k=1,i-1
     2810        do il=1,ncum
     2811         if (i.le.inb(il) .and. iflag(il) .le. 1) then
     2812          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2813          cpinv=1.0/cpn(il,i)
     2814          if (cvflag_grav) then
     2815           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     2816     :        *(traent(il,k,i,j)-tra(il,i,j))
     2817          else
     2818           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     2819     :        *(traent(il,k,i,j)-tra(il,i,j))
     2820          endif
     2821         endif
     2822        enddo
     2823       enddo
     2824      enddo
    25942825
    25952826      do 490 k=i,nl+1
     2827c
     2828      IF (iflag_mix .ne. 0) then
     2829       do il=1,ncum
     2830        if (i.le.inb(il) .and. k.le.inb(il)
     2831     $               .and. iflag(il) .le. 1) then
     2832         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2833         cpinv=1.0/cpn(il,i)
     2834       if (cvflag_grav) then
     2835      ft(il,i)=ft(il,i)
     2836     :       +0.01*grav*dpinv*ment(il,k,i)*(hent(il,k,i)-h(il,i)
     2837     :       +t(il,i)*(cpv-cpd)*(rr(il,i)-Qent(il,k,i)))*cpinv
     2838c
     2839c
     2840       else
     2841      ft(il,i)=ft(il,i)
     2842     :       +0.1*dpinv*ment(il,k,i)*(hent(il,k,i)-h(il,i)
     2843     :       +t(il,i)*(cpv-cpd)*(rr(il,i)-Qent(il,k,i)))*cpinv
     2844       endif  !cvflag_grav
     2845       endif  ! i
     2846       enddo
     2847      ENDIF
     2848c
    25962849       do 1380 il=1,ncum
    2597         if (i.le.inb(il) .and. k.le.inb(il)) then
     2850        if (i.le.inb(il) .and. k.le.inb(il)
     2851     $               .and. iflag(il) .le. 1) then
    25982852         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    25992853         cpinv=1.0/cpn(il,i)
    … …  
    26062860         fv(il,i)=fv(il,i)
    26072861     :         +0.01*grav*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
    2608          else  ! cvflag_grav 
     2862         else  ! cvflag_grav
    26092863         fr(il,i)=fr(il,i)
    26102864     :         +0.1*dpinv*ment(il,k,i)*(qent(il,k,i)-rr(il,i))
    … …  
    26132867         fv(il,i)=fv(il,i)
    26142868     :         +0.1*dpinv*ment(il,k,i)*(vent(il,k,i)-v(il,i))
    2615          endif ! cvflag_grav 
     2869         endif ! cvflag_grav
    26162870        endif ! i and k
    261728711380   continue
    26182872490   continue
    26192873
    2620 c      do j=1,ntra
    2621 c       do k=i,nl+1
    2622 c        do il=1,ncum
    2623 c         if (i.le.inb(il) .and. k.le.inb(il)) then
    2624 c          dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2625 c          cpinv=1.0/cpn(il,i)
    2626 c          if (cvflag_grav) then
    2627 c           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
    2628 c     :         *(traent(il,k,i,j)-tra(il,i,j))
    2629 c          else
    2630 c           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
    2631 c     :             *(traent(il,k,i,j)-tra(il,i,j))
    2632 c          endif
    2633 c         endif ! i and k
    2634 c        enddo
    2635 c       enddo
    2636 c      enddo
    2637 
    2638       do 1400 il=1,ncum
    2639        if (i.le.inb(il)) then
    2640         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2641         cpinv=1.0/cpn(il,i)
    2642 
    2643         if (cvflag_grav) then
    2644 c sb: on ne fait pas encore la correction permettant de mieux
    2645 c conserver l'eau:
    2646          fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1))
    2647      :        +0.01*grav*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
    2648      :               *(rp(il,i)-rr(il,i-1)))*dpinv
    2649 
    2650          fu(il,i)=fu(il,i)+0.01*grav*(mp(il,i+1)*(up(il,i+1)-u(il,i))
    2651      :             -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
    2652          fv(il,i)=fv(il,i)+0.01*grav*(mp(il,i+1)*(vp(il,i+1)-v(il,i))
    2653      :             -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
    2654         else  ! cvflag_grav
    2655          fr(il,i)=fr(il,i)+0.5*sigd*(evap(il,i)+evap(il,i+1))
    2656      :        +0.1*(mp(il,i+1)*(rp(il,i+1)-rr(il,i))-mp(il,i)
    2657      :               *(rp(il,i)-rr(il,i-1)))*dpinv
    2658          fu(il,i)=fu(il,i)+0.1*(mp(il,i+1)*(up(il,i+1)-u(il,i))
    2659      :             -mp(il,i)*(up(il,i)-u(il,i-1)))*dpinv
    2660          fv(il,i)=fv(il,i)+0.1*(mp(il,i+1)*(vp(il,i+1)-v(il,i))
    2661      :             -mp(il,i)*(vp(il,i)-v(il,i-1)))*dpinv
    2662         endif ! cvflag_grav
    2663 
    2664       endif ! i
    2665 1400  continue
     2874      do j=1,ntra
     2875       do k=i,nl+1
     2876        do il=1,ncum
     2877         if (i.le.inb(il) .and. k.le.inb(il)
     2878     $                .and. iflag(il) .le. 1) then
     2879          dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2880          cpinv=1.0/cpn(il,i)
     2881          if (cvflag_grav) then
     2882           ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv*ment(il,k,i)
     2883     :         *(traent(il,k,i,j)-tra(il,i,j))
     2884          else
     2885           ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv*ment(il,k,i)
     2886     :             *(traent(il,k,i,j)-tra(il,i,j))
     2887          endif
     2888         endif ! i and k
     2889        enddo
     2890       enddo
     2891      enddo
    26662892
    26672893c sb: interface with the cloud parameterization:          ! cld
    26682894
    26692895      do k=i+1,nl
    2670        do il=1,ncum
    2671         if (k.le.inb(il) .and. i.le.inb(il)) then         ! cld
     2896       do il=1,ncum
     2897        if (k.le.inb(il) .and. i.le.inb(il)
     2898     $               .and. iflag(il) .le. 1) then         ! cld
    26722899C (saturated downdrafts resulting from mixing)            ! cld
    26732900          qcond(il,i)=qcond(il,i)+elij(il,k,i)            ! cld
    … …  
    26792906C (particular case: no detraining level is found)         ! cld
    26802907      do il=1,ncum                                        ! cld
    2681        if (i.le.inb(il) .and. nent(il,i).eq.0) then       ! cld
     2908       if (i.le.inb(il) .and. nent(il,i).eq.0
     2909     $                 .and. iflag(il) .le. 1) then       ! cld
    26822910          qcond(il,i)=qcond(il,i)+(1.-ep(il,i))*clw(il,i) ! cld
    26832911          nqcond(il,i)=nqcond(il,i)+1.                    ! cld
    … …  
    26862914
    26872915      do il=1,ncum                                        ! cld
    2688        if (i.le.inb(il) .and. nqcond(il,i).ne.0.) then    ! cld
     2916       if (i.le.inb(il) .and. nqcond(il,i).ne.0
     2917     $                   .and. iflag(il) .le. 1) then     ! cld
    26892918          qcond(il,i)=qcond(il,i)/nqcond(il,i)            ! cld
    26902919       endif                                              ! cld
    26912920      enddo
    26922921
    2693 c      do j=1,ntra
    2694 c       do il=1,ncum
    2695 c        if (i.le.inb(il)) then
    2696 c         dpinv=1.0/(ph(il,i)-ph(il,i+1))
    2697 c         cpinv=1.0/cpn(il,i)
    2698 
    2699 c         if (cvflag_grav) then
    2700 c          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
    2701 c     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    2702 c     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
    2703 c         else
    2704 c          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
    2705 c     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
    2706 c     :     -mp(il,i)*(trap(il,i,j)-tra(il,i-1,j)))
    2707 c         endif
    2708 c        endif ! i
    2709 c       enddo
    2710 c      enddo
     2922      do j=1,ntra
     2923       do il=1,ncum
     2924        if (i.le.inb(il) .and. iflag(il) .le. 1) then
     2925         dpinv=1.0/(ph(il,i)-ph(il,i+1))
     2926         cpinv=1.0/cpn(il,i)
     2927
     2928         if (cvflag_grav) then
     2929          ftra(il,i,j)=ftra(il,i,j)+0.01*grav*dpinv
     2930     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     2931     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     2932         else
     2933          ftra(il,i,j)=ftra(il,i,j)+0.1*dpinv
     2934     :     *(mp(il,i+1)*(trap(il,i+1,j)-tra(il,i,j))
     2935     :     -mp(il,i)*(trap(il,i,j)-trap(il,i-1,j)))
     2936         endif
     2937        endif ! i
     2938       enddo
     2939      enddo
     2940
    27112941
    27122942500   continue
    … …  
    27182948c
    27192949      do 503 il=1,ncum
    2720 
     2950      IF (iflag(il) .le. 1) THEN
    27212951      ax=0.1*ment(il,inb(il),inb(il))*(hp(il,inb(il))-h(il,inb(il))
    27222952     : +t(il,inb(il))*(cpv-cpd)
    … …  
    27482978     :    +dx*(ph(il,inb(il))-ph(il,inb(il)+1))
    27492979     :       /(ph(il,inb(il)-1)-ph(il,inb(il)))
    2750 
     2980      ENDIF    !iflag
    27512981503   continue
    27522982
    2753 c      do j=1,ntra
    2754 c       do il=1,ncum
    2755 c        ex=0.1*ment(il,inb(il),inb(il))
    2756 c     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
    2757 c     :      /(ph(il,inb(il))-ph(il,inb(il)+1))
    2758 c        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
    2759 c        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
    2760 c     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
    2761 c     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
    2762 c       enddo
    2763 c      enddo
    2764 
    2765 c
    2766 c   ***    homoginize tendencies below cloud base    ***
     2983      do j=1,ntra
     2984       do il=1,ncum
     2985        IF (iflag(il) .le. 1) THEN
     2986        ex=0.1*ment(il,inb(il),inb(il))
     2987     :      *(traent(il,inb(il),inb(il),j)-tra(il,inb(il),j))
     2988     :      /(ph(i l,inb(il))-ph(il,inb(il)+1))
     2989        ftra(il,inb(il),j)=ftra(il,inb(il),j)-ex
     2990        ftra(il,inb(il)-1,j)=ftra(il,inb(il)-1,j)
     2991     :       +ex*(ph(il,inb(il))-ph(il,inb(il)+1))
     2992     :          /(ph(il,inb(il)-1)-ph(il,inb(il)))
     2993        ENDIF    !iflag
     2994       enddo
     2995      enddo
     2996
     2997c
     2998c   ***    homogenize tendencies below cloud base    ***
    27672999c
    27683000c
    … …  
    27763008      do i=1,nl
    27773009       do il=1,ncum
    2778         if (i.le.(icb(il)-1)) then
     3010        if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
    27793011      asum(il)=asum(il)+ft(il,i)*(ph(il,i)-ph(il,i+1))
    27803012      bsum(il)=bsum(il)+fr(il,i)*(lv(il,i)+(cl-cpd)*(t(il,i)-t(il,1)))
    … …  
    27833015     :                      *(ph(il,i)-ph(il,i+1))
    27843016      dsum(il)=dsum(il)+t(il,i)*(ph(il,i)-ph(il,i+1))/th(il,i)
    2785         endif 
     3017        endif
    27863018       enddo
    27873019      enddo
    … …  
    27903022      do i=1,nl
    27913023       do il=1,ncum
    2792         if (i.le.(icb(il)-1)) then
     3024        if (i.le.(icb(il)-1) .and. iflag(il) .le. 1) then
    27933025         ft(il,i)=asum(il)*t(il,i)/(th(il,i)*dsum(il))
    27943026         fr(il,i)=bsum(il)/csum(il)
    … …  
    28113043       enddo
    28123044      enddo
    2813      
     3045
    28143046      do i=1,nl
    28153047       do il=1,ncum
    … …  
    28553087      enddo
    28563088
     3089      do i=1,nl
     3090       do k=1,nl
     3091        do il=1,ncum
     3092         if(i.ge.icb(il)) then
     3093          if(k.ge.i.and. k.le.(inb(il))) then
     3094            upwd(il,i)=upwd(il,i)+m(il,k)
     3095          endif
     3096         else
     3097          if(k.lt.i) then
     3098            upwd(il,i)=upwd(il,i)+cbmf(il)*wghti(il,k)
     3099          endif
     3100        endif
     3101cc        print *,'cbmf',il,i,k,cbmf(il),wghti(il,k)
     3102        end do
     3103       end do
     3104      end do
     3105
    28573106      do i=2,nl
    28583107       do k=i,nl
    … …  
    28603109ctest         if (i.ge.icb(il).and.i.le.inb(il).and.k.le.inb(il)) then
    28613110         if (i.le.inb(il).and.k.le.inb(il)) then
    2862             upwd(il,i)=upwd(il,i)+m(il,k)+up1(il,k,i)
     3111            upwd(il,i)=upwd(il,i)+up1(il,k,i)
    28633112            dnwd(il,i)=dnwd(il,i)+dn1(il,k,i)
    28643113         endif
     3114cc         print *,'upwd',il,i,k,inb(il),upwd(il,i),m(il,k),up1(il,k,i)
    28653115        enddo
    28663116       enddo
    … …  
    28703120c!!!      DO il=1,ncum
    28713121c!!!      do i=icb(il),inb(il)
    2872 c!!!     
     3122c!!!
    28733123c!!!      upwd(il,i)=0.0
    28743124c!!!      dnwd(il,i)=0.0
    … …  
    28893139cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    28903140c        determination de la variation de flux ascendant entre
    2891 c        deux niveau non dilue mike
     3141c        deux niveau non dilue mip
    28923142cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    28933143
    28943144      do i=1,nl
    28953145       do il=1,ncum
    2896         mike(il,i)=m(il,i)
     3146        mip(il,i)=m(il,i)
    28973147       enddo
    28983148      enddo
    … …  
    29003150      do i=nl+1,nd
    29013151       do il=1,ncum
    2902         mike(il,i)=0.
     3152        mip(il,i)=0.
    29033153       enddo
    29043154      enddo
    … …  
    29693219        do k=i+1,nl+1                                       ! cld
    29703220         do il=1,ncum                                       ! cld
    2971           if (i.le.inb(il) .and. k.le.(inb(il)+1)) then     ! cld
     3221          if (i.le.inb(il) .and. k.le.(inb(il)+1)
     3222     $                     .and. iflag(il) .le. 1) then     ! cld
    29723223            mac(il,i)=mac(il,i)+m(il,k)                     ! cld
    29733224          endif                                             ! cld
    … …  
    29803231         do il=1,ncum                                       ! cld
    29813232          if (i.ge.icb(il) .and. i.le.(inb(il)-1)           ! cld
    2982      :      .and. j.ge.icb(il) ) then                       ! cld
     3233     :      .and. j.ge.icb(il)
     3234     :      .and. iflag(il) .le. 1 ) then                   ! cld
    29833235           sax(il,i)=sax(il,i)+rrd*(tvp(il,j)-tv(il,j))     ! cld
    29843236     :        *(ph(il,j)-ph(il,j+1))/p(il,j)                ! cld
    … …  
    29913243        do il=1,ncum                                        ! cld
    29923244         if (i.ge.icb(il) .and. i.le.(inb(il)-1)            ! cld
    2993      :       .and. sax(il,i).gt.0.0 ) then                  ! cld
     3245     :       .and. sax(il,i).gt.0.0
     3246     :       .and. iflag(il) .le. 1 ) then                  ! cld
    29943247           wa(il,i)=sqrt(2.*sax(il,i))                      ! cld
    29953248         endif                                              ! cld
    29963249        enddo                                               ! cld
    29973250       enddo                                                ! cld
    2998            
     3251
    29993252       do i=1,nl                                            ! cld
    30003253        do il=1,ncum                                        ! cld
    3001          if (wa(il,i).gt.0.0)                               ! cld
     3254         if (wa(il,i).gt.0.0 .and. iflag(il) .le. 1)        ! cld
    30023255     :     siga(il,i)=mac(il,i)/wa(il,i)                    ! cld
    30033256     :         *rrd*tvp(il,i)/p(il,i)/100./delta            ! cld
    … …  
    30083261     :           + (1.-siga(il,i))*qcond(il,i)              ! cld
    30093262         else if (iflag_clw.eq.1) then
    3010           qcondc(il,i)=qcond(il,i)              ! cld
     3263          qcondc(il,i)=qcond(il,i)                          ! cld
    30113264         endif
    30123265
    30133266        enddo                                               ! cld
    3014        enddo                                                ! cld
    3015 
     3267       enddo
     3268c        print*,'cv3_yield fin'       
     3269                                              ! cld
    30163270        return
    30173271        end
    30183272
    3019       SUBROUTINE cv3_tracer(nloc,len,ncum,nd,na,
    3020      &                        ment,sij,da,phi)
    3021         implicit none
    3022 c inputs:
    3023         integer ncum, nd, na, nloc,len
    3024         real ment(nloc,na,na),sij(nloc,na,na)
    3025 c ouputs:
    3026         real da(nloc,na),phi(nloc,na,na)
    3027 c local variables:
    3028         integer i,j,k
    3029 c       
    3030         da(:,:)=0.
    3031 c
    3032         do j=1,na
    3033           do k=1,na
    3034             do i=1,ncum
    3035             da(i,j)=da(i,j)+(1.-sij(i,k,j))*ment(i,k,j)
    3036             phi(i,j,k)=sij(i,k,j)*ment(i,k,j)
    3037 c            print *,'da',j,k,da(i,j),sij(i,k,j),ment(i,k,j)
    3038             end do
    3039           end do
    3040         end do
    3041    
    3042         return
    3043         end
    3044 
    30453273
    30463274      SUBROUTINE cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
    30473275     :         ,iflag
    3048      :         ,precip,VPrecip,sig,w0
     3276     :         ,precip,sig,w0
    30493277     :         ,ft,fq,fu,fv,ftra
    3050      :         ,inb
    30513278     :         ,Ma,upwd,dnwd,dnwd0,qcondc,wd,cape
    3052      :         ,da,phi,mp
    30533279     :         ,iflag1
    3054      :         ,precip1,VPrecip1,sig1,w01
     3280     :         ,precip1,sig1,w01
    30553281     :         ,ft1,fq1,fu1,fv1,ftra1
    3056      :         ,inb1
    30573282     :         ,Ma1,upwd1,dnwd1,dnwd01,qcondc1,wd1,cape1
    3058      :         ,da1,phi1,mp1)
     3283     :                               )
    30593284      implicit none
    30603285
    3061 #include "cvparam3.h"
     3286#include "cv3param.h"
    30623287
    30633288c inputs:
    … …  
    30653290      integer idcum(nloc)
    30663291      integer iflag(nloc)
    3067       integer inb(nloc)
    30683292      real precip(nloc)
    3069       real VPrecip(nloc,nd+1)
    30703293      real sig(nloc,nd), w0(nloc,nd)
    30713294      real ft(nloc,nd), fq(nloc,nd), fu(nloc,nd), fv(nloc,nd)
    … …  
    30753298      real qcondc(nloc,nd)
    30763299      real wd(nloc),cape(nloc)
    3077       real da(nloc,nd),phi(nloc,nd,nd),mp(nloc,nd)
    30783300
    30793301c outputs:
    30803302      integer iflag1(len)
    3081       integer inb1(len)
    30823303      real precip1(len)
    3083       real VPrecip1(len,nd+1)
    30843304      real sig1(len,nd), w01(len,nd)
    30853305      real ft1(len,nd), fq1(len,nd), fu1(len,nd), fv1(len,nd)
    … …  
    30893309      real qcondc1(nloc,nd)
    30903310      real wd1(nloc),cape1(nloc)
    3091       real da1(nloc,nd),phi1(nloc,nd,nd),mp1(nloc,nd)
    30923311
    30933312c local variables:
    … …  
    30983317         iflag1(idcum(i))=iflag(i)
    30993318         wd1(idcum(i))=wd(i)
    3100          inb1(idcum(i))=inb(i)
    31013319         cape1(idcum(i))=cape(i)
    31023320 2000   continue
    … …  
    31043322        do 2020 k=1,nl
    31053323          do 2010 i=1,ncum
    3106             VPrecip1(idcum(i),k)=VPrecip(i,k)
    31073324            sig1(idcum(i),k)=sig(i,k)
    31083325            w01(idcum(i),k)=w0(i,k)
    … …  
    31163333            dnwd01(idcum(i),k)=dnwd0(i,k)
    31173334            qcondc1(idcum(i),k)=qcondc(i,k)
    3118             da1(idcum(i),k)=da(i,k)
    3119             mp1(idcum(i),k)=mp(i,k)
    31203335 2010     continue
    31213336 2020   continue
    … …  
    31263341
    31273342
    3128 c        do 2100 j=1,ntra
    3129 c         do 2110 k=1,nd ! oct3
    3130 c          do 2120 i=1,ncum
    3131 c            ftra1(idcum(i),k,j)=ftra(i,k,j)
    3132 c 2120     continue
    3133 c 2110    continue
    3134 c 2100   continue
    3135         do j=1,nd
    3136          do k=1,nd
    3137           do i=1,ncum
    3138             phi1(idcum(i),k,j)=phi(i,k,j)
    3139           end do
    3140          end do
    3141         end do
    3142 
     3343        do 2100 j=1,ntra
     3344c oct3         do 2110 k=1,nl
     3345         do 2110 k=1,nd ! oct3
     3346          do 2120 i=1,ncum
     3347            ftra1(idcum(i),k,j)=ftra(i,k,j)
     3348 2120     continue
     3349 2110    continue
     3350 2100   continue
    31433351        return
    31443352        end

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

    r766 r879  
    312312c       (common cvparam)
    313313
    314       if (iflag_con.eq.3) then
    315        CALL cv3_param(nd,delt)
     314
     315      if (iflag_con.eq.30) then
     316       CALL cv30_param(nd,delt)
    316317      endif
    317318
    … …  
    361362 60   continue
    362363
    363       if (iflag_con.eq.3) then
     364      if (iflag_con.eq.30) then
    364365        do il=1,len
    365366         sig1(il,nd)=sig1(il,nd)+1.
    … …  
    372373!--------------------------------------------------------------------
    373374
    374       if (iflag_con.eq.3) then
    375        CALL cv3_prelim(len,nd,ndp1,t1,q1,p1,ph1            ! nd->na
     375      if (iflag_con.eq.30) then
     376
     377       print*,'Emanuel version 30 '
     378       CALL cv30_prelim(len,nd,ndp1,t1,q1,p1,ph1            ! nd->na
    376379     o               ,lv1,cpn1,tv1,gz1,h1,hm1,th1)
    377380      endif
    … …  
    386389!--------------------------------------------------------------------
    387390
    388       if (iflag_con.eq.3) then
    389        CALL cv3_feed(len,nd,t1,q1,qs1,p1,ph1,hm1,gz1           ! nd->na
     391      if (iflag_con.eq.30) then
     392       CALL cv30_feed(len,nd,t1,q1,qs1,p1,ph1,hm1,gz1           ! nd->na
    390393     o         ,nk1,icb1,icbmax,iflag1,tnk1,qnk1,gznk1,plcl1)
    391394      endif
    … …  
    403406!--------------------------------------------------------------------
    404407
    405       if (iflag_con.eq.3) then
    406        CALL cv3_undilute1(len,nd,t1,q1,qs1,gz1,plcl1,p1,nk1,icb1  ! nd->na
     408      if (iflag_con.eq.30) then
     409       CALL cv30_undilute1(len,nd,t1,q1,qs1,gz1,plcl1,p1,nk1,icb1  ! nd->na
    407410     o                        ,tp1,tvp1,clw1,icbs1)
    408411      endif
    … …  
    417420!-------------------------------------------------------------------
    418421
    419       if (iflag_con.eq.3) then
    420        CALL cv3_trigger(len,nd,icb1,plcl1,p1,th1,tv1,tvp1      ! nd->na
     422      if (iflag_con.eq.30) then
     423       CALL cv30_trigger(len,nd,icb1,plcl1,p1,th1,tv1,tvp1      ! nd->na
    421424     o                 ,pbase1,buoybase1,iflag1,sig1,w01)
    422425      endif
    … …  
    447450!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    448451
    449       if (iflag_con.eq.3) then
    450        CALL cv3_compress( len,nloc,ncum,nd,ntra
     452      if (iflag_con.eq.30) then
     453       CALL cv30_compress( len,nloc,ncum,nd,ntra
    451454     :    ,iflag1,nk1,icb1,icbs1
    452455     :    ,plcl1,tnk1,qnk1,gznk1,pbase1,buoybase1
    … …  
    485488!-------------------------------------------------------------------
    486489
    487       if (iflag_con.eq.3) then
    488        CALL cv3_undilute2(nloc,ncum,nd,icb,icbs,nk        !na->nd
     490      if (iflag_con.eq.30) then
     491       CALL cv30_undilute2(nloc,ncum,nd,icb,icbs,nk        !na->nd
    489492     :                        ,tnk,qnk,gznk,t,q,qs,gz
    490493     :                        ,p,h,tv,lv,pbase,buoybase,plcl
    … …  
    503506!-------------------------------------------------------------------
    504507
    505       if (iflag_con.eq.3) then
    506        CALL cv3_closure(nloc,ncum,nd,icb,inb              ! na->nd
     508      if (iflag_con.eq.30) then
     509       CALL cv30_closure(nloc,ncum,nd,icb,inb              ! na->nd
    507510     :                       ,pbase,p,ph,tv,buoy
    508511     o                       ,sig,w0,cape,m)
    … …  
    519522!-------------------------------------------------------------------
    520523
    521       if (iflag_con.eq.3) then
    522        CALL cv3_mixing(nloc,ncum,nd,nd,ntra,icb,nk,inb    ! na->nd
     524      if (iflag_con.eq.30) then
     525       CALL cv30_mixing(nloc,ncum,nd,nd,ntra,icb,nk,inb    ! na->nd
    523526     :                     ,ph,t,q,qs,u,v,tra,h,lv,qnk
    524527     :                     ,hp,tv,tvp,ep,clw,m,sig
    … …  
    537540!-------------------------------------------------------------------
    538541
    539       if (iflag_con.eq.3) then
    540        CALL cv3_unsat(nloc,ncum,nd,nd,ntra,icb,inb    ! na->nd
     542      if (iflag_con.eq.30) then
     543       CALL cv30_unsat(nloc,ncum,nd,nd,ntra,icb,inb    ! na->nd
    541544     :               ,t,q,qs,gz,u,v,tra,p,ph
    542545     :               ,th,tv,lv,cpn,ep,sigp,clw
    … …  
    557560!-------------------------------------------------------------------
    558561
    559       if (iflag_con.eq.3) then
    560        CALL cv3_yield(nloc,ncum,nd,nd,ntra            ! na->nd
     562      if (iflag_con.eq.30) then
     563       CALL cv30_yield(nloc,ncum,nd,nd,ntra            ! na->nd
    561564     :                     ,icb,inb,delt
    562565     :                     ,t,q,u,v,tra,gz,p,ph,h,hp,lv,cpn,th
    … …  
    584587!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    585588
    586       if (iflag_con.eq.3) then
    587        CALL cv3_tracer(nloc,len,ncum,nd,nd,
     589      if (iflag_con.eq.30) then
     590       CALL cv30_tracer(nloc,len,ncum,nd,nd,
    588591     :                  ment,sij,da,phi)
    589592      endif
    … …  
    597600      end do
    598601c
    599       if (iflag_con.eq.3) then
    600        CALL cv3_uncompress(nloc,len,ncum,nd,ntra,idcum
     602      if (iflag_con.eq.30) then
     603       CALL cv30_uncompress(nloc,len,ncum,nd,ntra,idcum
    601604     :          ,iflag
    602605     :          ,precip,VPrecip,sig,w0
    … …  
    670673       t0=273.15
    671674       grav=g
    672       endif
     675      else
    673676
    674677c constants consistent with LMDZ:
    675       if (iflag_con.eq.3) then
    676678       cpd = RCPD
    677679       cpv = RCPV

  • LMDZ4/trunk/libf/phylmd/diagphy.F

    r766 r879  
    4343C
    4444C J.L. Dufresne, July 2002
     45C Version prise sur ~rlmd833/LMDZOR_201102/modipsl/modeles/LMDZ.3.3/libf/phylmd
     46C  le 25 Novembre 2002.
    4547C======================================================================
    4648C

  • LMDZ4/trunk/libf/phylmd/ini_histday.h

    r878 r879  
    143143     .                "ave(X)", zstophy,zout)
    144144c
     145         CALL histdef(nid_day, "plul","Precipitation ls liq+sol"
     146     .                , "kg/(s*m2)",
     147     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     148     .                "ave(X)", zstophy,zout)
     149c
    145150         CALL histdef(nid_day, "snowf", "Snow fall", "kg/(m2*s)",
    146151     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    … …  
    333338     .                "ave(X)", zstophy,zout)
    334339c
     340         CALL histdef(nid_day, "rhum", "relative humidity", " ",
     341     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     342     .                "ave(X)", zstophy,zout)
     343c
    335344         CALL histdef(nid_day, "ovap", "Specific humidity", "kg/kg",
    336345     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    … …  
    599608     .                "ave(X)", zstophy,zout)
    600609c
     610         CALL histdef(nid_day, "dqpbl", "PBL dT", "K/s",
     611     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     612     .                "ave(X)", zstophy,zout)
     613c
    601614         CALL histdef(nid_day, "rnebcon", "Convective Cloud Fraction"
    602615     .                , "-",
     616     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     617     .                "ave(X)", zstophy,zout)
     618c
     619         CALL histdef(nid_day, "dtwak", "Wake dT", "K/s",
    603620     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    604621     .                "ave(X)", zstophy,zout)
    … …  
    621638     .                "ave(X)", zstophy,zout)
    622639c
     640
     641         CALL histdef(nid_day, "cape", "Conv avlbl pot ener", "J/kg",
     642     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     643     .                "ave(X)", zstophy,zout)
     644
     645         CALL histdef(nid_day, "wh", "Wake height", "m",
     646     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     647     .                "ave(X)", zstophy,zout)
     648c
     649         CALL histdef(nid_day, "ws", "Wake surface", "m2",
     650     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     651     .                "ave(X)", zstophy,zout)
     652c
    623653        CALL histdef(nid_day,"meantaucld",
    624654     .                "ISCCP mean cloud optical thickness","1",
    … …  
    634664     .                "ave(X)", zstophy,zout)
    635665c
     666c
     667      CALL histdef(nid_day, "cin", "Conv Inhibition", "J/kg",
     668     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     669     .                "ave(X)", zstophy,zout)     
     670c
     671         CALL histdef(nid_day, "wale", "Wake Available Energy", "J/kg",
     672     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     673     .                "ave(X)", zstophy,zout)
     674c
     675         CALL histdef(nid_day, "walp",
     676     .                "Available Lifting Energy due to wake", "W/m2",
     677     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     678     .                "ave(X)", zstophy,zout)
     679c
     680         CALL histdef(nid_day, "blale", "PBL Available Energy", "J/kg",
     681     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     682     .                "ave(X)", zstophy,zout)
     683c
     684         CALL histdef(nid_day, "blalp",
     685     .                "Available Lifting Energy due to PBL", "W/m2",
     686     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     687     .                "ave(X)", zstophy,zout)
     688c
     689         CALL histdef(nid_day, "wdt1",
     690     .                "Temp diff wake layer 1", "K",
     691     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     692     .                "ave(X)", zstophy,zout)
     693c
     694         CALL histdef(nid_day, "wdq1",
     695     .                "Temp diff wake layer 1", "g/kg",
     696     .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     697     .                "ave(X)", zstophy,zout)
    636698cIM rajout AMMA-MIP
    637699c

  • LMDZ4/trunk/libf/phylmd/ini_histmth.h

    r878 r879  
    3434     .                 klev, presnivs/100., nvert)
    3535c
     36         CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s",
     37     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     38     .                "ave(X)", zstophy,zout)
     39
     40         CALL histdef(nid_mth, "dtcon", "Convection dT", "K/s",
     41     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
     42     .                "ave(X)", zstophy,zout)
     43c
    3644      IF(type_run.EQ."CLIM".OR.type_run.EQ."ENSP") THEN
    3745c
    … …  
    7280c
    7381c ENSEMBLES BEG
    74          CALL histdef(nid_mth, "t2m_min", "Temp. 2m min.",
    75      .                "K",
    76      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    77      .                t2mincels, zstophy,zout)
    78 c
    79          CALL histdef(nid_mth, "t2m_max", "Temp. 2m max.",
    80      .                "K",
    81      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    82      .                t2maxcels, zstophy,zout)
     82c        CALL histdef(nid_mth, "t2m_min", "Temp. 2m min.",
     83c    .                "K",
     84c    .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     85c    .                t2mincels, zstophy,zout)
     86c
     87c        CALL histdef(nid_mth, "t2m_max", "Temp. 2m max.",
     88c    .                "K",
     89c    .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     90c    .                t2maxcels, zstophy,zout)
    8391c
    8492c        CALL histdef(nid_mth, "tsoil", "Sfce soil Temperature",
    … …  
    122130      endif
    123131c
    124          CALL histdef(nid_mth, "ndayrain",
    125      .                "Number of day with rain (liq+sol)", "-",
    126      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    127      .                "inst(X)", zstomth,zout)
     132c        CALL histdef(nid_mth, "ndayrain",
     133c    .                "Number of day with rain (liq+sol)", "-",
     134c    .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     135c    .                "inst(X)", zstomth,zout)
    128136c
    129137         CALL histdef(nid_mth, "precip", "Precipitation Totale liq+sol",
    … …  
    990998     .                "ave(X)", zstophy,zout)
    991999c
    992           CALL histdef(nid_mth, "dtcon","dtcon","K/s",
    993      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    994      .                "ave(X)", zstophy,zout)
    995 c
    996          CALL histdef(nid_mth, "dtthe", "Dry adjust. dT", "K/s",
    997      .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    998      .                "ave(X)", zstophy,zout)
    999 
    10001000         CALL histdef(nid_mth,"dqthe","Dry adjust. dQ","(kg/kg)/s",
    10011001     .                iim,jj_nb,nhori, klev,1,klev,nvert, 32,
    … …  
    10831083     .                "ave(X)", zstorad,zout)
    10841084c
    1085          CALL histdef(nid_mth, "t2m_min", "Temp. 2m min.",
    1086      .                "K",
    1087      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    1088      .                t2mincels, zstophy,zout)
    1089 c
    1090          CALL histdef(nid_mth, "t2m_max", "Temp. 2m max.",
    1091      .                "K",
    1092      .                iim,jj_nb,nhori, 1,1,1, -99, 32,
    1093      .                t2maxcels, zstophy,zout)
     1085c        CALL histdef(nid_mth, "t2m_min", "Temp. 2m min.",
     1086c    .                "K",
     1087c    .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     1088c    .                t2mincels, zstophy,zout)
     1089c
     1090c        CALL histdef(nid_mth, "t2m_max", "Temp. 2m max.",
     1091c    .                "K",
     1092c    .                iim,jj_nb,nhori, 1,1,1, -99, 32,
     1093c    .                t2maxcels, zstophy,zout)
    10941094c
    10951095c        CALL histdef(nid_mth, "tsoil", "Sfce soil Temperature",

  • LMDZ4/trunk/libf/phylmd/ini_paramLMDZ_phy.h

    r828 r879  
    1010c
    1111       CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlon,zx_lon)
     12       if (iim.gt.1) then
    1213       DO i = 1, iim
    1314         zx_lon(i,1) = rlon(i+1)
    1415         zx_lon(i,jjmp1) = rlon(i+1)
    1516       ENDDO
     17       endif
    1618       CALL gr_fi_ecrit(1,klon,iim,jjmp1,rlat,zx_lat)
    1719c

  • LMDZ4/trunk/libf/phylmd/iniphysiq.F

    r805 r879  
    8888      rlatd(1:klon_omp)=plat(offset+klon_omp_begin:offset+klon_omp_end)
    8989
    90       call suphec
     90      call suphel
    9191
    9292c$OMP END PARALLEL

  • LMDZ4/trunk/libf/phylmd/iophy.F90

    r793 r879  
     1!
     2! $Header$
     3!
    14module iophy
    25 
    … …  
    2225  include 'dimensions.h'   
    2326    real,dimension(iim),intent(in) :: lon
    24     real,dimension(jjm+1),intent(in) :: lat
     27    real,dimension(jjm+1-1/iim),intent(in) :: lat
    2528
    2629    INTEGER,DIMENSION(2) :: ddid
    … …  
    3336
    3437!$OMP MASTER 
    35     allocate(io_lat(jjm+1))
     38    allocate(io_lat(jjm+1-1/iim))
    3639    io_lat(:)=lat(:)
    3740    allocate(io_lon(iim))
    … …  
    4548   
    4649    ddid=(/ 1,2 /)
    47     dsg=(/ iim, jjm+1 /)
     50    dsg=(/ iim, jjm+1-1/iim /)
    4851    dsl=(/ iim, jj_nb /)
    4952    dpf=(/ 1,jj_begin /)
    … …  
    167170!         enddo
    168171!         
    169 !         if (jjphy_end==jjm+1) then
     172!         if (jjphy_end==jjm+1-1/iim) then
    170173!             field_dyn(:,jjphy_nb,l)=field_phy(klon_mpi,l)
    171174!         else

  • LMDZ4/trunk/libf/phylmd/mod_phys_lmdz_mpi_data.F90

    r775 r879  
    6464#endif
    6565   
    66     klon_glo=iim*(jjp1-2)+2
     66    if (iim.eq.1) then
     67       klon_glo=1
     68    else
     69       klon_glo=iim*(jjp1-2)+2
     70    endif
    6771   
    6872    COMM_LMDZ_PHY=COMM_LMDZ

  • LMDZ4/trunk/libf/phylmd/phyetat0.F

    r878 r879  
    138138      CHARACTER*7 str7
    139139      CHARACTER*2 str2
    140       real iolat(jjm+1)
     140
     141c FH1D
     142c     real iolat(jjm+1)
     143      real iolat(jjm+1-1/iim)
    141144c
    142145c Ouvrir le fichier contenant l'etat initial:
    … …  
    15491552cym  en attendant mieux
    15501553        iolat(1)=rlat(1)
    1551         iolat(jjm+1)=rlat(klon_glo)
     1554
     1555!FH1D   
     1556!iolat(jjm+1)=rlat(klon_glo)
     1557        iolat(jjm+1-1/iim)=rlat(klon_glo)
     1558        if (iim.gt.1) then
    15521559        do i=2,jjm
    15531560          iolat(i)=rlat(2+(i-2)*iim)
    15541561        enddo
     1562        endif
     1563
    15551564        CALL bcast_mpi(iolat)
    15561565        CALL bcast_mpi(rlon)
    1557         call init_iophy(iolat,rlon(2:iim+1))
     1566
     1567!FH1D
     1568!       call init_iophy(iolat,rlon(2:iim+1))
     1569        call init_iophy(iolat,rlon(2-1/iim:iim+1-1/iim))
    15581570       
    15591571c$OMP END MASTER

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

    r878 r879  
    108108      LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
    109109      PARAMETER (ok_gust=.FALSE.)
     110      integer iflag_radia     ! active ou non le rayonnement (MPL)
     111      save iflag_radia
    110112c======================================================================
    111113      LOGICAL check ! Verifier la conservation du modele en eau
    … …  
    131133      REAL fluxg(klon)    !flux turbulents ocean-atmosphere
    132134      REAL amn, amx
     135      INTEGER igout
    133136c======================================================================
    134137c Clef controlant l'activation du cycle diurne:
    … …  
    222225      REAL u(klon,klev)
    223226      REAL v(klon,klev)
    224       REAL t(klon,klev)
     227      REAL t(klon,klev),theta(klon,klev)
    225228      REAL qx(klon,klev,nqmax)
    226229
    … …  
    784787cym      SAVE wd       ! sb
    785788
     789      REAL,allocatable,save :: sigd(:)
     790c
     791c=================================================================================================
     792cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
     793c Variables liées à la poche froide (jyg)
     794
     795      REAL mip(klon,klev)  ! mass flux shed by the adiab ascent at each level
     796      REAL Vprecip(klon,klev)   ! precipitation vertical profile
     797      REAL,allocatable,save :: cin(:)            ! CIN
     798      REAL,allocatable,save :: ftd(:,:)       ! differential heating between wake and environment
     799      REAL,allocatable,save :: fqd(:,:)       ! differential moistening between wake and environment
     800c34EK
     801c -- Variables de controle de ALE et ALP
     802      REAL,allocatable,save :: ALE(:)     ! Energie disponible pour soulevement : utilisee par la         
     803c convection d'Emanuel pour le declenchement et la regulation
     804      REAL,allocatable,save :: ALP(:)     ! Puissance  disponible pour soulevement           !
     805c
     806      REAL wape_prescr, fip_prescr
     807      INTEGER it_wape_prescr
     808      SAVE wape_prescr, fip_prescr, it_wape_prescr
     809c
     810c variables supplementaires de concvl
     811      REAL Tconv(klon,klev)
     812      REAL ment(klon,klev,klev),sij(klon,klev,klev)
     813      REAL dd_t(klon,klev),dd_q(klon,klev)
     814cnouvelles variables pour le couplage convection-couche limite
     815      real,allocatable,save :: Ale_bl(:)
     816      real,allocatable,save :: Alp_bl(:)
     817      integer,allocatable,save :: lalim_conv(:)
     818      real,allocatable,save :: wght_th(:,:)
     819      real alp_bl_prescr
     820      save alp_bl_prescr
     821      real ale_bl_prescr
     822      save ale_bl_prescr
     823      real ale_wake(klon)
     824      real alp_wake(klon)
     825cRC
     826c Variables liées à la poche froide (jyg et rr)
     827c Version diagnostique pour l'instant : pas de rétroaction sur la convection
     828
     829      REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
     830
     831      REAL,allocatable,save :: wake_deltat(:,:)     ! Wake : ecart de temperature avec la
     832c                                              zone non perturbee
     833      REAL,allocatable,save :: wake_deltaq(:,:)     ! Wake : ecart d'humidite avec la
     834c                                              zone non perturbee
     835      REAL wake_dth(klon,klev)        ! wake : temp pot difference
     836
     837      REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
     838      REAL wake_omgbdth(klon,klev)    ! Wake : flux of Delta_Theta transported by LS omega
     839      REAL wake_dp_omgb(klon,klev)    ! Wake : vertical gradient of large scale omega
     840      REAL wake_dtKE(klon,klev)       ! Wake : differential heating (wake - unpertubed) CONV
     841      REAL wake_dqKE(klon,klev)       ! Wake : differential moistening (wake - unpertubed) CONV
     842      REAL wake_dtPBL(klon,klev)      ! Wake : differential heating (wake - unpertubed) PBL
     843      REAL wake_dqPBL(klon,klev)      ! Wake : differential moistening (wake - unpertubed) PBL
     844      REAL wake_omg(klon,klev)        ! Wake : velocity difference (wake - unpertubed)
     845      REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
     846      REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
     847      REAL wake_spread(klon,klev)     ! spreading term in wake_delt
     848      REAL,allocatable,save :: wake_Cstar(:)           ! vitesse d'etalement de la poche
     849cpourquoi y'a pas de save??
     850      REAL wake_h(klon)               ! Wake : hauteur de la poche froide
     851      REAL,allocatable,save :: wake_s(:)               ! Wake : fraction surfacique occupee par
     852c                                              la poche froide
     853      INTEGER wake_k(klon)            ! Wake sommet
     854c
     855      REAL t_undi(klon,klev)               ! temperature moyenne dans la zone non perturbee
     856      REAL q_undi(klon,klev)               ! humidite moyenne dans la zone non perturbee
     857c
     858      REAL wake_pe(klon)              ! Wake potential energy - WAPE
     859      REAL,allocatable,save :: wake_fip(:)             ! Gust Front Impinging power - ALP
     860
     861      REAL wake_gfl(klon)             ! Gust Front Length
     862      REAL wake_dens(klon)
     863c
     864      REAL,allocatable,save :: dt_wake(:,:)
     865      REAL,allocatable,save :: dq_wake(:,:) ! LS tendencies due to wake
     866c
     867      REAL dt_dwn(klon,klev)
     868      REAL dq_dwn(klon,klev)
     869      REAL wdt_PBL(klon,klev)
     870      REAL udt_PBL(klon,klev)
     871      REAL wdq_PBL(klon,klev)
     872      REAL udq_PBL(klon,klev)
     873      REAL M_dwn(klon,klev)
     874      REAL M_up(klon,klev)
     875      REAL dt_a(klon,klev)
     876      REAL dq_a(klon,klev)
     877c
     878cRR:fin declarations poches froides
     879c=======================================================================================================
     880
    786881c Variables locales pour la couche limite (al1):
    787882c
    … …  
    9091004      EXTERNAL phyredem  ! ecrire l'etat de redemarrage de la physique
    9101005      EXTERNAL radlwsw   ! rayonnements solaire et infrarouge
    911       EXTERNAL suphec    ! initialiser certaines constantes
     1006      EXTERNAL suphel    ! initialiser certaines constantes
    9121007      EXTERNAL transp    ! transport total de l'eau et de l'energie
    9131008      EXTERNAL ecribina  ! ecrire le fichier binaire global
    … …  
    13781473c
    13791474c======================================================================
     1475! Ecriture eventuelle d'un profil verticale en entree de la physique.
     1476! Utilise notamment en 1D mais peut etre active egalement en 3D
     1477! en imposant la valeur de igout.
     1478c======================================================================
     1479
     1480      if (klon.eq.1) then
     1481          print*,'WARNING !!!! omega=0'
     1482          omega=0.
     1483          igout=1
     1484         write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
     1485         write(lunout,*)
     1486     s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
     1487         write(lunout,*)
     1488     s  nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
     1489
     1490         write(lunout,*) 'papers, play, phi, u, v, t, omega'
     1491         do k=1,nlev
     1492            write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
     1493     s   u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
     1494         enddo
     1495         write(lunout,*) 'ovap (g/kg),  oliq (g/kg)'
     1496         do k=1,nlev
     1497            write(lunout,*) qx(igout,k,1)*1000,qx(igout,k,2)*1000.
     1498         enddo
     1499      endif
     1500
     1501c======================================================================
    13801502
    13811503cym => necessaire pour iflag_con != 2   
    … …  
    14171539      allocate( zuthe(klon),zvthe(klon))
    14181540      allocate( alb_neig(klon))
    1419       allocate( ema_workcbmf(klon))   
     1541      allocate( ema_workcbmf(klon))
     1542cCR:nvelles variables convection/poches froides
     1543      allocate( sigd(klon))
     1544      allocate( cin(klon))
     1545      allocate( ftd(klon,klev))
     1546      allocate( fqd(klon,klev))
     1547      allocate( ALE(klon))
     1548      allocate( ALP(klon))
     1549      allocate( Ale_bl(klon))
     1550      allocate( Alp_bl(klon))
     1551      allocate( lalim_conv(klon))
     1552      allocate( wght_th(klon,klev))
     1553      allocate( wake_deltat(klon,klev))
     1554      allocate( wake_deltaq(klon,klev))
     1555      allocate( wake_Cstar(klon))
     1556      allocate( wake_s(klon))
     1557      allocate( wake_fip(klon))
     1558      allocate( dt_wake(klon,klev))
     1559      allocate( dq_wake(klon,klev))
     1560cfinCR
    14201561      allocate( ema_cbmf(klon))
    14211562      allocate( ema_pcb(klon))
    … …  
    15281669      ENDIF
    15291670      IF (debut) THEN
    1530          CALL suphec ! initialiser constantes et parametres phys.
     1671         CALL suphel ! initialiser constantes et parametres phys.
    15311672      ENDIF
    15321673
    … …  
    15431684C
    15441685!rv
     1686cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
     1687cde la convection a partir des caracteristiques du thermique
     1688         wght_th(:,:)=1.
     1689         lalim_conv(:)=1
     1690cRC
    15451691         u10m(:,:)=0.
    15461692         v10m(:,:)=0.
    … …  
    15931739         call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
    15941740     .                  ok_instan, ok_hf, seuil_inversion,
    1595      .                  fact_cldcon, facttemps,ok_newmicro,
     1741     .                  fact_cldcon, facttemps,ok_newmicro,iflag_radia,
    15961742     .                  iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
    15971743     .                  ok_ade, ok_aie,
    15981744     .                  bl95_b0, bl95_b1,
    1599      .                  iflag_thermals,nsplit_thermals)
     1745     .                  iflag_thermals,nsplit_thermals,
     1746cnv flags pour la convection et les poches froides
     1747     .                   iflag_coupl,iflag_clos,iflag_wake)
     1748
     1749      print*,'iflag_coupl,iflag_clos,iflag_wake',
     1750     .   iflag_coupl,iflag_clos,iflag_wake
    16001751
    16011752c
    … …  
    17111862          ENDDO
    17121863cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
     1864c===============================================================================
     1865cCR:04.12.07: initialisations poches froides
     1866c Controle de ALE et ALP pour la fermeture convective (jyg)
     1867          if (iflag_wake.eq.1) then
     1868            CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
     1869     s                  ,alp_bl_prescr, ale_bl_prescr)
     1870c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
     1871c        print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
     1872          endif
     1873
     1874        do i = 1,klon
     1875         wake_s(i) = 0.
     1876         wake_fip(i) = 0.
     1877         wake_cstar(i) = 0.
     1878         DO k=1,klev
     1879          wake_deltat(i,k)=0.
     1880          wake_deltaq(i,k)=0.
     1881         ENDDO
     1882        enddo
     1883c================================================================================
    17131884
    17141885         ENDIF
    … …  
    22302401c supprimer les calculs / ftra.
    22312402              ntra = 1
     2403
     2404c=====================================================================================
     2405cajout pour la parametrisation des poches froides:
     2406ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
     2407      do k=1,klev
     2408            do i=1,klon
     2409             if (iflag_wake.eq.1) then
     2410             t_wake(i,k) = t_seri(i,k)
     2411     .           +(1-wake_s(i))*wake_deltat(i,k)
     2412             q_wake(i,k) = q_seri(i,k)
     2413     .           +(1-wake_s(i))*wake_deltaq(i,k)
     2414             t_undi(i,k) = t_seri(i,k)
     2415     .           -wake_s(i)*wake_deltat(i,k)
     2416             q_undi(i,k) = q_seri(i,k)
     2417     .           -wake_s(i)*wake_deltaq(i,k)
     2418             else
     2419             t_wake(i,k) = t_seri(i,k)
     2420             q_wake(i,k) = q_seri(i,k)
     2421             t_undi(i,k) = t_seri(i,k)
     2422             q_undi(i,k) = q_seri(i,k)
     2423             endif
     2424            enddo
     2425         enddo
     2426     
     2427cc--   Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
     2428cc--    pour le soulevement des particules dans le modele convectif
     2429c
     2430      do i = 1,klon
     2431        ALE(i) = 0.
     2432        ALP(i) = 0.
     2433      enddo
     2434c
     2435ccalcul de ale_wake et alp_wake
     2436       do i = 1,klon
     2437          if (iflag_wake.eq.1) then
     2438          ale_wake(i) = 0.5*wake_cstar(i)**2
     2439          alp_wake(i) = wake_fip(i)
     2440          else
     2441          ale_wake(i) = 0.
     2442          alp_wake(i) = 0.
     2443          endif
     2444       enddo
     2445ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
     2446cdans le thermique sinon
     2447       if (iflag_coupl.eq.0) then
     2448          print*,'ALE et ALP imposes'
     2449          do i = 1,klon
     2450con ne couple que ale
     2451c           ALE(i) = max(ale_wake(i),Ale_bl(i))
     2452            ALE(i) = max(ale_wake(i),ale_bl_prescr)
     2453con ne couple que alp
     2454c           ALP(i) = alp_wake(i) + Alp_bl(i)
     2455            ALP(i) = alp_wake(i) + alp_bl_prescr
     2456          enddo
     2457       else
     2458          print*,'ALE et ALP couples au thermique'
     2459          do i = 1,klon
     2460              ALE(i) = max(ale_wake(i),Ale_bl(i))
     2461              ALP(i) = alp_wake(i) + Alp_bl(i)
     2462c         write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
     2463c         write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
     2464          enddo
     2465       endif
     2466
     2467cfin calcul ale et alp
     2468c=================================================================================================
     2469
     2470
    22322471c sb, oct02:
    22332472c Schema de convection modularise et vectorise:
    … …  
    22362475          IF (ok_cvl) THEN ! new driver for convectL
    22372476
    2238           CALL concvl (iflag_con,
    2239      .        dtime,paprs,pplay,t_seri,q_seri,
    2240      .        u_seri,v_seri,tr_seri,ntra,
     2477          CALL concvl (iflag_con,iflag_clos,
     2478     .        dtime,paprs,pplay,t_undi,q_undi,
     2479     .        t_wake,q_wake,
     2480     .        u_seri,v_seri,tr_seri,nbtr,
     2481     .        ALE,ALP,
    22412482     .        ema_work1,ema_work2,
    22422483     .        d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
    2243      .        rain_con, snow_con, ibas_con, itop_con,
     2484     .        rain_con, snow_con, ibas_con, itop_con, sigd,
    22442485     .        upwd,dnwd,dnwd0,
    2245      .        Ma,cape,tvp,iflagctrl,
     2486     .        Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
    22462487     .        pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
    2247      .        pmflxr,pmflxs,
    2248      .        da,phi,mp)
     2488     .        pmflxr,pmflxs,da,phi,mp,
     2489     .        ftd,fqd,lalim_conv,wght_th)
    22492490
    22502491cIM cf. FH
    … …  
    23142555          ENDDO
    23152556          DO i = 1, klon
    2316             ema_pct(i)  = paprs(i,itop_con(i))
     2557
     2558! L'idicage de itop_con peut cacher un pb potentiel
     2559! FH sous la dictee de JYG, CR
     2560            ema_pct(i)  = paprs(i,itop_con(i)+1)
     2561
    23172562            if (itop_con(i).gt.klev-3) then
    23182563               print*,'La convection monte trop haut '
    … …  
    23972642      ENDIF
    23982643      zx_ajustq=.FALSE.
     2644
     2645c
     2646c=============================================================================
     2647cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
     2648cpour la couche limite diffuse pour l instant
     2649c
     2650      if (iflag_wake.eq.1) then
     2651      DO k=1,klev
     2652        DO i=1,klon
     2653          dt_dwn(i,k)  = ftd(i,k)
     2654          wdt_PBL(i,k) = 0.
     2655          dq_dwn(i,k)  = fqd(i,k)
     2656          wdq_PBL(i,k) = 0.
     2657          M_dwn(i,k)   = dnwd0(i,k)
     2658          M_up(i,k)    = upwd(i,k)
     2659          dt_a(i,k)    = d_t_con(i,k)/dtime - ftd(i,k)
     2660          udt_PBL(i,k) = 0.
     2661          dq_a(i,k)    = d_q_con(i,k)/dtime - fqd(i,k)
     2662          udq_PBL(i,k) = 0.
     2663        ENDDO
     2664      ENDDO
     2665c
     2666ccalcul caracteristiques de la poche froide
     2667      call calWAKE (paprs,pplay,dtime
     2668     :               ,t_seri,q_seri,omega,ibas_con
     2669     :               ,dt_dwn,dq_dwn,M_dwn,M_up
     2670     :               ,dt_a,dq_a,sigd
     2671     :               ,wdt_PBL,wdq_PBL
     2672     :               ,udt_PBL,udq_PBL
     2673     o               ,wake_deltat,wake_deltaq,wake_dth
     2674     o               ,wake_h,wake_s,wake_dens
     2675     o               ,wake_pe,wake_fip,wake_gfl
     2676     o               ,dt_wake,dq_wake
     2677     o               ,wake_k, t_undi,q_undi
     2678     o               ,wake_omgbdth,wake_dp_omgb
     2679     o               ,wake_dtKE,wake_dqKE
     2680     o               ,wake_dtPBL,wake_dqPBL
     2681     o               ,wake_omg,wake_dp_deltomg
     2682     o               ,wake_spread,wake_Cstar,wake_d_deltat_gw
     2683     o               ,wake_ddeltat,wake_ddeltaq)
     2684c
     2685cIncrementation des tendances de la poche froide
     2686      DO k = 1, klev
     2687        DO i = 1, klon
     2688          t_seri(i,k) = t_seri(i,k) + dt_wake(i,k)*dtime
     2689          q_seri(i,k) = q_seri(i,k) + dq_wake(i,k)*dtime
     2690        ENDDO
     2691      ENDDO
     2692
     2693      endif
     2694c      print*,'apres callwake iflag_cldcon=', iflag_cldcon
    23992695c
    24002696c===================================================================
    … …  
    24412737     s      ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
    24422738     s      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
    2443      s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth,
    2444      s       ratqsdiff,zqsatth)
     2739     s      ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
     2740     s      ,ratqsdiff,zqsatth
     2741con rajoute ale et alp, et les caracteristiques de la couche alim
     2742     s      ,Ale_bl,Alp_bl,lalim_conv,wght_th)
    24452743         endif
    24462744
    … …  
    30163314      ENDIF
    30173315      itaprad = itaprad + 1
     3316
     3317      if (iflag_radia.eq.0) then
     3318      print *,'--------------------------------------------------'
     3319      print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
     3320      print *,'>>>>           heat et cool mis a zero '
     3321      print *,'--------------------------------------------------'
     3322      heat=0.
     3323      cool=0.
     3324      endif
     3325
     3326
    30183327c
    30193328c Ajouter la tendance des rayonnements (tous les pas)
    … …  
    34203729      ENDDO
    34213730c
     3731!==========================================================================
     3732! Sorties des tendances pour un point particulier
     3733! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
     3734! pour le debug
     3735! La valeur de igout est attribuee plus haut dans le programme
     3736!==========================================================================
     3737
     3738      if (klon.eq.1) then
     3739      write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
     3740      write(lunout,*)
     3741     s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
     3742      write(lunout,*)
     3743     s  nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
     3744      write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
     3745      do k=1,nlev
     3746         write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
     3747     s   d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
     3748     s   d_t_eva(igout,k)
     3749      enddo
     3750      write(lunout,*) 'cool,heat'
     3751      do k=1,nlev
     3752         write(lunout,*) cool(igout,k),heat(igout,k)
     3753      enddo
     3754
     3755      write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
     3756      do k=1,nlev
     3757         write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
     3758     s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
     3759      enddo
     3760
     3761      write(lunout,*) 'd_ps ',d_ps(igout)
     3762      write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
     3763      do k=1,nlev
     3764         write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
     3765     s  d_qx(igout,k,1),d_qx(igout,k,2)
     3766      enddo
     3767      endif
     3768
     3769!==========================================================================
     3770
     3771c============================================================
     3772c   Calcul de la temperature potentielle
     3773c============================================================
     3774      DO k = 1, klev
     3775      DO i = 1, klon
     3776        theta(i,k)=t(i,k)*(100000./pplay(i,k))**(RD/RCPD)
     3777      ENDDO
     3778      ENDDO
     3779c
     3780
    34223781c 22.03.04 BEG
    34233782c=============================================================

  • LMDZ4/trunk/libf/phylmd/printflag.F

    r782 r879  
    3737           PRINT *,' *****           Shema  convection  Tiedtke 
    3838     ,          ******'
    39        ELSE IF ( iflag_con.EQ. 3 )   THEN
    40            PRINT *,' *****           Shema  convection    CCM     
     39       ELSE IF ( iflag_con.GE. 3 )   THEN
     40           PRINT *,' *****           Shema  convection    Emanuel     
    4141     ,          ******'
    4242       ENDIF

  • LMDZ4/trunk/libf/phylmd/thermcell.h

    r878 r879  
    22      real r_aspect_thermals,l_mix_thermals,tho_thermals
    33      integer w2di_thermals,isplit
     4      integer iflag_coupl,iflag_clos,iflag_wake
    45
    56      common/ctherm1/iflag_thermals,nsplit_thermals
    67      common/ctherm2/r_aspect_thermals,l_mix_thermals,tho_thermals
    78      common/ctherm3/w2di_thermals
    8 
     9      common/ctherm4/iflag_coupl,iflag_clos,iflag_wake

  • LMDZ4/trunk/libf/phylmd/thermcell_main.F90

    r878 r879  
    88     &                  ,fm0,entr0,zqla,lmax  &
    99     &                  ,ratqscth,ratqsdiff,zqsatth  &
    10      &                  ,r_aspect,l_mix,w2di,tho)
     10     &                  ,r_aspect,l_mix,w2di,tho &
     11     &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th)
    1112
    1213      IMPLICIT NONE
    … …  
    5859!   ------
    5960
    60 !      integer,save :: igout=4259
    61       integer,save :: igout=2928
     61      integer,save :: igout=1521
    6262      integer,save :: lunout=6
    63       integer,save :: lev_out=10
     63      integer,save :: lev_out=1
    6464
    6565      INTEGER ig,k,l,ll
    … …  
    117117      real ztla(klon,klev),zqla(klon,klev),zqta(klon,klev)
    118118!niveau de condensation
    119       real nivcon(klon)
     119      integer nivcon(klon)
    120120      real zcon(klon)
    121121      REAL CHI
    … …  
    136136
    137137!
     138      !nouvelles variables pour la convection
     139      real Ale_bl(klon)
     140      real Alp_bl(klon)
     141      real alp_int(klon)
     142      real ale_int(klon)
     143      integer n_int(klon)
     144      real fm_tot(klon)
     145      real wght_th(klon,klev)
     146      integer lalim_conv(klon)
     147      logical therm
     148      save therm
    138149
    139150      character*2 str2
    … …  
    166177            f0(ig)=1.e-5
    167178            zmax0(ig)=40.
     179            therm=.false.
    168180      endif
    169181      enddo
    … …  
    311323      endif
    312324
     325      do ig=1,klon
     326      if (alim_star(ig,1).gt.1.e-10) then
     327         therm=.true.
     328      endif
     329      enddo
    313330!-----------------------------------------------------------------------------
    314331!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
    … …  
    460477! initialisation
    461478      do ig=1,ngrid
    462          nivcon(ig)=0.
     479         nivcon(ig)=0
    463480         zcon(ig)=0.
    464481      enddo
    … …  
    512529         enddo
    513530      enddo
     531!calcul de ale_bl et alp_bl
     532!pour le calcul d'une valeur intégrée entre la surface et lmax
     533      do ig=1,ngrid
     534      alp_int(ig)=0.
     535      ale_int(ig)=0.
     536      n_int(ig)=0
     537      enddo
     538      do ig=1,ngrid
     539!      do l=nivcon(ig),lmax(ig)
     540      do l=1,lmax(ig)
     541      alp_int(ig)=alp_int(ig)+0.5*rhobarz(ig,l)*wth3(ig,l)
     542      ale_int(ig)=ale_int(ig)+0.5*zw2(ig,l)**2
     543      n_int(ig)=n_int(ig)+1
     544      enddo
     545      enddo
     546!      print*,'avant calcul ale et alp'
     547!calcul de ALE et ALP pour la convection
     548      do ig=1,ngrid
     549!      Alp_bl(ig)=0.5*rhobarz(ig,lmix_bis(ig))*wth3(ig,lmix(ig))
     550!          Alp_bl(ig)=0.5*rhobarz(ig,nivcon(ig))*wth3(ig,nivcon(ig))
     551!      Alp_bl(ig)=0.5*rhobarz(ig,lmix(ig))*wth3(ig,lmix(ig))
     552!     &           *0.1
     553!valeur integree de alp_bl * 0.5:
     554       if (n_int(ig).gt.0) then
     555       Alp_bl(ig)=0.5*alp_int(ig)/n_int(ig)
     556!       if (Alp_bl(ig).lt.0.) then
     557!       Alp_bl(ig)=0.
     558       endif
     559!       endif
     560!         write(18,*),'rhobarz,wth3,Alp',rhobarz(ig,nivcon(ig)),
     561!     s               wth3(ig,nivcon(ig)),Alp_bl(ig)
     562!            write(18,*),'ALP_BL',Alp_bl(ig),lmix(ig)
     563!      Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
     564!      if (nivcon(ig).eq.1) then
     565!       Ale_bl(ig)=0.
     566!       else
     567!valeur max de ale_bl:
     568       Ale_bl(ig)=0.5*zw2(ig,lmix(ig))**2
     569!     & /2.
     570!     & *0.1
     571!        Ale_bl(ig)=0.5*zw2(ig,lmix_bis(ig))**2
     572!       if (n_int(ig).gt.0) then
     573!       Ale_bl(ig)=ale_int(ig)/n_int(ig)
     574!        Ale_bl(ig)=4.
     575!       endif
     576!       endif
     577!            Ale_bl(ig)=0.5*wth2(ig,lmix_bis(ig))
     578!          Ale_bl(ig)=wth2(ig,nivcon(ig))
     579!          write(19,*),'wth2,ALE_BL',wth2(ig,nivcon(ig)),Ale_bl(ig)
     580      enddo
     581!test:calcul de la ponderation des couches pour KE
     582!initialisations
     583!      print*,'ponderation'
     584      do ig=1,ngrid
     585           fm_tot(ig)=0.
     586      enddo
     587       do ig=1,ngrid
     588        do k=1,klev
     589           wght_th(ig,k)=1.
     590        enddo
     591       enddo
     592       do ig=1,ngrid
     593!         lalim_conv(ig)=lmix_bis(ig)
     594!la hauteur de la couche alim_conv = hauteur couche alim_therm
     595         lalim_conv(ig)=lalim(ig)
     596!         zentr(ig)=zlev(ig,lalim(ig))
     597      enddo
     598      do ig=1,ngrid
     599        do k=1,lalim_conv(ig)
     600           fm_tot(ig)=fm_tot(ig)+fm(ig,k)
     601        enddo
     602      enddo
     603      do ig=1,ngrid
     604        do k=1,lalim_conv(ig)
     605           if (fm_tot(ig).gt.1.e-10) then
     606!           wght_th(ig,k)=fm(ig,k)/fm_tot(ig)
     607           endif
     608!on pondere chaque couche par a*
     609             if (alim_star(ig,k).gt.1.e-10) then
     610             wght_th(ig,k)=alim_star(ig,k)
     611             else
     612             wght_th(ig,k)=1.
     613             endif
     614        enddo
     615      enddo
     616!      print*,'apres wght_th'
     617!test pour prolonger la convection
     618      do ig=1,ngrid
     619      if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
     620      lalim_conv(ig)=1
     621      wght_th(ig,1)=1.
     622!      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
     623      endif
     624      enddo
     625
    514626!calcul du ratqscdiff
    515627      print*,'14e OK convect8'
    … …  
    583695      character*21 comment
    584696
    585       print*,'TEST ',comment
     697      print*,'WARNING !!! TEST ',comment
     698      return
     699
    586700!  test sur la hauteur des thermiques ...
    587701         do i=1,klon

  • LMDZ4/trunk/libf/phylmd/write_histday.h

    r878 r879  
    1111         itau_w = itau_phy + itap
    1212
     13c  temperature tendency due to moist convective processes
     14       DO l=1, klev
     15       DO i=1, klon
     16        zx_tmp_fi3d(i,l)=d_t_con(i,l)/pdtphys
     17       ENDDO !i
     18       ENDDO !l
    1319c
    1420      IF(lev_histday.GE.1) THEN
    … …  
    3440      CALL histwrite_phy(nid_day,"weakinv",itau_w,weak_inversion)
    3541
     42cIM: 101003 : K/30min ==> K/s
    3643
    3744cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, zxtsol,zx_tmp_2d)
    … …  
    320327      CALL histwrite_phy(nid_day,"pres",itau_w,pplay)
    321328c
     329      CALL histwrite_phy(nid_day,"rhum",itau_w,zx_rh)
     330
     331      CALL histwrite_phy(nid_day,"wh",itau_w,wake_h)
     332c
     333      CALL histwrite_phy(nid_day,"ws",itau_w,wake_s)
     334c
     335      CALL histwrite_phy(nid_day,"cin",itau_w,cin)
     336
     337      CALL histwrite_phy(nid_day,"wale",itau_w,ale_wake)
     338c
     339      CALL histwrite_phy(nid_day,"walp",itau_w,alp_wake)
     340c
     341      CALL histwrite_phy(nid_day,"blale",itau_w,Ale_bl)
     342c
     343      CALL histwrite_phy(nid_day,"blalp",itau_w,Alp_bl)
     344
     345cFin Wake parameters 2D     
     346c Wake parameters 3D
     347c
     348          print*,'SSS test2'
     349c  if (1.eq.0) then
     350      CALL histwrite_phy(nid_day,"wdt1",itau_w,wake_deltat(:,1))
     351c
     352      CALL histwrite_phy(nid_day,"wdq1",itau_w,wake_deltaq(:,1))
     353c
     354
     355c
    322356      ENDIF !lev_histday.GE.3
    323357c=================================================================
    … …  
    554588c
    555589cIM: 101003 : K/30min ==> K/s
     590      zx_tmp_fi3d(1:klon,1:klev)=
     591     s   (d_q_vdf(1:klon,1:klev)+d_q_ajs(1:klon,1:klev))/pdtphys
     592      CALL histwrite_phy(nid_day,"dqpbl",itau_w,zx_tmp_fi3d)
     593
     594      CALL histwrite_phy(nid_day,"dtwak",itau_w,dt_wake)
     595
     596          CALL histwrite_phy(nid_day,"cape",itau_w,cape)
     597c
    556598      zx_tmp_fi3d(1:klon,1:klev)=d_t_con(1:klon,1:klev)/pdtphys
    557599cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)

  • LMDZ4/trunk/libf/phylmd/write_histmth.h

    r878 r879  
    99         itau_w = itau_phy + itap
    1010
     11c  temperature tendency due to moist convective processes
     12       DO l=1, klev
     13       DO i=1, klon
     14        zx_tmp_fi3d(i,l)=d_t_con(i,l)/pdtphys
     15       ENDDO !i
     16       ENDDO !l
     17c
     18
     19cIM: 101003 : K/30min ==> K/s
     20      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
     21cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     22      CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d)
     23c
     24cIM: 101003 : K/30min ==> K/s
     25      zx_tmp_fi3d(1:klon,1:klev)=d_t_con(1:klon,1:klev)/pdtphys
     26cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     27      CALL histwrite_phy(nid_mth,"dtcon",itau_w,zx_tmp_fi3d)
     28c
    1129c
    1230      IF(type_run.EQ."CLIM".OR.type_run.EQ."ENSP") THEN
    … …  
    4462c ENSEMBLES BEG
    4563cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zt2m,zx_tmp_2d)
    46       CALL histwrite_phy(nid_mth,"t2m_min",itau_w,zt2m)
     64c     CALL histwrite_phy(nid_mth,"t2m_min",itau_w,zt2m)
    4765c
    4866cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zt2m,zx_tmp_2d)
    49       CALL histwrite_phy(nid_mth,"t2m_max",itau_w,zt2m)
     67c     CALL histwrite_phy(nid_mth,"t2m_max",itau_w,zt2m)
    5068c
    5169c     CALL gr_fi_ecrit(1,klon,iim,jjmp1,ftsoil(:,1,is_ter),zx_tmp_2d)
    … …  
    103121c
    104122cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, nday_rain,zx_tmp_2d)
    105       CALL histwrite_phy(nid_mth,"ndayrain",itau_w,nday_rain)
     123c     CALL histwrite_phy(nid_mth,"ndayrain",itau_w,nday_rain)
    106124c
    107125      DO i = 1, klon
    … …  
    968986      CALL histwrite_phy(nid_mth,"dtlw",itau_w,zx_tmp_fi3d)
    969987c
    970 c  temperature tendency due to moist convective processes
    971        DO l=1, klev
    972        DO i=1, klon
    973         zx_tmp_fi3d(i,l)=d_t_con(i,l)/pdtphys
    974        ENDDO !i
    975        ENDDO !l
    976 c
    977 cym      CALL gr_fi_ecrit(klev, klon,iim,jjmp1, zx_tmp_fi3d,zx_tmp_3d)
    978       CALL histwrite_phy(nid_mth,"dtcon",itau_w,zx_tmp_fi3d)
     988
     989cIM: 101003 : K/30min ==> K/s
     990      zx_tmp_fi3d(1:klon,1:klev)=d_t_ajsb(1:klon,1:klev)/pdtphys
     991cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     992      CALL histwrite_phy(nid_mth,"dtajs",itau_w,zx_tmp_fi3d)
     993
     994c
     995      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
     996cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
     997      CALL histwrite_phy(nid_mth,"dqthe",itau_w,zx_tmp_fi3d)
    979998c
    980999
    … …  
    10781097c ENSEMBLES BEG
    10791098cym      CALL gr_fi_ecrit(1,klon,iim,jjmp1,zt2m,zx_tmp_2d)
    1080       CALL histwrite_phy(nid_mth,"t2m_min",itau_w,zt2m)
     1099c     CALL histwrite_phy(nid_mth,"t2m_min",itau_w,zt2m)
    10811100c
    10821101cym     CALL gr_fi_ecrit(1,klon,iim,jjmp1,zt2m,zx_tmp_2d)
    1083       CALL histwrite_phy(nid_mth,"t2m_max",itau_w,zt2m)
     1102c     CALL histwrite_phy(nid_mth,"t2m_max",itau_w,zt2m)
    10841103c
    10851104c     CALL gr_fi_ecrit(1,klon,iim,jjmp1,ftsoil(:,1,is_ter),zx_tmp_2d)
    … …  
    11371156c
    11381157cym      CALL gr_fi_ecrit(1, klon,iim,jjmp1, nday_rain,zx_tmp_2d)
    1139       CALL histwrite_phy(nid_mth,"ndayrain",itau_w,nday_rain)
     1158c     CALL histwrite_phy(nid_mth,"ndayrain",itau_w,nday_rain)
    11401159c
    11411160      DO i = 1, klon
    … …  
    18061825c
    18071826cIM: 101003 : K/30min ==> K/s
    1808       zx_tmp_fi3d(1:klon,1:klev)=d_t_ajs(1:klon,1:klev)/pdtphys
    1809 cym      CALL gr_fi_ecrit(klev,klon,iim,jjmp1,zx_tmp_fi3d,zx_tmp_3d)
    1810       CALL histwrite_phy(nid_mth,"dtthe",itau_w,zx_tmp_fi3d)
    18111827c
    18121828      zx_tmp_fi3d(1:klon,1:klev)=d_q_ajs(1:klon,1:klev)/pdtphys
Note: See TracChangeset for help on using the changeset viewer.