Changeset 1943


Ignore:
Timestamp:
Jan 22, 2014, 10:51:36 AM (10 years ago)
Author:
idelkadi
Message:

Changes concerinng the Hourdin et al. 2002 version of the thermals and
Ayotte cases

calltherm.F90

call to thermcell_2002 modified :
1) iflag_thermals changed from 1 to >= 1000 ; iflag_thermals-1000
controls sub options.
2) thermals w and fraction in output files.

hbtm.F

Singularity in the 1/heat dependency of the Monin Obukov length L when
heat=0. Since 1/L is used rather than L, it is preferable to compute
directly L. There is a dependency in 1/u* then which is treated with a
threshold value.
(+ some cleaning in the syntax).

lmdz1d.F

If nday<0, -nday is used as the total number of time steps of the
simulation.
The option with imposed wtsurf and wqsurf read in lmdz1d.def was not
active anymore.
< IF(.NOT.ok_flux_surf) THEN
changed to

IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN

before

fsens=-wtsurf*rcpd*rho(1)
flat=-wqsurf*rlvtt*rho(1)

phys_output_write.F90 et phys_local_var_mod.F90

Removing the d_u_ajsb contribution to duthe (same for dv).
Those tendencies are not computed by the model ...
< zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
< d_u_ajsb(1:klon,1:klev)/pdtphys
---

zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys

! d_u_ajsb(1:klon,1:klev)/pdtphys

thermcell_dq.F90, thermcell_main.F90

Some cleaning

thermcell_old.F

Control by iflag_thermals >= 1000
wa_moy and fraca in outputs
+ cleaning

Location:
LMDZ5/trunk/libf
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phy1d/lmdz1d.F

    r1914 r1943  
    357357c      Le numero du jour est dans "day". L heure est traitee separement.
    358358c      La date complete est dans "daytime" (l'unite est le jour).
    359       fnday=nday
     359      if (nday>0) then
     360         fnday=nday
     361      else
     362         fnday=-nday/float(day_step)
     363      endif
     364
    360365c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
    361366      IF(forcing_type .EQ. 61) fnday=53100./86400.
    … …  
    370375      call ymds2ju(annee_ref,mois,day_ref,heure,day)
    371376      day_ini = day
    372       day_end = day_ini + nday
     377      day_end = day_ini + fnday
    373378
    374379      IF (forcing_type .eq.2) THEN
    … …  
    463468!! mpl et jyg le 22/08/2012 :
    464469!!  pour que les cas a flux de surface imposes marchent
    465       IF(.NOT.ok_flux_surf) THEN
     470      IF(.NOT.ok_flux_surf.or.max(abs(wtsurf),abs(wqsurf))>0.) THEN
    466471       fsens=-wtsurf*rcpd*rho(1)
    467472       flat=-wqsurf*rlvtt*rho(1)
    468473       print *,'Flux: ok_flux wtsurf wqsurf',ok_flux_surf,wtsurf,wqsurf
    469474      ENDIF
     475      print*,'Flux sol ',fsens,flat
    470476!!      ok_flux_surf=.false.
    471477!!      fsens=-wtsurf*rcpd*rho(1)

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

    r1907 r1943  
    174174         do isplit=1,nsplit_thermals
    175175
    176           if (iflag_thermals.eq.1) then
    177             CALL thermcell_2002(klon,klev,zdt   &
     176          if (iflag_thermals>=1000) then
     177            CALL thermcell_2002(klon,klev,zdt,iflag_thermals   &
    178178     &      ,pplay,paprs,pphi  &
    179179     &      ,u_seri,v_seri,t_seri,q_seri  &
    180180     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
    181      &      ,zfm_therm,zentr_therm  &
     181     &      ,zfm_therm,zentr_therm,fraca,zw2  &
    182182     &      ,r_aspect_thermals,30.,w2di_thermals  &
    183183     &      ,tau_thermals)
    … …  
    273273       flag_bidouille_stratocu=iflag_thermals<=12.or.iflag_thermals==14.or.iflag_thermals==16.or.iflag_thermals==18
    274274
    275       if (iflag_thermals<=12) then
    276          lmax=1
     275! Calcul a posteriori du niveau max des thermiques pour les schémas qui
     276! ne la sortent pas.
     277      if (iflag_thermals<=12.or.iflag_thermals>=1000) then
     278         lmax(:)=1
    277279         do k=1,klev-1
    278280            zdetr_therm(:,k)=zentr_therm(:,k)+zfm_therm(:,k)-zfm_therm(:,k+1)
     281         enddo
     282         do k=1,klev-1
     283            do i=1,klon
     284               if (zfm_therm(i,k+1)>0.) lmax(i)=k
     285            enddo
    279286         enddo
    280287      endif

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

    r1907 r1943  
    6767      INTEGER isommet
    6868cum      PARAMETER (isommet=klev) ! limite max sommet pbl
    69       REAL vk
    70       PARAMETER (vk=0.35)     ! Von Karman => passer a .41 ! cf U.Olgstrom
    71       REAL ricr
    72       PARAMETER (ricr=0.4)
    73       REAL fak
    74       PARAMETER (fak=8.5)     ! b calcul du Prandtl et de dTetas
    75       REAL fakn
    76       PARAMETER (fakn=7.2)    ! a
    77       REAL onet
    78       PARAMETER (onet=1.0/3.0)
    79       REAL t_coup
    80       PARAMETER(t_coup=273.15)
    81       REAL zkmin
    82       PARAMETER (zkmin=0.01)
    83       REAL betam
    84       PARAMETER (betam=15.0)  ! pour Phim / h dans la S.L stable
    85       REAL betah
    86       PARAMETER (betah=15.0)
    87       REAL betas
    88       PARAMETER (betas=5.0)   ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
    89       REAL sffrac
    90       PARAMETER (sffrac=0.1)  ! S.L. = z/h < .1
    91       REAL binm
    92       PARAMETER (binm=betam*sffrac)
    93       REAL binh
    94       PARAMETER (binh=betah*sffrac)
    95       REAL ccon
    96       PARAMETER (ccon=fak*sffrac*vk)
    97 c
     69      REAL, PARAMETER :: vk=0.35  ! Von Karman => passer a .41 ! cf U.Olgstrom
     70      REAL, PARAMETER :: ricr=0.4
     71      REAL, PARAMETER :: fak=8.5     ! b calcul du Prandtl et de dTetas
     72      REAL, PARAMETER :: fakn=7.2    ! a
     73      REAL, PARAMETER :: onet=1.0/3.0
     74      REAL, PARAMETER:: t_coup=273.15
     75      REAL, PARAMETER :: zkmin=0.01
     76      REAL, PARAMETER :: betam=15.0  ! pour Phim / h dans la S.L stable
     77      REAL, PARAMETER :: betah=15.0
     78      REAL, PARAMETER :: betas=5.0 ! Phit dans la S.L. stable (mais 2 formes / z/OBL<>1
     79      REAL, PARAMETER :: sffrac=0.1  ! S.L. = z/h < .1
     80      REAL, PARAMETER :: usmin=1.e-12
     81      REAL, PARAMETER :: binm=betam*sffrac
     82      REAL, PARAMETER :: binh=betah*sffrac
     83      REAL, PARAMETER :: ccon=fak*sffrac*vk
     84      REAL, PARAMETER :: b1=70.,b2=20.
     85      REAL, PARAMETER :: zref=2.  ! Niveau de ref a 2m peut eventuellement
     86c                              etre choisi a 10m
    9887      REAL q_star,t_star
    99       REAL b1,b2,b212,b2sr     ! Lambert correlations T' q' avec T* q*
    100       PARAMETER (b1=70.,b2=20.)
     88      REAL b212,b2sr     ! Lambert correlations T' q' avec T* q*
    10189c
    10290      REAL z(klon,klev)
    10391cAM      REAL pcfm(klon,klev), pcfh(klon,klev)
    104 cAM
    105       REAL zref
    106       PARAMETER (zref=2.)    ! Niveau de ref a 2m peut eventuellement
    107 c                              etre choisi a 10m
    108 cMA
    109 c
    11092      INTEGER i, k, j
    11193      REAL zxt
    … …  
    129111cAM      REAL cgq(klon,2:klev) ! counter-gradient term for constituents
    130112cAM      REAL cgs(klon,2:klev) ! counter-gradient star (cg/flux)
    131       REAL obklen(klon)     ! Monin-Obukhov lengh
     113      REAL unsobklen(klon)     ! Monin-Obukhov lengh
    132114cAM      REAL ztvd, ztvu,
    133115      REAL zdu2
    … …  
    345327         plcl(i) = 6000.
    346328c Lambda = -u*^3 / (alpha.g.kvon.<w'Theta'v>
    347          obklen(i) = -t(i,1)*ustar(i)**3/(RG*vk*heatv(i))
     329         unsobklen(i) = -RG*vk*heatv(i)/(t(i,1)*max(ustar(i),usmin)**3)
    348330         trmb1(i)   = 0.
    349331         trmb2(i)   = 0.
    … …  
    420402      DO i = 1, knon
    421403      IF (check(i)) THEN
    422         phiminv(i) = (1.-binm*pblh(i)/obklen(i))**onet
     404        phiminv(i) = (1.-binm*pblh(i)*unsobklen(i))**onet
    423405c ***************************************************
    424406c Wm ? et W* ? c'est la formule pour z/h < .1
    … …  
    621603          zxt=(Th_th(i)-zref*0.5*RG/RCPD/(1.+RVTMP2*qT_th(i)))
    622604     .         *(1.+RETV*qT_th(i))
    623           phiminv(i) = (1. - binm*pblh(i)/obklen(i))**onet
    624           phihinv(i) = sqrt(1. - binh*pblh(i)/obklen(i))
     605          phiminv(i) = (1. - binm*pblh(i)*unsobklen(i))**onet
     606          phihinv(i) = sqrt(1. - binh*pblh(i)*unsobklen(i))
    625607          wm(i)      = ustar(i)*phiminv(i)
    626608          fak2(i)    = wm(i)*pblh(i)*vk
    … …  
    657639
    658640            zh(i) = zmzp/pblh(i)
    659             zl(i) = zmzp/obklen(i)
     641            zl(i) = zmzp*unsobklen(i)
    660642            zzh(i) = 0.
    661643            IF (zh(i).LE.1.0) zzh(i) = (1. - zh(i))**2

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

    r1924 r1943  
    4141      REAL, SAVE, ALLOCATABLE :: d_u_ajs(:,:), d_v_ajs(:,:)
    4242      !$OMP THREADPRIVATE(d_u_ajs, d_v_ajs)
    43       REAL, SAVE, ALLOCATABLE :: d_u_ajsb(:,:), d_v_ajsb(:,:)
    44       !$OMP THREADPRIVATE(d_u_ajsb, d_v_ajsb)
    4543      REAL, SAVE, ALLOCATABLE :: d_t_eva(:,:),d_q_eva(:,:)
    4644      !$OMP THREADPRIVATE(d_t_eva,d_q_eva)
    … …  
    303301      allocate(d_t_ajs(klon,klev),d_q_ajs(klon,klev))
    304302      allocate(d_u_ajs(klon,klev),d_v_ajs(klon,klev))
    305       allocate(d_u_ajsb(klon,klev),d_v_ajsb(klon,klev))
    306303      allocate(d_t_eva(klon,klev),d_q_eva(klon,klev))
    307304      allocate(d_t_lscst(klon,klev),d_q_lscst(klon,klev))
    … …  
    456453      deallocate(d_t_ajs,d_q_ajs)
    457454      deallocate(d_u_ajs,d_v_ajs)
    458       deallocate(d_u_ajsb,d_v_ajsb)
    459455      deallocate(d_t_eva,d_q_eva)
    460456      deallocate(d_t_lscst,d_q_lscst)

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

    r1938 r1943  
    206206         rneb, rnebjn, zx_rh, frugs, agesno, d_t_dyn, d_q_dyn, &
    207207         d_u_dyn, d_v_dyn, d_t_con, d_t_ajsb, d_t_ajs, &
    208          d_u_ajsb, d_u_ajs, d_v_ajsb, d_v_ajs, &
     208         d_u_ajs, d_v_ajs, &
    209209         d_u_con, d_v_con, d_q_con, d_q_ajs, d_t_lsc, &
    210210         d_t_eva, d_q_lsc, beta_prec, d_t_lscth, &
    … …  
    10461046       CALL histwrite_phy(o_dtthe, zx_tmp_fi3d)
    10471047       IF (vars_defined) THEN
    1048           zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys - &
    1049                d_u_ajsb(1:klon,1:klev)/pdtphys
     1048          zx_tmp_fi3d(1:klon,1:klev)=d_u_ajs(1:klon,1:klev)/pdtphys
    10501049       ENDIF
    10511050       CALL histwrite_phy(o_duthe, zx_tmp_fi3d)
    10521051       IF (vars_defined) THEN
    1053           zx_tmp_fi3d(1:klon,1:klev)=d_v_ajs(1:klon,1:klev)/pdtphys - &
    1054                d_v_ajsb(1:klon,1:klev)/pdtphys
     1052          zx_tmp_fi3d(1:klon,1:klev)=d_v_ajs(1:klon,1:klev)/pdtphys
    10551053       ENDIF
    10561054       CALL histwrite_phy(o_dvthe, zx_tmp_fi3d)

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

    r1907 r1943  
    5353            cfl=max(cfl,fm(ig,k)/zzm)
    5454            if (entr(ig,k).gt.zzm) then
    55                print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
    56                abort_message = ''
     55               print*,'entr*dt>m,1',k,entr(ig,k)*ptimestep,masse(ig,k)
     56               abort_message = 'entr dt > m, 1st'
    5757               CALL abort_gcm (modname,abort_message,1)
    5858            endif
    … …  
    188188            cfl=max(cfl,fm(ig,k)/zzm)
    189189            if (entr(ig,k).gt.zzm) then
    190                print*,'entr dt > m ',entr(ig,k)*ptimestep,masse(ig,k)
    191                abort_message = ''
     190               print*,'entr*dt>m,2',k,entr(ig,k)*ptimestep,masse(ig,k)
     191               abort_message = 'entr dt > m, 2nd'
    192192               CALL abort_gcm (modname,abort_message,1)
    193193            endif

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

    r1907 r1943  
    8383      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    8484      real pphi(ngrid,nlay)
     85      LOGICAL debut
    8586
    8687!   local:
    … …  
    178179!FH/IM     save f0
    179180      real zlevinter(klon)
    180       logical debut
    181181       real seuil
    182182      real csc(klon,klev)
    … …  
    256256!
    257257
    258       seuil=0.25
    259 
    260       if (debut)  then
    261 !        call getin('dvdq',dvdq)
    262 !        call getin('dqimpl',dqimpl)
    263 
    264          if (iflag_thermals==15.or.iflag_thermals==16) then
    265             dvdq=0
    266             dqimpl=-1
    267          else
    268             dvdq=1
    269             dqimpl=1
    270          endif
    271 
    272          fm0=0.
    273          entr0=0.
    274          detr0=0.
    275 
    276 
    277 #undef wrgrads_thermcell
    278 #ifdef wrgrads_thermcell
    279 ! Initialisation des sorties grads pour les thermiques.
    280 ! Pour l'instant en 1D sur le point igout.
    281 ! Utilise par thermcell_out3d.h
    282          str10='therm'
    283          call inigrads(1,1,rlond(igout),1.,-180.,180.,jjm, &
    284      &   rlatd(igout),-90.,90.,1.,llm,pplay(igout,:),1.,   &
    285      &   ptimestep,str10,'therm ')
    286 #endif
    287 
    288 
    289 
    290       endif
    291 
    292       fm=0. ; entr=0. ; detr=0.
    293 
    294 
    295       icount=icount+1
     258   seuil=0.25
     259
     260   if (debut) then
     261      if (iflag_thermals==15.or.iflag_thermals==16) then
     262         dvdq=0
     263         dqimpl=-1
     264      else
     265         dvdq=1
     266         dqimpl=1
     267      endif
     268
     269      fm0=0.
     270      entr0=0.
     271      detr0=0.
     272   endif
     273   fm=0. ; entr=0. ; detr=0.
     274   icount=icount+1
    296275
    297276!IM 090508 beg

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

    r1907 r1943  
    1       SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep
     1      SUBROUTINE thermcell_2002(ngrid,nlay,ptimestep,iflag_thermals
    22     s                  ,pplay,pplev,pphi
    33     s                  ,pu,pv,pt,po
    44     s                  ,pduadj,pdvadj,pdtadj,pdoadj
    5      s                  ,fm0,entr0
    6 c    s                  ,pu_therm,pv_therm
     5     s                  ,fm0,entr0,fraca,wa_moy
    76     s                  ,r_aspect,l_mix,w2di,tho)
    87
    98      USE dimphy
     9      USE write_field_phy
    1010      IMPLICIT NONE
    1111
    … …  
    3535
    3636#include "dimensions.h"
    37 cccc#include "dimphy.h"
    3837#include "YOMCST.h"
    3938
    … …  
    4140c   ----------
    4241
    43       INTEGER ngrid,nlay,w2di
     42      INTEGER ngrid,nlay,w2di,iflag_thermals
    4443      REAL tho
    4544      real ptimestep,l_mix,r_aspect
    … …  
    5049      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    5150      real pphi(ngrid,nlay)
    52 
    53       integer idetr
    54       save idetr
    55       data idetr/3/
    56 c$OMP THREADPRIVATE(idetr)
     51      real fraca(ngrid,nlay+1),zw2(ngrid,nlay+1)
     52
     53      integer,save :: idetr=3,lev_out=1
     54c$OMP THREADPRIVATE(idetr,lev_out)
    5755
    5856c   local:
    5957c   ------
    6058
     59      INTEGER, SAVE :: dvdq=0,flagdq=0,dqimpl=1
     60      LOGICAL, SAVE :: debut=.true.
     61!$OMP THREADPRIVATE(dvdq,flagdq,debut,dqimpl)
     62
    6163      INTEGER ig,k,l,lmax(klon,klev+1),lmaxa(klon),lmix(klon)
    62       real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
     64      real zmax(klon),zw,zz,ztva(klon,klev),zzz
    6365
    6466      real zlev(klon,klev+1),zlay(klon,klev)
    … …  
    7981      real zha(klon,klev)
    8082      real wa_moy(klon,klev+1)
    81       real fraca(klon,klev+1)
    8283      real fracc(klon,klev+1)
    8384      real zf,zf2
    … …  
    8687
    8788      real count_time
    88       integer isplit,nsplit,ialt
    89       parameter (nsplit=10)
    90       data isplit/0/
    91       save isplit
    92 c$OMP THREADPRIVATE(isplit)
    9389
    9490      logical sorties
    … …  
    142138c   ---------------------------------------------------
    143139
    144       print*,'0 OK convect8'
     140!     print*,'0 OK convect8'
    145141
    146142      DO 1010 l=1,nlay
    … …  
    1551511010  CONTINUE
    156152
    157 c     print*,'1 OK convect8'
     153!     print*,'1 OK convect8'
    158154c                       --------------------
    159155c
    … …  
    177173c-----------------------------------------------------------------------
    178174
     175      if (debut)  then
     176        flagdq=(iflag_thermals-1000)/100
     177        dvdq=(iflag_thermals-(1000+flagdq*100))/10
     178        if (flagdq==2) dqimpl=-1
     179        if (flagdq==3) dqimpl=1
     180        debut=.false.
     181      endif
     182      print*,'TH flag th ',iflag_thermals,flagdq,dvdq,dqimpl
     183
    179184      do l=2,nlay
    180185         do ig=1,ngrid
    … …  
    192197      enddo
    193198
    194 c     print*,'2 OK convect8'
     199!     print*,'2 OK convect8'
    195200c-----------------------------------------------------------------------
    196201c   Calcul des densites
    … …  
    217222      enddo
    218223
    219 c     print*,'3 OK convect8'
     224!     print*,'3 OK convect8'
    220225c------------------------------------------------------------------
    221226c   Calcul de w2, quarre de w a partir de la cape
    … …  
    274279      enddo
    275280
    276 c     print*,'4 OK convect8'
     281!     print*,'4 OK convect8'
    277282c Calcul de la couche correspondant a la hauteur du thermique
    278283      do k=1,nlay-1
    … …  
    287292      enddo
    288293
    289 c     print*,'5 OK convect8'
     294!     print*,'5 OK convect8'
    290295c   Calcule du w max du thermique
    291296      do k=1,nlay
    … …  
    315320      enddo
    316321
    317 c     print*,'6 OK convect8'
     322!     print*,'6 OK convect8'
    318323c   Longueur caracteristique correspondant a la hauteur des thermiques.
    319324      do  ig=1,ngrid
    … …  
    330335c      call dump2d(iim,jjm-1,zmax(2:ngrid-1),'ZMAX      ')
    331336
     337!     print*,'OKl336'
    332338c   Calcul de l'entrainement.
    333339c   Le rapport d'aspect relie la largeur de l'ascendance a l'epaisseur
    … …  
    348354      enddo
    349355
    350 c     print*,'7 OK convect8'
     356
     357!     print*,'7 OK convect8'
    351358      do k=1,klev+1
    352359         do ig=1,ngrid
    … …  
    359366      enddo
    360367
    361 c     print*,'8 OK convect8'
     368!     print*,'8 OK convect8'
    362369      do ig=1,ngrid
    363370         lmaxa(ig)=1
    … …  
    367374
    368375
     376!     print*,'OKl372'
    369377      do l=1,nlay-2
    370378         do ig=1,ngrid
    … …  
    422430      enddo
    423431
    424 c     print*,'9 OK convect8'
     432!     print*,'9 OK convect8'
    425433c     print*,'WA1 ',wa_moy
    426434
    … …  
    433441c   de la vitesse d'entrainement horizontal dans la couche alimentante.
    434442
     443!     print*,'OKl439'
    435444      do l=2,nlay
    436445         do ig=1,ngrid
    … …  
    463472       enddo
    464473
    465 c     print*,'10 OK convect8'
     474!     print*,'10 OK convect8'
    466475c     print*,'WA2 ',wa_moy
    467476c   calcul de la fraction de la maille concernée par l'ascendance en tenant
    … …  
    500509      enddo
    501510
    502 c     print*,'11 OK convect8'
     511!     print*,'11 OK convect8'
    503512c     print*,'Ea3 ',wa_moy
    504513c------------------------------------------------------------------
    … …  
    535544      enddo
    536545
    537 c     print*,'12 OK convect8'
     546!     print*,'12 OK convect8'
    538547c     print*,'WA4 ',wa_moy
    539548cc------------------------------------------------------------------
    … …  
    589598
    5905994444   continue
     600!     print*,'OK 444 '
    591601
    592602      if (w2di.eq.1) then
    … …  
    598608      endif
    599609
    600       if (1.eq.1) then
     610      if (flagdq==0) then
    601611         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    602612     .    ,zh,zdhadj,zha)
    603613         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    604614     .    ,zo,pdoadj,zoa)
    605       else
     615         print*,'THERMALS OPT 1'
     616      else if (flagdq==1) then
    606617         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
    607618     .    ,zh,zdhadj,zha)
    608619         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
    609620     .    ,zo,pdoadj,zoa)
     621         print*,'THERMALS OPT 2'
     622      else
     623         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,
     624     .                    zh,zdhadj,zha,lev_out)
     625         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,
     626     .                   zo,pdoadj,zoa,lev_out)
     627         print*,'THERMALS OPT 3',dqimpl
    610628      endif
    611629
    612       if (1.eq.0) then
    613          call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
    614      .    ,fraca,zmax
    615      .    ,zu,zv,pduadj,pdvadj,zua,zva)
    616       else
     630      print*,'TH VENT ',dvdq
     631      if (dvdq==0) then
     632!     print*,'TH VENT OK ',dvdq
    617633         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    618634     .    ,zu,pduadj,zua)
    619635         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    620636     .    ,zv,pdvadj,zva)
     637      else if (dvdq==1) then
     638         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
     639     .    ,fraca,zmax
     640     .    ,zu,zv,pduadj,pdvadj,zua,zva)
     641      else if (dvdq==2) then
     642         call thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse
     643     &    ,fraca,zmax
     644     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
     645      else if (dvdq==3) then
     646         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse
     647     &    ,zu,pduadj,zua,lev_out)
     648         call thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse
     649     &    ,zv,pdvadj,zva,lev_out)
    621650      endif
     651
     652!     CALL writefield_phy('duadj',pduadj,klev)
    622653
    623654      do l=1,nlay
    … …  
    632663
    633664
    634 c     print*,'13 OK convect8'
     665!     print*,'13 OK convect8'
    635666c     print*,'WA5 ',wa_moy
    636667      do l=1,nlay
    … …  
    654685c      enddo
    655686
    656 c     print*,'14 OK convect8'
     687!     print*,'14 OK convect8'
    657688c------------------------------------------------------------------
    658689c   Calculs pour les sorties
    … …  
    679710         enddo
    680711      enddo
    681 
    682 c     print*,'15 OK convect8'
    683 
    684       isplit=isplit+1
    685 
    686 
    687 c #define und
    688         goto 123
    689 #ifdef und
    690       CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
    691       CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
    692       CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
    693       CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
    694       CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
    695       CALL writeg1d(1,nlay,zla,'la      ','la      ')
    696       CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
    697       CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
    698       CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
    699       CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
    700       CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
    701       CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
    702       CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
    703       CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
    704       CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
    705       CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
    706       CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
    707       CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
    708       CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
    709       CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
    710       CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
    711       CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
    712       CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
    713       CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
    714 
    715       CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
    716       CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
    717       CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
    718 c   recalcul des flux en diagnostique...
    719 c     print*,'PAS DE TEMPS ',ptimestep
    720        call dt2F(pplev,pplay,pt,pdtadj,wh)
    721       CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
    722 #endif
    723 123   continue
    724 ! #define troisD
    725 #ifdef troisD
    726 c       if (sorties) then
    727       print*,'Debut des wrgradsfi'
    728 
    729 c      print*,'16 OK convect8'
    730          call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
    731          call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
    732          call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
    733          call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
    734          call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
    735          call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
    736 c      print*,'WA6 ',wa_moy
    737          call wrgradsfi(1,nlay,zla,'la        ','la        ')
    738          call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
    739          call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
    740          call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
    741          call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
    742          call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
    743          call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
    744          call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
    745          call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
    746          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    747          call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
    748          call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
    749          call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
    750          call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
    751          call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    752          call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
    753          call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
    754          call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
    755          call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
    756          call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
    757          call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
    758          call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
    759          call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
    760          call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
    761          call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
    762          call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
    763 
    764          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    765          call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
    766          call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
    767 
    768 
    769 c      print*,'17 OK convect8'
    770 
    771          do k=1,klev/10
    772             write(str2,'(i2.2)') k
    773             str10='wa'//str2
    774             do l=1,nlay
    775                do ig=1,ngrid
    776                   zsortie(ig,l)=wa(ig,k,l)
    777                enddo
    778             enddo   
    779             CALL wrgradsfi(1,nlay,zsortie,str10,str10)
    780             do l=1,nlay
    781                do ig=1,ngrid
    782                   zsortie(ig,l)=larg_part(ig,k,l)
    783                enddo
    784             enddo
    785             str10='la'//str2
    786             CALL wrgradsfi(1,nlay,zsortie,str10,str10)
    787          enddo
    788 
    789 
    790 c     print*,'18 OK convect8'
    791 c      endif
    792       print*,'Fin des wrgradsfi'
    793 #endif
    794 
    795712      endif
    796713
     714!     print*,'15 OK convect8'
     715
     716
    797717c     if(wa_moy(1,4).gt.1.e-10) stop
    798718
    799 c     print*,'19 OK convect8'
     719!     print*,'19 OK convect8'
    800720      return
    801721      end
    … …  
    925845      real ratqsdiff(klon,klev)
    926846      real count_time
    927       integer isplit,nsplit,ialt
    928       parameter (nsplit=10)
    929       data isplit/0/
    930       save isplit
    931 c$OMP THREADPRIVATE(isplit)
     847      integer ialt
    932848
    933849      logical sorties
    … …  
    21972113      enddo
    21982114
    2199       print*,'12 OK convect8'
     2115!     print*,'12 OK convect8'
    22002116c     print*,'WA4 ',wa_moy
    22012117cc------------------------------------------------------------------
    … …  
    23472263c      enddo
    23482264
    2349       print*,'14 OK convect8'
     2265!     print*,'14 OK convect8'
    23502266c------------------------------------------------------------------
    23512267c   Calculs pour les sorties
    … …  
    24332349      enddo
    24342350      enddo     
    2435 cdeja fait
    2436 c      do l=1,nlay
    2437 c         do ig=1,ngrid
    2438 c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
    2439 c            if (detr(ig,l).lt.0.) then
    2440 c                entr(ig,l)=entr(ig,l)-detr(ig,l)
    2441 c                detr(ig,l)=0.
    2442 c     print*,'WARNING !!! detrainement negatif ',ig,l
    2443 c            endif
    2444 c         enddo
    2445 c      enddo
    2446 
    2447 c     print*,'15 OK convect8'
    2448 
    2449       isplit=isplit+1       
    2450 
    2451 
    2452 c #define und
    2453         goto 123
    2454 #ifdef und
    2455       CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
    2456       CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
    2457       CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
    2458       CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
    2459       CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
    2460       CALL writeg1d(1,nlay,zla,'la      ','la      ')
    2461       CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
    2462       CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
    2463       CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
    2464       CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
    2465       CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
    2466       CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
    2467       CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
    2468       CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
    2469       CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
    2470       CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
    2471       CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
    2472       CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
    2473       CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
    2474       CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
    2475       CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
    2476       CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
    2477       CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
    2478       CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
    2479 
    2480       CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
    2481       CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
    2482       CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
    2483 
    2484 c   recalcul des flux en diagnostique...
    2485 c     print*,'PAS DE TEMPS ',ptimestep
    2486        call dt2F(pplev,pplay,pt,pdtadj,wh)
    2487       CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
    2488 #endif
    2489 123   continue
    2490 ! #define troisD
    2491 #ifdef troisD
    2492 c       if (sorties) then
    2493       print*,'Debut des wrgradsfi'
    2494 
    2495 c      print*,'16 OK convect8'
    2496 c         call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
    2497 c         call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
    2498          call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
    2499          call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
    2500 c         call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
    2501 c         call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
    2502 c      print*,'WA6 ',wa_moy
    2503 c         call wrgradsfi(1,nlay,zla,'la        ','la        ')
    2504 c         call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
    2505          call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
    2506          call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
    2507          call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
    2508          call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
    2509          call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
    2510          call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
    2511          call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
    2512          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    2513          call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
    2514          call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
    2515          call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
    2516          call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
    2517          call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    2518          call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
    2519          call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
    2520          call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
    2521          call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
    2522          call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
    2523          call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
    2524          call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
    2525          call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
    2526          call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
    2527          call wrgradsfi(1,nlay,w_est,'w_est      ','w_est      ')
    2528 con sort les moments
    2529          call wrgradsfi(1,nlay,thetath2,'zh2       ','zh2       ')
    2530          call wrgradsfi(1,nlay,wth2,'w2       ','w2       ')
    2531          call wrgradsfi(1,nlay,wth3,'w3       ','w3       ')
    2532          call wrgradsfi(1,nlay,q2,'q2       ','q2       ')
    2533          call wrgradsfi(1,nlay,dtheta,'dT       ','dT       ')
    2534 c
    2535          call wrgradsfi(1,nlay,zw_sec,'zw_sec       ','zw_sec       ')
    2536          call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
    2537          call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
    2538          call wrgradsfi(1,nlay,nu,'nu       ','nu       ')
    2539 
    2540          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    2541          call wrgradsfi(1,nlay,zoa,'zoa        ','zoa        ')
    2542 c         call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
    2543 c         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
    2544 
    2545 cAM:nouveaux diagnostiques
    2546          call wrgradsfi(1,nlay,zthl,'zthl        ','zthl        ')
    2547          call wrgradsfi(1,nlay,zta,'zta        ','zta        ')
    2548          call wrgradsfi(1,nlay,zl,'zl        ','zl        ')
    2549          call wrgradsfi(1,nlay,zdthladj,'zdthladj    ',
    2550      s        'zdthladj    ')
    2551          call wrgradsfi(1,nlay,ztla,'ztla      ','ztla      ')
    2552          call wrgradsfi(1,nlay,zqta,'zqta      ','zqta      ')
    2553          call wrgradsfi(1,nlay,zqla,'zqla      ','zqla      ')
    2554          call wrgradsfi(1,nlay,deltaz,'deltaz      ','deltaz      ')
    2555 cCR:nouveaux diagnostiques
    2556       call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')
    2557       call wrgradsfi(1,nlay,detr_star  ,'detr_star   ','detr_star   ')     
    2558       call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
    2559       call wrgradsfi(1,nlay,zqsat    ,'zqsat   ','zqsat   ')
    2560       call wrgradsfi(1,nlay,zqsatth    ,'qsath   ','qsath   ')
    2561       call wrgradsfi(1,nlay,alim_star    ,'alim_star   ','alim_star   ')
    2562       call wrgradsfi(1,nlay,alim    ,'alim   ','alim   ')
    2563       call wrgradsfi(1,1,f,'f      ','f      ')
    2564       call wrgradsfi(1,1,alim_star_tot,'a_s_t      ','a_s_t      ')
    2565       call wrgradsfi(1,1,alim_star2,'a_2      ','a_2      ')
    2566       call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    2567       call wrgradsfi(1,1,zmax_sec,'z_sec      ','z_sec      ')
    2568 c      call wrgradsfi(1,1,zmax_sec2,'zz_se      ','zz_se      ')
    2569       call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
    2570       call wrgradsfi(1,1,nivcon,'nivcon      ','nivcon      ')
    2571       call wrgradsfi(1,1,zcon,'zcon      ','zcon      ')
    2572       zsortie1d(:)=lmax(:)
    2573       call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
    2574       call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
    2575       call wrgradsfi(1,1,wmax_sec,'w_sec      ','w_sec      ')
    2576       zsortie1d(:)=lmix(:)
    2577       call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
    2578       zsortie1d(:)=lentr(:)
    2579       call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
    2580 
    2581 c      print*,'17 OK convect8'
    2582 
    2583          do k=1,klev/10
    2584             write(str2,'(i2.2)') k
    2585             str10='wa'//str2
    2586             do l=1,nlay
    2587                do ig=1,ngrid
    2588                   zsortie(ig,l)=wa(ig,k,l)
    2589                enddo
    2590             enddo   
    2591             CALL wrgradsfi(1,nlay,zsortie,str10,str10)
    2592             do l=1,nlay
    2593                do ig=1,ngrid
    2594                   zsortie(ig,l)=larg_part(ig,k,l)
    2595                enddo
    2596             enddo
    2597            str10='la'//str2
    2598             CALL wrgradsfi(1,nlay,zsortie,str10,str10)
    2599          enddo
    2600 
    2601 
    2602 c     print*,'18 OK convect8'
    2603 c      endif
    2604       print*,'Fin des wrgradsfi'
    2605 #endif
    26062351
    26072352      endif
    26082353
    2609 c     if(wa_moy(1,4).gt.1.e-10) stop
    2610 
    2611       print*,'19 OK convect8'
     2354!     print*,'19 OK convect8'
    26122355      return
    26132356      end
     2357
    26142358      SUBROUTINE thermcell_eau(ngrid,nlay,ptimestep
    26152359     s                  ,pplay,pplev,pphi
    … …  
    27102454
    27112455      real count_time
    2712       integer isplit,nsplit,ialt
    2713       parameter (nsplit=10)
    2714       data isplit/0/
    2715       save isplit
    2716 c$OMP THREADPRIVATE(isplit)
     2456      integer ialt
    27172457
    27182458      logical sorties
    … …  
    28612601c   ---------------------------------------------------
    28622602
    2863       print*,'0 OK convect8'
     2603!     print*,'0 OK convect8'
    28642604
    28652605      DO 1010 l=1,nlay
    … …  
    35683308c------------------------------------------------------------------
    35693309
    3570       if(sorties) then
     3310      return
     3311      end
     3312
     3313      SUBROUTINE thermcell(ngrid,nlay,ptimestep
     3314     s                  ,pplay,pplev,pphi
     3315     s                  ,pu,pv,pt,po
     3316     s                  ,pduadj,pdvadj,pdtadj,pdoadj
     3317     s                  ,fm0,entr0
     3318c    s                  ,pu_therm,pv_therm
     3319     s                  ,r_aspect,l_mix,w2di,tho)
     3320
     3321      USE dimphy
     3322      IMPLICIT NONE
     3323
     3324c=======================================================================
     3325c
     3326c   Calcul du transport verticale dans la couche limite en presence
     3327c   de "thermiques" explicitement representes
     3328c
     3329c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
     3330c
     3331c   le thermique est supposé homogène et dissipé par mélange avec
     3332c   son environnement. la longueur l_mix contrôle l'efficacité du
     3333c   mélange
     3334c
     3335c   Le calcul du transport des différentes espèces se fait en prenant
     3336c   en compte:
     3337c     1. un flux de masse montant
     3338c     2. un flux de masse descendant
     3339c     3. un entrainement
     3340c     4. un detrainement
     3341c
     3342c=======================================================================
     3343
     3344c-----------------------------------------------------------------------
     3345c   declarations:
     3346c   -------------
     3347
     3348#include "dimensions.h"
     3349cccc#include "dimphy.h"
     3350#include "YOMCST.h"
     3351
     3352c   arguments:
     3353c   ----------
     3354
     3355      INTEGER ngrid,nlay,w2di
     3356      REAL tho
     3357      real ptimestep,l_mix,r_aspect
     3358      REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
     3359      REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
     3360      REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
     3361      REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
     3362      REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
     3363      real pphi(ngrid,nlay)
     3364
     3365      integer idetr
     3366      save idetr
     3367      data idetr/3/
     3368c$OMP THREADPRIVATE(idetr)
     3369
     3370c   local:
     3371c   ------
     3372
     3373      INTEGER ig,k,l,lmaxa(klon),lmix(klon)
     3374      real zsortie1d(klon)
     3375c CR: on remplace lmax(klon,klev+1)
     3376      INTEGER lmax(klon),lmin(klon),lentr(klon)
     3377      real linter(klon)
     3378      real zmix(klon), fracazmix(klon)
     3379c RC
     3380      real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
     3381
     3382      real zlev(klon,klev+1),zlay(klon,klev)
     3383      REAL zh(klon,klev),zdhadj(klon,klev)
     3384      REAL ztv(klon,klev)
     3385      real zu(klon,klev),zv(klon,klev),zo(klon,klev)
     3386      REAL wh(klon,klev+1)
     3387      real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
     3388      real zla(klon,klev+1)
     3389      real zwa(klon,klev+1)
     3390      real zld(klon,klev+1)
     3391      real zwd(klon,klev+1)
     3392      real zsortie(klon,klev)
     3393      real zva(klon,klev)
     3394      real zua(klon,klev)
     3395      real zoa(klon,klev)
     3396
     3397      real zha(klon,klev)
     3398      real wa_moy(klon,klev+1)
     3399      real fraca(klon,klev+1)
     3400      real fracc(klon,klev+1)
     3401      real zf,zf2
     3402      real thetath2(klon,klev),wth2(klon,klev)
     3403!      common/comtherm/thetath2,wth2
     3404
     3405      real count_time
     3406      integer ialt
     3407
     3408      logical sorties
     3409      real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
     3410      real zpspsk(klon,klev)
     3411
     3412c     real wmax(klon,klev),wmaxa(klon)
     3413      real wmax(klon),wmaxa(klon)
     3414      real wa(klon,klev,klev+1)
     3415      real wd(klon,klev+1)
     3416      real larg_part(klon,klev,klev+1)
     3417      real fracd(klon,klev+1)
     3418      real xxx(klon,klev+1)
     3419      real larg_cons(klon,klev+1)
     3420      real larg_detr(klon,klev+1)
     3421      real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
     3422      real pu_therm(klon,klev),pv_therm(klon,klev)
     3423      real fm(klon,klev+1),entr(klon,klev)
     3424      real fmc(klon,klev+1)
     3425
     3426cCR:nouvelles variables
     3427      real f_star(klon,klev+1),entr_star(klon,klev)
     3428      real entr_star_tot(klon),entr_star2(klon)
     3429      real f(klon), f0(klon)
     3430      real zlevinter(klon)
     3431      logical first
     3432      data first /.false./
     3433      save first
     3434c$OMP THREADPRIVATE(first)
     3435cRC
     3436
     3437      character*2 str2
     3438      character*10 str10
     3439
     3440      character (len=20) :: modname='thermcell'
     3441      character (len=80) :: abort_message
     3442
     3443      LOGICAL vtest(klon),down
     3444
     3445      EXTERNAL SCOPY
     3446
     3447      integer ncorrec,ll
     3448      save ncorrec
     3449      data ncorrec/0/
     3450c$OMP THREADPRIVATE(ncorrec)
     3451     
     3452c
     3453c-----------------------------------------------------------------------
     3454c   initialisation:
     3455c   ---------------
     3456c
     3457       sorties=.true.
     3458      IF(ngrid.NE.klon) THEN
     3459         PRINT*
     3460         PRINT*,'STOP dans convadj'
     3461         PRINT*,'ngrid    =',ngrid
     3462         PRINT*,'klon  =',klon
     3463      ENDIF
     3464c
     3465c-----------------------------------------------------------------------
     3466c   incrementation eventuelle de tendances precedentes:
     3467c   ---------------------------------------------------
     3468
     3469!      print*,'0 OK convect8'
     3470
     3471      DO 1010 l=1,nlay
     3472         DO 1015 ig=1,ngrid
     3473            zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
     3474            zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
     3475            zu(ig,l)=pu(ig,l)
     3476            zv(ig,l)=pv(ig,l)
     3477            zo(ig,l)=po(ig,l)
     3478            ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
     34791015     CONTINUE
     34801010  CONTINUE
     3481
     3482!      print*,'1 OK convect8'
     3483c                       --------------------
     3484c
     3485c
     3486c                       + + + + + + + + + + +
     3487c
     3488c
     3489c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
     3490c  wh,wt,wo ...
     3491c
     3492c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
     3493c
     3494c
     3495c                       --------------------   zlev(1)
     3496c                       \\\\\\\\\\\\\\\\\\\\
     3497c
     3498c
     3499
     3500c-----------------------------------------------------------------------
     3501c   Calcul des altitudes des couches
     3502c-----------------------------------------------------------------------
     3503
     3504      do l=2,nlay
     3505         do ig=1,ngrid
     3506            zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
     3507         enddo
     3508      enddo
     3509      do ig=1,ngrid
     3510         zlev(ig,1)=0.
     3511         zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
     3512      enddo
    35713513      do l=1,nlay
    35723514         do ig=1,ngrid
    3573             zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
    3574             zld(ig,l)=fracd(ig,l)*zmax(ig)
    3575             if(1.-fracd(ig,l).gt.1.e-10)
    3576      s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
    3577          enddo
    3578       enddo
     3515            zlay(ig,l)=pphi(ig,l)/RG
     3516         enddo
     3517      enddo
     3518
     3519c      print*,'2 OK convect8'
     3520c-----------------------------------------------------------------------
     3521c   Calcul des densites
     3522c-----------------------------------------------------------------------
    35793523
    35803524      do l=1,nlay
     3525         do ig=1,ngrid
     3526            rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
     3527         enddo
     3528      enddo
     3529
     3530      do l=2,nlay
     3531         do ig=1,ngrid
     3532            rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
     3533         enddo
     3534      enddo
     3535
     3536      do k=1,nlay
     3537         do l=1,nlay+1
     3538            do ig=1,ngrid
     3539               wa(ig,k,l)=0.
     3540            enddo
     3541         enddo
     3542      enddo
     3543
     3544c      print*,'3 OK convect8'
     3545c------------------------------------------------------------------
     3546c   Calcul de w2, quarre de w a partir de la cape
     3547c   a partir de w2, on calcule wa, vitesse de l'ascendance
     3548c
     3549c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
     3550c   w2 est stoke dans wa
     3551c
     3552c   ATTENTION: dans convect8, on n'utilise le calcule des wa
     3553c   independants par couches que pour calculer l'entrainement
     3554c   a la base et la hauteur max de l'ascendance.
     3555c
     3556c   Indicages:
     3557c   l'ascendance provenant du niveau k traverse l'interface l avec
     3558c   une vitesse wa(k,l).
     3559c
     3560c                       --------------------
     3561c
     3562c                       + + + + + + + + + +
     3563c
     3564c  wa(k,l)   ----       --------------------    l
     3565c             /\
     3566c            /||\       + + + + + + + + + +
     3567c             ||
     3568c             ||        --------------------
     3569c             ||
     3570c             ||        + + + + + + + + + +
     3571c             ||
     3572c             ||        --------------------
     3573c             ||__
     3574c             |___      + + + + + + + + + +     k
     3575c
     3576c                       --------------------
     3577c
     3578c
     3579c
     3580c------------------------------------------------------------------
     3581
     3582cCR: ponderation entrainement des couches instables
     3583cdef des entr_star tels que entr=f*entr_star     
     3584      do l=1,klev
     3585         do ig=1,ngrid
     3586            entr_star(ig,l)=0.
     3587         enddo
     3588      enddo
     3589c determination de la longueur de la couche d entrainement
     3590      do ig=1,ngrid
     3591         lentr(ig)=1
     3592      enddo
     3593
     3594con ne considere que les premieres couches instables
     3595      do k=nlay-2,1,-1
     3596         do ig=1,ngrid
     3597            if (ztv(ig,k).gt.ztv(ig,k+1).and.
     3598     s          ztv(ig,k+1).le.ztv(ig,k+2)) then
     3599               lentr(ig)=k
     3600            endif
     3601          enddo
     3602      enddo
     3603   
     3604c determination du lmin: couche d ou provient le thermique
     3605      do ig=1,ngrid
     3606         lmin(ig)=1
     3607      enddo
     3608      do ig=1,ngrid
     3609         do l=nlay,2,-1
     3610            if (ztv(ig,l-1).gt.ztv(ig,l)) then
     3611               lmin(ig)=l-1
     3612            endif
     3613         enddo
     3614      enddo
     3615c
     3616c definition de l'entrainement des couches
     3617      do l=1,klev-1
     3618         do ig=1,ngrid
     3619            if (ztv(ig,l).gt.ztv(ig,l+1).and.
     3620     s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
     3621                 entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
     3622     s                           (zlev(ig,l+1)-zlev(ig,l))
     3623            endif
     3624         enddo
     3625      enddo
     3626c pas de thermique si couches 1->5 stables
     3627      do ig=1,ngrid
     3628         if (lmin(ig).gt.5) then
     3629            do l=1,klev
     3630               entr_star(ig,l)=0.
     3631            enddo
     3632         endif
     3633      enddo
     3634c calcul de l entrainement total
     3635      do ig=1,ngrid
     3636         entr_star_tot(ig)=0.
     3637      enddo
     3638      do ig=1,ngrid
     3639         do k=1,klev
     3640            entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
     3641         enddo
     3642      enddo
     3643c
     3644      print*,'fin calcul entr_star'
     3645      do k=1,klev
     3646         do ig=1,ngrid
     3647            ztva(ig,k)=ztv(ig,k)
     3648         enddo
     3649      enddo
     3650cRC
     3651c      print*,'7 OK convect8'
     3652      do k=1,klev+1
     3653         do ig=1,ngrid
     3654            zw2(ig,k)=0.
     3655            fmc(ig,k)=0.
     3656cCR
     3657            f_star(ig,k)=0.
     3658cRC
     3659            larg_cons(ig,k)=0.
     3660            larg_detr(ig,k)=0.
     3661            wa_moy(ig,k)=0.
     3662         enddo
     3663      enddo
     3664
     3665c      print*,'8 OK convect8'
     3666      do ig=1,ngrid
     3667         linter(ig)=1.
     3668         lmaxa(ig)=1
     3669         lmix(ig)=1
     3670         wmaxa(ig)=0.
     3671      enddo
     3672
     3673cCR:
     3674      do l=1,nlay-2
     3675         do ig=1,ngrid
     3676            if (ztv(ig,l).gt.ztv(ig,l+1)
     3677     s         .and.entr_star(ig,l).gt.1.e-10
     3678     s         .and.zw2(ig,l).lt.1e-10) then
     3679               f_star(ig,l+1)=entr_star(ig,l)
     3680ctest:calcul de dteta
     3681               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
     3682     s                     *(zlev(ig,l+1)-zlev(ig,l))
     3683     s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
     3684               larg_detr(ig,l)=0.
     3685            else if ((zw2(ig,l).ge.1e-10).and.
     3686     s               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
     3687               f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
     3688               ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
     3689     s                    *ztv(ig,l))/f_star(ig,l+1)
     3690               zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
     3691     s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     3692     s                     *(zlev(ig,l+1)-zlev(ig,l))
     3693            endif
     3694c determination de zmax continu par interpolation lineaire
     3695            if (zw2(ig,l+1).lt.0.) then
     3696ctest
     3697               if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
     3698                  print*,'pb linter'
     3699               endif
     3700               linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
     3701     s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
     3702               zw2(ig,l+1)=0.
     3703               lmaxa(ig)=l
     3704            else
     3705               if (zw2(ig,l+1).lt.0.) then
     3706                  print*,'pb1 zw2<0'
     3707               endif
     3708               wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
     3709            endif
     3710            if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
     3711c   lmix est le niveau de la couche ou w (wa_moy) est maximum
     3712               lmix(ig)=l+1
     3713               wmaxa(ig)=wa_moy(ig,l+1)
     3714            endif
     3715         enddo
     3716      enddo
     3717      print*,'fin calcul zw2'
     3718c
     3719c Calcul de la couche correspondant a la hauteur du thermique
     3720      do ig=1,ngrid
     3721         lmax(ig)=lentr(ig)
     3722      enddo
     3723      do ig=1,ngrid
     3724         do l=nlay,lentr(ig)+1,-1
     3725            if (zw2(ig,l).le.1.e-10) then
     3726               lmax(ig)=l-1
     3727            endif
     3728         enddo
     3729      enddo
     3730c pas de thermique si couches 1->5 stables
     3731      do ig=1,ngrid
     3732         if (lmin(ig).gt.5) then
     3733            lmax(ig)=1
     3734            lmin(ig)=1
     3735         endif
     3736      enddo
     3737c   
     3738c Determination de zw2 max
     3739      do ig=1,ngrid
     3740         wmax(ig)=0.
     3741      enddo
     3742
     3743      do l=1,nlay
     3744         do ig=1,ngrid
     3745            if (l.le.lmax(ig)) then
     3746                if (zw2(ig,l).lt.0.)then
     3747                  print*,'pb2 zw2<0'
     3748                endif
     3749                zw2(ig,l)=sqrt(zw2(ig,l))
     3750                wmax(ig)=max(wmax(ig),zw2(ig,l))
     3751            else
     3752                 zw2(ig,l)=0.
     3753            endif
     3754          enddo
     3755      enddo
     3756
     3757c   Longueur caracteristique correspondant a la hauteur des thermiques.
     3758      do  ig=1,ngrid
     3759         zmax(ig)=0.
     3760         zlevinter(ig)=zlev(ig,1)
     3761      enddo
     3762      do  ig=1,ngrid
     3763c calcul de zlevinter
     3764          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
     3765     s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
     3766     s    -zlev(ig,lmax(ig)))
     3767       zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
     3768      enddo
     3769
     3770      print*,'avant fermeture'
     3771c Fermeture,determination de f
     3772      do ig=1,ngrid
     3773         entr_star2(ig)=0.
     3774      enddo
     3775      do ig=1,ngrid
     3776         if (entr_star_tot(ig).LT.1.e-10) then
     3777            f(ig)=0.
     3778         else
     3779             do k=lmin(ig),lentr(ig)
     3780                entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
     3781     s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
     3782             enddo
     3783c Nouvelle fermeture
     3784             f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
     3785     s             *entr_star2(ig))*entr_star_tot(ig)
     3786ctest
     3787c             if (first) then
     3788c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
     3789c     s             *wmax(ig))
     3790c             endif
     3791         endif
     3792c         f0(ig)=f(ig)
     3793c         first=.true.
     3794      enddo
     3795      print*,'apres fermeture'
     3796
     3797c Calcul de l'entrainement
     3798       do k=1,klev
     3799         do ig=1,ngrid
     3800            entr(ig,k)=f(ig)*entr_star(ig,k)
     3801         enddo
     3802      enddo
     3803c Calcul des flux
     3804      do ig=1,ngrid
     3805         do l=1,lmax(ig)-1
     3806            fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
     3807         enddo
     3808      enddo
     3809
     3810cRC
     3811
     3812
     3813c      print*,'9 OK convect8'
     3814c     print*,'WA1 ',wa_moy
     3815
     3816c   determination de l'indice du debut de la mixed layer ou w decroit
     3817
     3818c   calcul de la largeur de chaque ascendance dans le cas conservatif.
     3819c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
     3820c   d'une couche est égale à la hauteur de la couche alimentante.
     3821c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
     3822c   de la vitesse d'entrainement horizontal dans la couche alimentante.
     3823
     3824      do l=2,nlay
     3825         do ig=1,ngrid
     3826            if (l.le.lmaxa(ig)) then
     3827               zw=max(wa_moy(ig,l),1.e-10)
     3828               larg_cons(ig,l)=zmax(ig)*r_aspect
     3829     s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
     3830            endif
     3831         enddo
     3832      enddo
     3833
     3834      do l=2,nlay
     3835         do ig=1,ngrid
     3836            if (l.le.lmaxa(ig)) then
     3837c              if (idetr.eq.0) then
     3838c  cette option est finalement en dur.
     3839                  if ((l_mix*zlev(ig,l)).lt.0.)then
     3840                   print*,'pb l_mix*zlev<0'
     3841                  endif
     3842                  larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
     3843c              else if (idetr.eq.1) then
     3844c                 larg_detr(ig,l)=larg_cons(ig,l)
     3845c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
     3846c              else if (idetr.eq.2) then
     3847c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
     3848c    s            *sqrt(wa_moy(ig,l))
     3849c              else if (idetr.eq.4) then
     3850c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
     3851c    s            *wa_moy(ig,l)
     3852c              endif
     3853            endif
     3854         enddo
     3855       enddo
     3856
     3857c      print*,'10 OK convect8'
     3858c     print*,'WA2 ',wa_moy
     3859c   calcul de la fraction de la maille concernée par l'ascendance en tenant
     3860c   compte de l'epluchage du thermique.
     3861c
     3862cCR def de  zmix continu (profil parabolique des vitesses)
     3863      do ig=1,ngrid
     3864           if (lmix(ig).gt.1.) then
     3865c test
     3866              if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
     3867     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
     3868     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
     3869     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
     3870     s        then
     3871c             
     3872            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
     3873     s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
     3874     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
     3875     s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
     3876     s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
     3877     s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
     3878     s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
     3879     s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
     3880            else
     3881            zmix(ig)=zlev(ig,lmix(ig))
     3882            print*,'pb zmix'
     3883            endif
     3884         else
     3885         zmix(ig)=0.
     3886         endif
     3887ctest
     3888         if ((zmax(ig)-zmix(ig)).lt.0.) then
     3889            zmix(ig)=0.99*zmax(ig)
     3890c            print*,'pb zmix>zmax'
     3891         endif
     3892      enddo
     3893c
     3894c calcul du nouveau lmix correspondant
     3895      do ig=1,ngrid
     3896         do l=1,klev
     3897            if (zmix(ig).ge.zlev(ig,l).and.
     3898     s          zmix(ig).lt.zlev(ig,l+1)) then
     3899              lmix(ig)=l
     3900             endif
     3901          enddo
     3902      enddo
     3903c
     3904      do l=2,nlay
     3905         do ig=1,ngrid
     3906            if(larg_cons(ig,l).gt.1.) then
     3907c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
     3908               fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
     3909     s            /(r_aspect*zmax(ig))
     3910c test
     3911               fraca(ig,l)=max(fraca(ig,l),0.)
     3912               fraca(ig,l)=min(fraca(ig,l),0.5)
     3913               fracd(ig,l)=1.-fraca(ig,l)
     3914               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
     3915            else
     3916c              wa_moy(ig,l)=0.
     3917               fraca(ig,l)=0.
     3918               fracc(ig,l)=0.
     3919               fracd(ig,l)=1.
     3920            endif
     3921         enddo
     3922      enddo                 
     3923cCR: calcul de fracazmix
     3924       do ig=1,ngrid
     3925          fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
     3926     s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
     3927     s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
     3928     s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
     3929       enddo
     3930c
     3931       do l=2,nlay
     3932          do ig=1,ngrid
     3933             if(larg_cons(ig,l).gt.1.) then
     3934               if (l.gt.lmix(ig)) then
     3935ctest
     3936                 if (zmax(ig)-zmix(ig).lt.1.e-10) then
     3937c                   print*,'pb xxx'
     3938                   xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
     3939                 else
     3940                 xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
     3941                 endif
     3942           if (idetr.eq.0) then
     3943               fraca(ig,l)=fracazmix(ig)
     3944           else if (idetr.eq.1) then
     3945               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
     3946           else if (idetr.eq.2) then
     3947               fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
     3948           else
     3949               fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
     3950           endif
     3951c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
     3952               fraca(ig,l)=max(fraca(ig,l),0.)
     3953               fraca(ig,l)=min(fraca(ig,l),0.5)
     3954               fracd(ig,l)=1.-fraca(ig,l)
     3955               fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
     3956             endif
     3957            endif
     3958         enddo
     3959      enddo
     3960     
     3961      print*,'fin calcul fraca'
     3962c      print*,'11 OK convect8'
     3963c     print*,'Ea3 ',wa_moy
     3964c------------------------------------------------------------------
     3965c   Calcul de fracd, wd
     3966c   somme wa - wd = 0
     3967c------------------------------------------------------------------
     3968
     3969
     3970      do ig=1,ngrid
     3971         fm(ig,1)=0.
     3972         fm(ig,nlay+1)=0.
     3973      enddo
     3974
     3975      do l=2,nlay
     3976           do ig=1,ngrid
     3977              fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
     3978cCR:test
     3979              if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
     3980     s            .and.l.gt.lmix(ig)) then
     3981                 fm(ig,l)=fm(ig,l-1)
     3982c                 write(1,*)'ajustement fm, l',l
     3983              endif
     3984c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
     3985cRC
     3986           enddo
     3987         do ig=1,ngrid
     3988            if(fracd(ig,l).lt.0.1) then
     3989              abort_message = 'fracd trop petit'
     3990              CALL abort_gcm (modname,abort_message,1)
     3991            else
     3992c    vitesse descendante "diagnostique"
     3993               wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
     3994            endif
     3995         enddo
     3996      enddo
     3997
     3998      do l=1,nlay
     3999         do ig=1,ngrid
     4000c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
     4001            masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
     4002         enddo
     4003      enddo
     4004
     4005!     print*,'12 OK convect8'
     4006c     print*,'WA4 ',wa_moy
     4007cc------------------------------------------------------------------
     4008c   calcul du transport vertical
     4009c------------------------------------------------------------------
     4010
     4011      go to 4444
     4012c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
     4013      do l=2,nlay-1
     4014         do ig=1,ngrid
     4015            if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
     4016     s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
     4017c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
     4018c    s         ,fm(ig,l+1)*ptimestep
     4019c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
     4020             endif
     4021         enddo
     4022      enddo
     4023
     4024      do l=1,nlay
     4025         do ig=1,ngrid
     4026            if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
     4027c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
     4028c    s         ,entr(ig,l)*ptimestep
     4029c    s         ,'   M=',masse(ig,l)
     4030             endif
     4031         enddo
     4032      enddo
     4033
     4034      do l=1,nlay
     4035         do ig=1,ngrid
     4036            if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
     4037c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
     4038c    s         ,'   FM=',fm(ig,l)
     4039            endif
     4040            if(.not.masse(ig,l).ge.1.e-10
     4041     s         .or..not.masse(ig,l).le.1.e4) then
     4042c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
     4043c    s         ,'   M=',masse(ig,l)
     4044c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
     4045c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
     4046c     print*,'zlev(ig,l+1),zlev(ig,l)'
     4047c    s                ,zlev(ig,l+1),zlev(ig,l)
     4048c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
     4049c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
     4050            endif
     4051            if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
     4052c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
     4053c    s         ,'   E=',entr(ig,l)
     4054            endif
     4055         enddo
     4056      enddo
     4057
     40584444   continue
     4059
     4060cCR:redefinition du entr
     4061       do l=1,nlay
    35814062         do ig=1,ngrid
    35824063            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
    … …  
    35884069         enddo
    35894070      enddo
     4071cRC
     4072      if (w2di.eq.1) then
     4073         fm0=fm0+ptimestep*(fm-fm0)/tho
     4074         entr0=entr0+ptimestep*(entr-entr0)/tho
     4075      else
     4076         fm0=fm
     4077         entr0=entr
     4078      endif
     4079
     4080      if (1.eq.1) then
     4081         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
     4082     .    ,zh,zdhadj,zha)
     4083         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
     4084     .    ,zo,pdoadj,zoa)
     4085      else
     4086         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
     4087     .    ,zh,zdhadj,zha)
     4088         call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
     4089     .    ,zo,pdoadj,zoa)
     4090      endif
     4091
     4092      if (1.eq.0) then
     4093         call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
     4094     .    ,fraca,zmax
     4095     .    ,zu,zv,pduadj,pdvadj,zua,zva)
     4096      else
     4097         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
     4098     .    ,zu,pduadj,zua)
     4099         call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
     4100     .    ,zv,pdvadj,zva)
     4101      endif
     4102
     4103      do l=1,nlay
     4104         do ig=1,ngrid
     4105            zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
     4106            zf2=zf/(1.-zf)
     4107            thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
     4108            wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
     4109         enddo
     4110      enddo
     4111
     4112
     4113
     4114c     print*,'13 OK convect8'
     4115c     print*,'WA5 ',wa_moy
     4116      do l=1,nlay
     4117         do ig=1,ngrid
     4118            pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
     4119         enddo
     4120      enddo
     4121
     4122
     4123c     do l=1,nlay
     4124c        do ig=1,ngrid
     4125c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
     4126c     print*,'WARN!!! ig=',ig,'  l=',l
     4127c    s         ,'   pdtadj=',pdtadj(ig,l)
     4128c           endif
     4129c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
     4130c     print*,'WARN!!! ig=',ig,'  l=',l
     4131c    s         ,'   pdoadj=',pdoadj(ig,l)
     4132c           endif
     4133c        enddo
     4134c      enddo
     4135
     4136!     print*,'14 OK convect8'
     4137c------------------------------------------------------------------
     4138c   Calculs pour les sorties
     4139c------------------------------------------------------------------
     4140
     4141      if(sorties) then
     4142      do l=1,nlay
     4143         do ig=1,ngrid
     4144            zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
     4145            zld(ig,l)=fracd(ig,l)*zmax(ig)
     4146            if(1.-fracd(ig,l).gt.1.e-10)
     4147     s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
     4148         enddo
     4149      enddo
     4150
     4151cdeja fait
     4152c      do l=1,nlay
     4153c         do ig=1,ngrid
     4154c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
     4155c            if (detr(ig,l).lt.0.) then
     4156c                entr(ig,l)=entr(ig,l)-detr(ig,l)
     4157c                detr(ig,l)=0.
     4158c     print*,'WARNING !!! detrainement negatif ',ig,l
     4159c            endif
     4160c         enddo
     4161c      enddo
    35904162
    35914163c     print*,'15 OK convect8'
    35924164
    3593       isplit=isplit+1
    3594 
    35954165
    35964166c #define und
    3597         goto 123
     4167      goto 123
    35984168#ifdef und
    35994169      CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
    … …  
    36324202#endif
    36334203123   continue
    3634 ! #define troisD
     4204#define troisD
    36354205#ifdef troisD
    36364206c       if (sorties) then
    … …  
    36764246         call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
    36774247
    3678 cAM:nouveaux diagnostiques
    3679          call wrgradsfi(1,nlay,zthl,'zthl        ','zthl        ')
    3680          call wrgradsfi(1,nlay,zta,'zta        ','zta        ')
    3681          call wrgradsfi(1,nlay,zl,'zl        ','zl        ')
    3682          call wrgradsfi(1,nlay,zdthladj,'zdthladj    ',
    3683      s        'zdthladj    ')
    3684          call wrgradsfi(1,nlay,ztla,'ztla      ','ztla      ')
    3685          call wrgradsfi(1,nlay,zqta,'zqta      ','zqta      ')
    3686          call wrgradsfi(1,nlay,zqla,'zqla      ','zqla      ')
    36874248cCR:nouveaux diagnostiques
    36884249      call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
    … …  
    37284289c     if(wa_moy(1,4).gt.1.e-10) stop
    37294290
    3730 c     print*,'19 OK convect8'
    3731       return
    3732       end
    3733 
    3734       SUBROUTINE thermcell(ngrid,nlay,ptimestep
    3735      s                  ,pplay,pplev,pphi
    3736      s                  ,pu,pv,pt,po
    3737      s                  ,pduadj,pdvadj,pdtadj,pdoadj
    3738      s                  ,fm0,entr0
    3739 c    s                  ,pu_therm,pv_therm
    3740      s                  ,r_aspect,l_mix,w2di,tho)
    3741 
    3742       USE dimphy
    3743       IMPLICIT NONE
    3744 
    3745 c=======================================================================
    3746 c
    3747 c   Calcul du transport verticale dans la couche limite en presence
    3748 c   de "thermiques" explicitement representes
    3749 c
    3750 c   Réécriture à partir d'un listing papier à Habas, le 14/02/00
    3751 c
    3752 c   le thermique est supposé homogène et dissipé par mélange avec
    3753 c   son environnement. la longueur l_mix contrôle l'efficacité du
    3754 c   mélange
    3755 c
    3756 c   Le calcul du transport des différentes espèces se fait en prenant
    3757 c   en compte:
    3758 c     1. un flux de masse montant
    3759 c     2. un flux de masse descendant
    3760 c     3. un entrainement
    3761 c     4. un detrainement
    3762 c
    3763 c=======================================================================
    3764 
    3765 c-----------------------------------------------------------------------
    3766 c   declarations:
    3767 c   -------------
    3768 
    3769 #include "dimensions.h"
    3770 cccc#include "dimphy.h"
    3771 #include "YOMCST.h"
    3772 
    3773 c   arguments:
    3774 c   ----------
    3775 
    3776       INTEGER ngrid,nlay,w2di
    3777       REAL tho
    3778       real ptimestep,l_mix,r_aspect
    3779       REAL pt(ngrid,nlay),pdtadj(ngrid,nlay)
    3780       REAL pu(ngrid,nlay),pduadj(ngrid,nlay)
    3781       REAL pv(ngrid,nlay),pdvadj(ngrid,nlay)
    3782       REAL po(ngrid,nlay),pdoadj(ngrid,nlay)
    3783       REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    3784       real pphi(ngrid,nlay)
    3785 
    3786       integer idetr
    3787       save idetr
    3788       data idetr/3/
    3789 c$OMP THREADPRIVATE(idetr)
    3790 
    3791 c   local:
    3792 c   ------
    3793 
    3794       INTEGER ig,k,l,lmaxa(klon),lmix(klon)
    3795       real zsortie1d(klon)
    3796 c CR: on remplace lmax(klon,klev+1)
    3797       INTEGER lmax(klon),lmin(klon),lentr(klon)
    3798       real linter(klon)
    3799       real zmix(klon), fracazmix(klon)
    3800 c RC
    3801       real zmax(klon),zw,zz,zw2(klon,klev+1),ztva(klon,klev),zzz
    3802 
    3803       real zlev(klon,klev+1),zlay(klon,klev)
    3804       REAL zh(klon,klev),zdhadj(klon,klev)
    3805       REAL ztv(klon,klev)
    3806       real zu(klon,klev),zv(klon,klev),zo(klon,klev)
    3807       REAL wh(klon,klev+1)
    3808       real wu(klon,klev+1),wv(klon,klev+1),wo(klon,klev+1)
    3809       real zla(klon,klev+1)
    3810       real zwa(klon,klev+1)
    3811       real zld(klon,klev+1)
    3812       real zwd(klon,klev+1)
    3813       real zsortie(klon,klev)
    3814       real zva(klon,klev)
    3815       real zua(klon,klev)
    3816       real zoa(klon,klev)
    3817 
    3818       real zha(klon,klev)
    3819       real wa_moy(klon,klev+1)
    3820       real fraca(klon,klev+1)
    3821       real fracc(klon,klev+1)
    3822       real zf,zf2
    3823       real thetath2(klon,klev),wth2(klon,klev)
    3824 !      common/comtherm/thetath2,wth2
    3825 
    3826       real count_time
    3827       integer isplit,nsplit,ialt
    3828       parameter (nsplit=10)
    3829       data isplit/0/
    3830       save isplit
    3831 c$OMP THREADPRIVATE(isplit)
    3832 
    3833       logical sorties
    3834       real rho(klon,klev),rhobarz(klon,klev+1),masse(klon,klev)
    3835       real zpspsk(klon,klev)
    3836 
    3837 c     real wmax(klon,klev),wmaxa(klon)
    3838       real wmax(klon),wmaxa(klon)
    3839       real wa(klon,klev,klev+1)
    3840       real wd(klon,klev+1)
    3841       real larg_part(klon,klev,klev+1)
    3842       real fracd(klon,klev+1)
    3843       real xxx(klon,klev+1)
    3844       real larg_cons(klon,klev+1)
    3845       real larg_detr(klon,klev+1)
    3846       real fm0(klon,klev+1),entr0(klon,klev),detr(klon,klev)
    3847       real pu_therm(klon,klev),pv_therm(klon,klev)
    3848       real fm(klon,klev+1),entr(klon,klev)
    3849       real fmc(klon,klev+1)
    3850 
    3851 cCR:nouvelles variables
    3852       real f_star(klon,klev+1),entr_star(klon,klev)
    3853       real entr_star_tot(klon),entr_star2(klon)
    3854       real f(klon), f0(klon)
    3855       real zlevinter(klon)
    3856       logical first
    3857       data first /.false./
    3858       save first
    3859 c$OMP THREADPRIVATE(first)
    3860 cRC
    3861 
    3862       character*2 str2
    3863       character*10 str10
    3864 
    3865       character (len=20) :: modname='thermcell'
    3866       character (len=80) :: abort_message
    3867 
    3868       LOGICAL vtest(klon),down
    3869 
    3870       EXTERNAL SCOPY
    3871 
    3872       integer ncorrec,ll
    3873       save ncorrec
    3874       data ncorrec/0/
    3875 c$OMP THREADPRIVATE(ncorrec)
    3876      
    3877 c
    3878 c-----------------------------------------------------------------------
    3879 c   initialisation:
    3880 c   ---------------
    3881 c
    3882        sorties=.true.
    3883       IF(ngrid.NE.klon) THEN
    3884          PRINT*
    3885          PRINT*,'STOP dans convadj'
    3886          PRINT*,'ngrid    =',ngrid
    3887          PRINT*,'klon  =',klon
    3888       ENDIF
    3889 c
    3890 c-----------------------------------------------------------------------
    3891 c   incrementation eventuelle de tendances precedentes:
    3892 c   ---------------------------------------------------
    3893 
    3894        print*,'0 OK convect8'
    3895 
    3896       DO 1010 l=1,nlay
    3897          DO 1015 ig=1,ngrid
    3898             zpspsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**RKAPPA
    3899             zh(ig,l)=pt(ig,l)/zpspsk(ig,l)
    3900             zu(ig,l)=pu(ig,l)
    3901             zv(ig,l)=pv(ig,l)
    3902             zo(ig,l)=po(ig,l)
    3903             ztv(ig,l)=zh(ig,l)*(1.+0.61*zo(ig,l))
    3904 1015     CONTINUE
    3905 1010  CONTINUE
    3906 
    3907        print*,'1 OK convect8'
    3908 c                       --------------------
    3909 c
    3910 c
    3911 c                       + + + + + + + + + + +
    3912 c
    3913 c
    3914 c  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
    3915 c  wh,wt,wo ...
    3916 c
    3917 c                       + + + + + + + + + + +  zh,zu,zv,zo,rho
    3918 c
    3919 c
    3920 c                       --------------------   zlev(1)
    3921 c                       \\\\\\\\\\\\\\\\\\\\
    3922 c
    3923 c
    3924 
    3925 c-----------------------------------------------------------------------
    3926 c   Calcul des altitudes des couches
    3927 c-----------------------------------------------------------------------
    3928 
    3929       do l=2,nlay
    3930          do ig=1,ngrid
    3931             zlev(ig,l)=0.5*(pphi(ig,l)+pphi(ig,l-1))/RG
    3932          enddo
    3933       enddo
    3934       do ig=1,ngrid
    3935          zlev(ig,1)=0.
    3936          zlev(ig,nlay+1)=(2.*pphi(ig,klev)-pphi(ig,klev-1))/RG
    3937       enddo
    3938       do l=1,nlay
    3939          do ig=1,ngrid
    3940             zlay(ig,l)=pphi(ig,l)/RG
    3941          enddo
    3942       enddo
    3943 
    3944 c      print*,'2 OK convect8'
    3945 c-----------------------------------------------------------------------
    3946 c   Calcul des densites
    3947 c-----------------------------------------------------------------------
    3948 
    3949       do l=1,nlay
    3950          do ig=1,ngrid
    3951             rho(ig,l)=pplay(ig,l)/(zpspsk(ig,l)*RD*zh(ig,l))
    3952          enddo
    3953       enddo
    3954 
    3955       do l=2,nlay
    3956          do ig=1,ngrid
    3957             rhobarz(ig,l)=0.5*(rho(ig,l)+rho(ig,l-1))
    3958          enddo
    3959       enddo
    3960 
    3961       do k=1,nlay
    3962          do l=1,nlay+1
    3963             do ig=1,ngrid
    3964                wa(ig,k,l)=0.
    3965             enddo
    3966          enddo
    3967       enddo
    3968 
    3969 c      print*,'3 OK convect8'
    3970 c------------------------------------------------------------------
    3971 c   Calcul de w2, quarre de w a partir de la cape
    3972 c   a partir de w2, on calcule wa, vitesse de l'ascendance
    3973 c
    3974 c   ATTENTION: Dans cette version, pour cause d'economie de memoire,
    3975 c   w2 est stoke dans wa
    3976 c
    3977 c   ATTENTION: dans convect8, on n'utilise le calcule des wa
    3978 c   independants par couches que pour calculer l'entrainement
    3979 c   a la base et la hauteur max de l'ascendance.
    3980 c
    3981 c   Indicages:
    3982 c   l'ascendance provenant du niveau k traverse l'interface l avec
    3983 c   une vitesse wa(k,l).
    3984 c
    3985 c                       --------------------
    3986 c
    3987 c                       + + + + + + + + + +
    3988 c
    3989 c  wa(k,l)   ----       --------------------    l
    3990 c             /\
    3991 c            /||\       + + + + + + + + + +
    3992 c             ||
    3993 c             ||        --------------------
    3994 c             ||
    3995 c             ||        + + + + + + + + + +
    3996 c             ||
    3997 c             ||        --------------------
    3998 c             ||__
    3999 c             |___      + + + + + + + + + +     k
    4000 c
    4001 c                       --------------------
    4002 c
    4003 c
    4004 c
    4005 c------------------------------------------------------------------
    4006 
    4007 cCR: ponderation entrainement des couches instables
    4008 cdef des entr_star tels que entr=f*entr_star     
    4009       do l=1,klev
    4010          do ig=1,ngrid
    4011             entr_star(ig,l)=0.
    4012          enddo
    4013       enddo
    4014 c determination de la longueur de la couche d entrainement
    4015       do ig=1,ngrid
    4016          lentr(ig)=1
    4017       enddo
    4018 
    4019 con ne considere que les premieres couches instables
    4020       do k=nlay-2,1,-1
    4021          do ig=1,ngrid
    4022             if (ztv(ig,k).gt.ztv(ig,k+1).and.
    4023      s          ztv(ig,k+1).le.ztv(ig,k+2)) then
    4024                lentr(ig)=k
    4025             endif
    4026           enddo
    4027       enddo
    4028    
    4029 c determination du lmin: couche d ou provient le thermique
    4030       do ig=1,ngrid
    4031          lmin(ig)=1
    4032       enddo
    4033       do ig=1,ngrid
    4034          do l=nlay,2,-1
    4035             if (ztv(ig,l-1).gt.ztv(ig,l)) then
    4036                lmin(ig)=l-1
    4037             endif
    4038          enddo
    4039       enddo
    4040 c
    4041 c definition de l'entrainement des couches
    4042       do l=1,klev-1
    4043          do ig=1,ngrid
    4044             if (ztv(ig,l).gt.ztv(ig,l+1).and.
    4045      s          l.ge.lmin(ig).and.l.le.lentr(ig)) then
    4046                  entr_star(ig,l)=(ztv(ig,l)-ztv(ig,l+1))*
    4047      s                           (zlev(ig,l+1)-zlev(ig,l))
    4048             endif
    4049          enddo
    4050       enddo
    4051 c pas de thermique si couches 1->5 stables
    4052       do ig=1,ngrid
    4053          if (lmin(ig).gt.5) then
    4054             do l=1,klev
    4055                entr_star(ig,l)=0.
    4056             enddo
    4057          endif
    4058       enddo
    4059 c calcul de l entrainement total
    4060       do ig=1,ngrid
    4061          entr_star_tot(ig)=0.
    4062       enddo
    4063       do ig=1,ngrid
    4064          do k=1,klev
    4065             entr_star_tot(ig)=entr_star_tot(ig)+entr_star(ig,k)
    4066          enddo
    4067       enddo
    4068 c
    4069       print*,'fin calcul entr_star'
    4070       do k=1,klev
    4071          do ig=1,ngrid
    4072             ztva(ig,k)=ztv(ig,k)
    4073          enddo
    4074       enddo
    4075 cRC
    4076 c      print*,'7 OK convect8'
    4077       do k=1,klev+1
    4078          do ig=1,ngrid
    4079             zw2(ig,k)=0.
    4080             fmc(ig,k)=0.
    4081 cCR
    4082             f_star(ig,k)=0.
    4083 cRC
    4084             larg_cons(ig,k)=0.
    4085             larg_detr(ig,k)=0.
    4086             wa_moy(ig,k)=0.
    4087          enddo
    4088       enddo
    4089 
    4090 c      print*,'8 OK convect8'
    4091       do ig=1,ngrid
    4092          linter(ig)=1.
    4093          lmaxa(ig)=1
    4094          lmix(ig)=1
    4095          wmaxa(ig)=0.
    4096       enddo
    4097 
    4098 cCR:
    4099       do l=1,nlay-2
    4100          do ig=1,ngrid
    4101             if (ztv(ig,l).gt.ztv(ig,l+1)
    4102      s         .and.entr_star(ig,l).gt.1.e-10
    4103      s         .and.zw2(ig,l).lt.1e-10) then
    4104                f_star(ig,l+1)=entr_star(ig,l)
    4105 ctest:calcul de dteta
    4106                zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)
    4107      s                     *(zlev(ig,l+1)-zlev(ig,l))
    4108      s                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
    4109                larg_detr(ig,l)=0.
    4110             else if ((zw2(ig,l).ge.1e-10).and.
    4111      s               (f_star(ig,l)+entr_star(ig,l).gt.1.e-10)) then
    4112                f_star(ig,l+1)=f_star(ig,l)+entr_star(ig,l)
    4113                ztva(ig,l)=(f_star(ig,l)*ztva(ig,l-1)+entr_star(ig,l)
    4114      s                    *ztv(ig,l))/f_star(ig,l+1)
    4115                zw2(ig,l+1)=zw2(ig,l)*(f_star(ig,l)/f_star(ig,l+1))**2+
    4116      s                     2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    4117      s                     *(zlev(ig,l+1)-zlev(ig,l))
    4118             endif
    4119 c determination de zmax continu par interpolation lineaire
    4120             if (zw2(ig,l+1).lt.0.) then
    4121 ctest
    4122                if (abs(zw2(ig,l+1)-zw2(ig,l)).lt.1e-10) then
    4123                   print*,'pb linter'
    4124                endif
    4125                linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))
    4126      s           -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    4127                zw2(ig,l+1)=0.
    4128                lmaxa(ig)=l
    4129             else
    4130                if (zw2(ig,l+1).lt.0.) then
    4131                   print*,'pb1 zw2<0'
    4132                endif
    4133                wa_moy(ig,l+1)=sqrt(zw2(ig,l+1))
    4134             endif
    4135             if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
    4136 c   lmix est le niveau de la couche ou w (wa_moy) est maximum
    4137                lmix(ig)=l+1
    4138                wmaxa(ig)=wa_moy(ig,l+1)
    4139             endif
    4140          enddo
    4141       enddo
    4142       print*,'fin calcul zw2'
    4143 c
    4144 c Calcul de la couche correspondant a la hauteur du thermique
    4145       do ig=1,ngrid
    4146          lmax(ig)=lentr(ig)
    4147       enddo
    4148       do ig=1,ngrid
    4149          do l=nlay,lentr(ig)+1,-1
    4150             if (zw2(ig,l).le.1.e-10) then
    4151                lmax(ig)=l-1
    4152             endif
    4153          enddo
    4154       enddo
    4155 c pas de thermique si couches 1->5 stables
    4156       do ig=1,ngrid
    4157          if (lmin(ig).gt.5) then
    4158             lmax(ig)=1
    4159             lmin(ig)=1
    4160          endif
    4161       enddo
    4162 c   
    4163 c Determination de zw2 max
    4164       do ig=1,ngrid
    4165          wmax(ig)=0.
    4166       enddo
    4167 
    4168       do l=1,nlay
    4169          do ig=1,ngrid
    4170             if (l.le.lmax(ig)) then
    4171                 if (zw2(ig,l).lt.0.)then
    4172                   print*,'pb2 zw2<0'
    4173                 endif
    4174                 zw2(ig,l)=sqrt(zw2(ig,l))
    4175                 wmax(ig)=max(wmax(ig),zw2(ig,l))
    4176             else
    4177                  zw2(ig,l)=0.
    4178             endif
    4179           enddo
    4180       enddo
    4181 
    4182 c   Longueur caracteristique correspondant a la hauteur des thermiques.
    4183       do  ig=1,ngrid
    4184          zmax(ig)=0.
    4185          zlevinter(ig)=zlev(ig,1)
    4186       enddo
    4187       do  ig=1,ngrid
    4188 c calcul de zlevinter
    4189           zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*
    4190      s    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)
    4191      s    -zlev(ig,lmax(ig)))
    4192        zmax(ig)=max(zmax(ig),zlevinter(ig)-zlev(ig,lmin(ig)))
    4193       enddo
    4194 
    4195       print*,'avant fermeture'
    4196 c Fermeture,determination de f
    4197       do ig=1,ngrid
    4198          entr_star2(ig)=0.
    4199       enddo
    4200       do ig=1,ngrid
    4201          if (entr_star_tot(ig).LT.1.e-10) then
    4202             f(ig)=0.
    4203          else
    4204              do k=lmin(ig),lentr(ig)
    4205                 entr_star2(ig)=entr_star2(ig)+entr_star(ig,k)**2
    4206      s                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
    4207              enddo
    4208 c Nouvelle fermeture
    4209              f(ig)=wmax(ig)/(max(500.,zmax(ig))*r_aspect
    4210      s             *entr_star2(ig))*entr_star_tot(ig)
    4211 ctest
    4212 c             if (first) then
    4213 c             f(ig)=f(ig)+(f0(ig)-f(ig))*exp(-ptimestep/zmax(ig)
    4214 c     s             *wmax(ig))
    4215 c             endif
    4216          endif
    4217 c         f0(ig)=f(ig)
    4218 c         first=.true.
    4219       enddo
    4220       print*,'apres fermeture'
    4221 
    4222 c Calcul de l'entrainement
    4223        do k=1,klev
    4224          do ig=1,ngrid
    4225             entr(ig,k)=f(ig)*entr_star(ig,k)
    4226          enddo
    4227       enddo
    4228 c Calcul des flux
    4229       do ig=1,ngrid
    4230          do l=1,lmax(ig)-1
    4231             fmc(ig,l+1)=fmc(ig,l)+entr(ig,l)
    4232          enddo
    4233       enddo
    4234 
    4235 cRC
    4236 
    4237 
    4238 c      print*,'9 OK convect8'
    4239 c     print*,'WA1 ',wa_moy
    4240 
    4241 c   determination de l'indice du debut de la mixed layer ou w decroit
    4242 
    4243 c   calcul de la largeur de chaque ascendance dans le cas conservatif.
    4244 c   dans ce cas simple, on suppose que la largeur de l'ascendance provenant
    4245 c   d'une couche est égale à la hauteur de la couche alimentante.
    4246 c   La vitesse maximale dans l'ascendance est aussi prise comme estimation
    4247 c   de la vitesse d'entrainement horizontal dans la couche alimentante.
    4248 
    4249       do l=2,nlay
    4250          do ig=1,ngrid
    4251             if (l.le.lmaxa(ig)) then
    4252                zw=max(wa_moy(ig,l),1.e-10)
    4253                larg_cons(ig,l)=zmax(ig)*r_aspect
    4254      s         *fmc(ig,l)/(rhobarz(ig,l)*zw)
    4255             endif
    4256          enddo
    4257       enddo
    4258 
    4259       do l=2,nlay
    4260          do ig=1,ngrid
    4261             if (l.le.lmaxa(ig)) then
    4262 c              if (idetr.eq.0) then
    4263 c  cette option est finalement en dur.
    4264                   if ((l_mix*zlev(ig,l)).lt.0.)then
    4265                    print*,'pb l_mix*zlev<0'
    4266                   endif
    4267                   larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
    4268 c              else if (idetr.eq.1) then
    4269 c                 larg_detr(ig,l)=larg_cons(ig,l)
    4270 c    s            *sqrt(l_mix*zlev(ig,l))/larg_cons(ig,lmix(ig))
    4271 c              else if (idetr.eq.2) then
    4272 c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
    4273 c    s            *sqrt(wa_moy(ig,l))
    4274 c              else if (idetr.eq.4) then
    4275 c                 larg_detr(ig,l)=sqrt(l_mix*zlev(ig,l))
    4276 c    s            *wa_moy(ig,l)
    4277 c              endif
    4278             endif
    4279          enddo
    4280        enddo
    4281 
    4282 c      print*,'10 OK convect8'
    4283 c     print*,'WA2 ',wa_moy
    4284 c   calcul de la fraction de la maille concernée par l'ascendance en tenant
    4285 c   compte de l'epluchage du thermique.
    4286 c
    4287 cCR def de  zmix continu (profil parabolique des vitesses)
    4288       do ig=1,ngrid
    4289            if (lmix(ig).gt.1.) then
    4290 c test
    4291               if (((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
    4292      s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
    4293      s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
    4294      s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))).gt.1e-10)
    4295      s        then
    4296 c             
    4297             zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
    4298      s        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)
    4299      s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
    4300      s        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))
    4301      s        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))
    4302      s        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))
    4303      s        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))
    4304      s        *((zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))))))
    4305             else
    4306             zmix(ig)=zlev(ig,lmix(ig))
    4307             print*,'pb zmix'
    4308             endif
    4309          else
    4310          zmix(ig)=0.
    4311          endif
    4312 ctest
    4313          if ((zmax(ig)-zmix(ig)).lt.0.) then
    4314             zmix(ig)=0.99*zmax(ig)
    4315 c            print*,'pb zmix>zmax'
    4316          endif
    4317       enddo
    4318 c
    4319 c calcul du nouveau lmix correspondant
    4320       do ig=1,ngrid
    4321          do l=1,klev
    4322             if (zmix(ig).ge.zlev(ig,l).and.
    4323      s          zmix(ig).lt.zlev(ig,l+1)) then
    4324               lmix(ig)=l
    4325              endif
    4326           enddo
    4327       enddo
    4328 c
    4329       do l=2,nlay
    4330          do ig=1,ngrid
    4331             if(larg_cons(ig,l).gt.1.) then
    4332 c     print*,ig,l,lmix(ig),lmaxa(ig),larg_cons(ig,l),'  KKK'
    4333                fraca(ig,l)=(larg_cons(ig,l)-larg_detr(ig,l))
    4334      s            /(r_aspect*zmax(ig))
    4335 c test
    4336                fraca(ig,l)=max(fraca(ig,l),0.)
    4337                fraca(ig,l)=min(fraca(ig,l),0.5)
    4338                fracd(ig,l)=1.-fraca(ig,l)
    4339                fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
    4340             else
    4341 c              wa_moy(ig,l)=0.
    4342                fraca(ig,l)=0.
    4343                fracc(ig,l)=0.
    4344                fracd(ig,l)=1.
    4345             endif
    4346          enddo
    4347       enddo                 
    4348 cCR: calcul de fracazmix
    4349        do ig=1,ngrid
    4350           fracazmix(ig)=(fraca(ig,lmix(ig)+1)-fraca(ig,lmix(ig)))/
    4351      s     (zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))*zmix(ig)
    4352      s    +fraca(ig,lmix(ig))-zlev(ig,lmix(ig))*(fraca(ig,lmix(ig)+1)
    4353      s    -fraca(ig,lmix(ig)))/(zlev(ig,lmix(ig)+1)-zlev(ig,lmix(ig)))
    4354        enddo
    4355 c
    4356        do l=2,nlay
    4357           do ig=1,ngrid
    4358              if(larg_cons(ig,l).gt.1.) then
    4359                if (l.gt.lmix(ig)) then
    4360 ctest
    4361                  if (zmax(ig)-zmix(ig).lt.1.e-10) then
    4362 c                   print*,'pb xxx'
    4363                    xxx(ig,l)=(lmaxa(ig)+1.-l)/(lmaxa(ig)+1.-lmix(ig))
    4364                  else
    4365                  xxx(ig,l)=(zmax(ig)-zlev(ig,l))/(zmax(ig)-zmix(ig))
    4366                  endif
    4367            if (idetr.eq.0) then
    4368                fraca(ig,l)=fracazmix(ig)
    4369            else if (idetr.eq.1) then
    4370                fraca(ig,l)=fracazmix(ig)*xxx(ig,l)
    4371            else if (idetr.eq.2) then
    4372                fraca(ig,l)=fracazmix(ig)*(1.-(1.-xxx(ig,l))**2)
    4373            else
    4374                fraca(ig,l)=fracazmix(ig)*xxx(ig,l)**2
    4375            endif
    4376 c     print*,ig,l,lmix(ig),lmaxa(ig),xxx(ig,l),'LLLLLLL'
    4377                fraca(ig,l)=max(fraca(ig,l),0.)
    4378                fraca(ig,l)=min(fraca(ig,l),0.5)
    4379                fracd(ig,l)=1.-fraca(ig,l)
    4380                fracc(ig,l)=larg_cons(ig,l)/(r_aspect*zmax(ig))
    4381              endif
    4382             endif
    4383          enddo
    4384       enddo
    4385      
    4386       print*,'fin calcul fraca'
    4387 c      print*,'11 OK convect8'
    4388 c     print*,'Ea3 ',wa_moy
    4389 c------------------------------------------------------------------
    4390 c   Calcul de fracd, wd
    4391 c   somme wa - wd = 0
    4392 c------------------------------------------------------------------
    4393 
    4394 
    4395       do ig=1,ngrid
    4396          fm(ig,1)=0.
    4397          fm(ig,nlay+1)=0.
    4398       enddo
    4399 
    4400       do l=2,nlay
    4401            do ig=1,ngrid
    4402               fm(ig,l)=fraca(ig,l)*wa_moy(ig,l)*rhobarz(ig,l)
    4403 cCR:test
    4404               if (entr(ig,l-1).lt.1e-10.and.fm(ig,l).gt.fm(ig,l-1)
    4405      s            .and.l.gt.lmix(ig)) then
    4406                  fm(ig,l)=fm(ig,l-1)
    4407 c                 write(1,*)'ajustement fm, l',l
    4408               endif
    4409 c              write(1,*)'ig,l,fm(ig,l)',ig,l,fm(ig,l)
    4410 cRC
    4411            enddo
    4412          do ig=1,ngrid
    4413             if(fracd(ig,l).lt.0.1) then
    4414               abort_message = 'fracd trop petit'
    4415               CALL abort_gcm (modname,abort_message,1)
    4416             else
    4417 c    vitesse descendante "diagnostique"
    4418                wd(ig,l)=fm(ig,l)/(fracd(ig,l)*rhobarz(ig,l))
    4419             endif
    4420          enddo
    4421       enddo
    4422 
    4423       do l=1,nlay
    4424          do ig=1,ngrid
    4425 c           masse(ig,l)=rho(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    4426             masse(ig,l)=(pplev(ig,l)-pplev(ig,l+1))/RG
    4427          enddo
    4428       enddo
    4429 
    4430       print*,'12 OK convect8'
    4431 c     print*,'WA4 ',wa_moy
    4432 cc------------------------------------------------------------------
    4433 c   calcul du transport vertical
    4434 c------------------------------------------------------------------
    4435 
    4436       go to 4444
    4437 c     print*,'XXXXXXXXXXXXXXX ptimestep= ',ptimestep
    4438       do l=2,nlay-1
    4439          do ig=1,ngrid
    4440             if(fm(ig,l+1)*ptimestep.gt.masse(ig,l)
    4441      s      .and.fm(ig,l+1)*ptimestep.gt.masse(ig,l+1)) then
    4442 c     print*,'WARN!!! FM>M ig=',ig,' l=',l,'  FM='
    4443 c    s         ,fm(ig,l+1)*ptimestep
    4444 c    s         ,'   M=',masse(ig,l),masse(ig,l+1)
    4445              endif
    4446          enddo
    4447       enddo
    4448 
    4449       do l=1,nlay
    4450          do ig=1,ngrid
    4451             if(entr(ig,l)*ptimestep.gt.masse(ig,l)) then
    4452 c     print*,'WARN!!! E>M ig=',ig,' l=',l,'  E=='
    4453 c    s         ,entr(ig,l)*ptimestep
    4454 c    s         ,'   M=',masse(ig,l)
    4455              endif
    4456          enddo
    4457       enddo
    4458 
    4459       do l=1,nlay
    4460          do ig=1,ngrid
    4461             if(.not.fm(ig,l).ge.0..or..not.fm(ig,l).le.10.) then
    4462 c     print*,'WARN!!! fm exagere ig=',ig,'   l=',l
    4463 c    s         ,'   FM=',fm(ig,l)
    4464             endif
    4465             if(.not.masse(ig,l).ge.1.e-10
    4466      s         .or..not.masse(ig,l).le.1.e4) then
    4467 c     print*,'WARN!!! masse exagere ig=',ig,'   l=',l
    4468 c    s         ,'   M=',masse(ig,l)
    4469 c     print*,'rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)',
    4470 c    s                 rho(ig,l),pplay(ig,l),zpspsk(ig,l),RD,zh(ig,l)
    4471 c     print*,'zlev(ig,l+1),zlev(ig,l)'
    4472 c    s                ,zlev(ig,l+1),zlev(ig,l)
    4473 c     print*,'pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)'
    4474 c    s                ,pphi(ig,l-1),pphi(ig,l),pphi(ig,l+1)
    4475             endif
    4476             if(.not.entr(ig,l).ge.0..or..not.entr(ig,l).le.10.) then
    4477 c     print*,'WARN!!! entr exagere ig=',ig,'   l=',l
    4478 c    s         ,'   E=',entr(ig,l)
    4479             endif
    4480          enddo
    4481       enddo
    4482 
    4483 4444   continue
    4484 
    4485 cCR:redefinition du entr
    4486        do l=1,nlay
    4487          do ig=1,ngrid
    4488             detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
    4489             if (detr(ig,l).lt.0.) then
    4490                 entr(ig,l)=entr(ig,l)-detr(ig,l)
    4491                 detr(ig,l)=0.
    4492 c     print*,'WARNING !!! detrainement negatif ',ig,l
    4493             endif
    4494          enddo
    4495       enddo
    4496 cRC
    4497       if (w2di.eq.1) then
    4498          fm0=fm0+ptimestep*(fm-fm0)/tho
    4499          entr0=entr0+ptimestep*(entr-entr0)/tho
    4500       else
    4501          fm0=fm
    4502          entr0=entr
    4503       endif
    4504 
    4505       if (1.eq.1) then
    4506          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    4507      .    ,zh,zdhadj,zha)
    4508          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    4509      .    ,zo,pdoadj,zoa)
    4510       else
    4511          call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
    4512      .    ,zh,zdhadj,zha)
    4513          call dqthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca
    4514      .    ,zo,pdoadj,zoa)
    4515       endif
    4516 
    4517       if (1.eq.0) then
    4518          call dvthermcell2(ngrid,nlay,ptimestep,fm0,entr0,masse
    4519      .    ,fraca,zmax
    4520      .    ,zu,zv,pduadj,pdvadj,zua,zva)
    4521       else
    4522          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    4523      .    ,zu,pduadj,zua)
    4524          call dqthermcell(ngrid,nlay,ptimestep,fm0,entr0,masse
    4525      .    ,zv,pdvadj,zva)
    4526       endif
    4527 
    4528       do l=1,nlay
    4529          do ig=1,ngrid
    4530             zf=0.5*(fracc(ig,l)+fracc(ig,l+1))
    4531             zf2=zf/(1.-zf)
    4532             thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l))**2
    4533             wth2(ig,l)=zf2*(0.5*(wa_moy(ig,l)+wa_moy(ig,l+1)))**2
    4534          enddo
    4535       enddo
    4536 
    4537 
    4538 
    4539 c     print*,'13 OK convect8'
    4540 c     print*,'WA5 ',wa_moy
    4541       do l=1,nlay
    4542          do ig=1,ngrid
    4543             pdtadj(ig,l)=zdhadj(ig,l)*zpspsk(ig,l)
    4544          enddo
    4545       enddo
    4546 
    4547 
    4548 c     do l=1,nlay
    4549 c        do ig=1,ngrid
    4550 c           if(abs(pdtadj(ig,l))*86400..gt.500.) then
    4551 c     print*,'WARN!!! ig=',ig,'  l=',l
    4552 c    s         ,'   pdtadj=',pdtadj(ig,l)
    4553 c           endif
    4554 c           if(abs(pdoadj(ig,l))*86400..gt.1.) then
    4555 c     print*,'WARN!!! ig=',ig,'  l=',l
    4556 c    s         ,'   pdoadj=',pdoadj(ig,l)
    4557 c           endif
    4558 c        enddo
    4559 c      enddo
    4560 
    4561       print*,'14 OK convect8'
    4562 c------------------------------------------------------------------
    4563 c   Calculs pour les sorties
    4564 c------------------------------------------------------------------
    4565 
    4566       if(sorties) then
    4567       do l=1,nlay
    4568          do ig=1,ngrid
    4569             zla(ig,l)=(1.-fracd(ig,l))*zmax(ig)
    4570             zld(ig,l)=fracd(ig,l)*zmax(ig)
    4571             if(1.-fracd(ig,l).gt.1.e-10)
    4572      s      zwa(ig,l)=wd(ig,l)*fracd(ig,l)/(1.-fracd(ig,l))
    4573          enddo
    4574       enddo
    4575 
    4576 cdeja fait
    4577 c      do l=1,nlay
    4578 c         do ig=1,ngrid
    4579 c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
    4580 c            if (detr(ig,l).lt.0.) then
    4581 c                entr(ig,l)=entr(ig,l)-detr(ig,l)
    4582 c                detr(ig,l)=0.
    4583 c     print*,'WARNING !!! detrainement negatif ',ig,l
    4584 c            endif
    4585 c         enddo
    4586 c      enddo
    4587 
    4588 c     print*,'15 OK convect8'
    4589 
    4590       isplit=isplit+1
    4591 
    4592 
    4593 c #define und
    4594         goto 123
    4595 #ifdef und
    4596       CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
    4597       CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
    4598       CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
    4599       CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
    4600       CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
    4601       CALL writeg1d(1,nlay,zla,'la      ','la      ')
    4602       CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
    4603       CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
    4604       CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
    4605       CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
    4606       CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
    4607       CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
    4608       CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
    4609       CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
    4610       CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
    4611       CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
    4612       CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
    4613       CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
    4614       CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
    4615       CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
    4616       CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
    4617       CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
    4618       CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
    4619       CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
    4620 
    4621       CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
    4622       CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
    4623       CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
    4624 
    4625 c   recalcul des flux en diagnostique...
    4626 c     print*,'PAS DE TEMPS ',ptimestep
    4627        call dt2F(pplev,pplay,pt,pdtadj,wh)
    4628       CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
    4629 #endif
    4630 123   continue
    4631 #define troisD
    4632 #ifdef troisD
    4633 c       if (sorties) then
    4634       print*,'Debut des wrgradsfi'
    4635 
    4636 c      print*,'16 OK convect8'
    4637          call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
    4638          call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
    4639          call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
    4640          call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
    4641          call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
    4642          call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
    4643 c      print*,'WA6 ',wa_moy
    4644          call wrgradsfi(1,nlay,zla,'la        ','la        ')
    4645          call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
    4646          call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
    4647          call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
    4648          call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
    4649          call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
    4650          call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
    4651          call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
    4652          call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
    4653          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    4654          call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
    4655          call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
    4656          call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
    4657          call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
    4658          call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    4659          call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
    4660          call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
    4661          call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
    4662          call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
    4663          call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
    4664          call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
    4665          call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
    4666          call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
    4667          call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
    4668          call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
    4669          call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
    4670 
    4671          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    4672          call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
    4673          call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
    4674 
    4675 cCR:nouveaux diagnostiques
    4676       call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
    4677       call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
    4678       call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    4679       call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
    4680       zsortie1d(:)=lmax(:)
    4681       call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
    4682       call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
    4683       zsortie1d(:)=lmix(:)
    4684       call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
    4685       zsortie1d(:)=lentr(:)
    4686       call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
    4687 
    4688 c      print*,'17 OK convect8'
    4689 
    4690          do k=1,klev/10
    4691             write(str2,'(i2.2)') k
    4692             str10='wa'//str2
    4693             do l=1,nlay
    4694                do ig=1,ngrid
    4695                   zsortie(ig,l)=wa(ig,k,l)
    4696                enddo
    4697             enddo   
    4698             CALL wrgradsfi(1,nlay,zsortie,str10,str10)
    4699             do l=1,nlay
    4700                do ig=1,ngrid
    4701                   zsortie(ig,l)=larg_part(ig,k,l)
    4702                enddo
    4703             enddo
    4704             str10='la'//str2
    4705             CALL wrgradsfi(1,nlay,zsortie,str10,str10)
    4706          enddo
    4707 
    4708 
    4709 c     print*,'18 OK convect8'
    4710 c      endif
    4711       print*,'Fin des wrgradsfi'
    4712 #endif
    4713 
    4714       endif
    4715 
    4716 c     if(wa_moy(1,4).gt.1.e-10) stop
    4717 
    4718       print*,'19 OK convect8'
     4291!     print*,'19 OK convect8'
    47194292      return
    47204293      end
    … …  
    52394812
    52404813      real count_time
    5241       integer isplit,nsplit,ialt
    5242       parameter (nsplit=10)
    5243       data isplit/0/
    5244       save isplit
    5245 c$OMP THREADPRIVATE(isplit)
     4814      integer ialt
    52464815
    52474816      logical sorties
    … …  
    60075576      enddo
    60085577
    6009 cdeja fait
    6010 c      do l=1,nlay
    6011 c         do ig=1,ngrid
    6012 c            detr(ig,l)=fm(ig,l)+entr(ig,l)-fm(ig,l+1)
    6013 c            if (detr(ig,l).lt.0.) then
    6014 c                entr(ig,l)=entr(ig,l)-detr(ig,l)
    6015 c                detr(ig,l)=0.
    6016 c     print*,'WARNING !!! detrainement negatif ',ig,l
    6017 c            endif
    6018 c         enddo
    6019 c      enddo
    6020 
    6021 c     print*,'15 OK convect8'
    6022 
    6023       isplit=isplit+1
    6024 
    6025 
    6026 c #define und
    6027         goto 123
    6028 #ifdef und
    6029       CALL writeg1d(1,nlay,wd,'wd      ','wd      ')
    6030       CALL writeg1d(1,nlay,zwa,'wa      ','wa      ')
    6031       CALL writeg1d(1,nlay,fracd,'fracd      ','fracd      ')
    6032       CALL writeg1d(1,nlay,fraca,'fraca      ','fraca      ')
    6033       CALL writeg1d(1,nlay,wa_moy,'wam         ','wam         ')
    6034       CALL writeg1d(1,nlay,zla,'la      ','la      ')
    6035       CALL writeg1d(1,nlay,zld,'ld      ','ld      ')
    6036       CALL writeg1d(1,nlay,pt,'pt      ','pt      ')
    6037       CALL writeg1d(1,nlay,zh,'zh      ','zh      ')
    6038       CALL writeg1d(1,nlay,zha,'zha      ','zha      ')
    6039       CALL writeg1d(1,nlay,zu,'zu      ','zu      ')
    6040       CALL writeg1d(1,nlay,zv,'zv      ','zv      ')
    6041       CALL writeg1d(1,nlay,zo,'zo      ','zo      ')
    6042       CALL writeg1d(1,nlay,wh,'wh      ','wh      ')
    6043       CALL writeg1d(1,nlay,wu,'wu      ','wu      ')
    6044       CALL writeg1d(1,nlay,wv,'wv      ','wv      ')
    6045       CALL writeg1d(1,nlay,wo,'w15uo     ','wXo     ')
    6046       CALL writeg1d(1,nlay,zdhadj,'zdhadj      ','zdhadj      ')
    6047       CALL writeg1d(1,nlay,pduadj,'pduadj      ','pduadj      ')
    6048       CALL writeg1d(1,nlay,pdvadj,'pdvadj      ','pdvadj      ')
    6049       CALL writeg1d(1,nlay,pdoadj,'pdoadj      ','pdoadj      ')
    6050       CALL writeg1d(1,nlay,entr  ,'entr        ','entr        ')
    6051       CALL writeg1d(1,nlay,detr  ,'detr        ','detr        ')
    6052       CALL writeg1d(1,nlay,fm    ,'fm          ','fm          ')
    6053 
    6054       CALL writeg1d(1,nlay,pdtadj,'pdtadj    ','pdtadj    ')
    6055       CALL writeg1d(1,nlay,pplay,'pplay     ','pplay     ')
    6056       CALL writeg1d(1,nlay,pplev,'pplev     ','pplev     ')
    6057 
    6058 c   recalcul des flux en diagnostique...
    6059 c     print*,'PAS DE TEMPS ',ptimestep
    6060        call dt2F(pplev,pplay,pt,pdtadj,wh)
    6061       CALL writeg1d(1,nlay,wh,'wh2     ','wh2     ')
    6062 #endif
    6063 123   continue
    6064 ! #define troisD
    6065 #ifdef troisD
    6066 c       if (sorties) then
    6067       print*,'Debut des wrgradsfi'
    6068 
    6069 c      print*,'16 OK convect8'
    6070          call wrgradsfi(1,nlay,wd,'wd        ','wd        ')
    6071          call wrgradsfi(1,nlay,zwa,'wa        ','wa        ')
    6072          call wrgradsfi(1,nlay,fracd,'fracd     ','fracd     ')
    6073          call wrgradsfi(1,nlay,fraca,'fraca     ','fraca     ')
    6074          call wrgradsfi(1,nlay,xxx,'xxx       ','xxx       ')
    6075          call wrgradsfi(1,nlay,wa_moy,'wam       ','wam       ')
    6076 c      print*,'WA6 ',wa_moy
    6077          call wrgradsfi(1,nlay,zla,'la        ','la        ')
    6078          call wrgradsfi(1,nlay,zld,'ld        ','ld        ')
    6079          call wrgradsfi(1,nlay,pt,'pt        ','pt        ')
    6080          call wrgradsfi(1,nlay,zh,'zh        ','zh        ')
    6081          call wrgradsfi(1,nlay,zha,'zha        ','zha        ')
    6082          call wrgradsfi(1,nlay,zua,'zua        ','zua        ')
    6083          call wrgradsfi(1,nlay,zva,'zva        ','zva        ')
    6084          call wrgradsfi(1,nlay,zu,'zu        ','zu        ')
    6085          call wrgradsfi(1,nlay,zv,'zv        ','zv        ')
    6086          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    6087          call wrgradsfi(1,nlay,wh,'wh        ','wh        ')
    6088          call wrgradsfi(1,nlay,wu,'wu        ','wu        ')
    6089          call wrgradsfi(1,nlay,wv,'wv        ','wv        ')
    6090          call wrgradsfi(1,nlay,wo,'wo        ','wo        ')
    6091          call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    6092          call wrgradsfi(1,nlay,zdhadj,'zdhadj    ','zdhadj    ')
    6093          call wrgradsfi(1,nlay,pduadj,'pduadj    ','pduadj    ')
    6094          call wrgradsfi(1,nlay,pdvadj,'pdvadj    ','pdvadj    ')
    6095          call wrgradsfi(1,nlay,pdoadj,'pdoadj    ','pdoadj    ')
    6096          call wrgradsfi(1,nlay,entr,'entr      ','entr      ')
    6097          call wrgradsfi(1,nlay,detr,'detr      ','detr      ')
    6098          call wrgradsfi(1,nlay,fm,'fm        ','fm        ')
    6099          call wrgradsfi(1,nlay,fmc,'fmc       ','fmc       ')
    6100          call wrgradsfi(1,nlay,zw2,'zw2       ','zw2       ')
    6101          call wrgradsfi(1,nlay,ztva,'ztva      ','ztva      ')
    6102          call wrgradsfi(1,nlay,ztv,'ztv       ','ztv       ')
    6103 
    6104          call wrgradsfi(1,nlay,zo,'zo        ','zo        ')
    6105          call wrgradsfi(1,nlay,larg_cons,'Lc        ','Lc        ')
    6106          call wrgradsfi(1,nlay,larg_detr,'Ldetr     ','Ldetr     ')
    6107 
    6108 cCR:nouveaux diagnostiques
    6109       call wrgradsfi(1,nlay,entr_star  ,'entr_star   ','entr_star   ')     
    6110       call wrgradsfi(1,nlay,f_star    ,'f_star   ','f_star   ')
    6111       call wrgradsfi(1,1,zmax,'zmax      ','zmax      ')
    6112       call wrgradsfi(1,1,zmix,'zmix      ','zmix      ')
    6113       zsortie1d(:)=lmax(:)
    6114       call wrgradsfi(1,1,zsortie1d,'lmax      ','lmax      ')
    6115       call wrgradsfi(1,1,wmax,'wmax      ','wmax      ')
    6116       zsortie1d(:)=lmix(:)
    6117       call wrgradsfi(1,1,zsortie1d,'lmix      ','lmix      ')
    6118       zsortie1d(:)=lentr(:)
    6119       call wrgradsfi(1,1,zsortie1d,'lentr      ','lentr     ')
    6120 
    6121 c      print*,'17 OK convect8'
     5578
    61225579
    61235580         do k=1,klev/10
    … …  
    61395596         enddo
    61405597
    6141 
    6142 c     print*,'18 OK convect8'
    6143 c      endif
    6144       print*,'Fin des wrgradsfi'
    6145 #endif
    6146 
    61475598      endif
    61485599
    6149 c     if(wa_moy(1,4).gt.1.e-10) stop
    6150 
    6151 c      print*,'19 OK convect8'
     5600
     5601
    61525602      return
    61535603      end
Note: See TracChangeset for help on using the changeset viewer.