Changeset 1294 for LMDZ4


Ignore:
Timestamp:
Jan 14, 2010, 3:35:30 PM (14 years ago)
Author:
Laurent Fairhead
Message:

Modifications pour la nouvelle version des thermiques (2009/2010) CR et FH

Location:
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/calltherm.F90

    r1146 r1294  
    4545!********************************************************
    4646!     declarations
     47      LOGICAL flag_bidouille_stratocu
    4748      real fmc_therm(klon,klev+1),zqasc(klon,klev)
    4849      real zqla(klon,klev)
    … …  
    187188     &      ,r_aspect_thermals,l_mix_thermals,w2di_thermals  &
    188189     &      ,tau_thermals)
    189           else if (iflag_thermals.ge.13) then
    190             CALL thermcell_main(itap,klon,klev,zdt  &
     190          else if (iflag_thermals==13.or.iflag_thermals==14) then
     191            CALL thermcellV0_main(itap,klon,klev,zdt  &
    191192     &      ,pplay,paprs,pphi,debut  &
    192193     &      ,u_seri,v_seri,t_seri,q_seri  &
    … …  
    197198     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
    198199     &      ,zmax0,f0,zw2,fraca)
     200          else if (iflag_thermals==15.or.iflag_thermals==16) then
     201            CALL thermcell_main(itap,klon,klev,zdt  &
     202     &      ,pplay,paprs,pphi,debut  &
     203     &      ,u_seri,v_seri,t_seri,q_seri  &
     204     &      ,d_u_the,d_v_the,d_t_the,d_q_the  &
     205     &      ,zfm_therm,zentr_therm,zdetr_therm,zqasc,zqla,lmax  &
     206     &      ,ratqscth,ratqsdiff,zqsatth  &
     207     &      ,r_aspect_thermals,l_mix_thermals  &
     208     &      ,tau_thermals,Ale,Alp,lalim_conv,wght_th &
     209     &      ,zmax0,f0,zw2,fraca)
     210           if (prt_level.gt.10) write(lunout,*)'Apres thermcell_main OK'
     211         else
     212            STOP'Cas des thermiques non prevu'
    199213         endif
    200214
     215       flag_bidouille_stratocu=iflag_thermals.lt.14.or.iflag_thermals.lt.16
    201216
    202217      fact(:)=0.
    203218      DO i=1,klon
    204        logexpr1(i)=iflag_thermals.lt.14.or.weak_inversion(i).gt.0.5
     219       logexpr1(i)=flag_bidouille_stratocu.or.weak_inversion(i).gt.0.5
    205220       IF(logexpr1(i)) fact(i)=1./float(nsplit_thermals)
    206221      ENDDO
    … …  
    235250            qmemoire(:,:)=q_seri(:,:)
    236251            q_seri(:,:) = q_seri(:,:) + d_q_the(:,:)
     252           if (prt_level.gt.10) write(lunout,*)'Apres apres thermcell_main OK'
    237253
    238254       DO i=1,klon

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/physiq.F

    r1287 r1294  
    25482548
    25492549!   les ratqs sont une combinaison de ratqss et ratqsc
    2550        print*,'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
     2550       write(lunout,*)'PHYLMD NOUVEAU TAU_RATQS ',tau_ratqs
    25512551
    25522552         if (tau_ratqs>1.e-10) then

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_closure.F90

    r1146 r1294  
     1!
     2! $Header$
     3!
    14      SUBROUTINE thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
    2      &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
     5     &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
    36
    47!-------------------------------------------------------------------------
    58!thermcell_closure: fermeture, determination de f
     9!
     10! Modification 7 septembre 2009
     11! 1. On enleve alim_star_tot des arguments pour le recalculer et etre ainis
     12! coherent avec l'integrale au numerateur.
     13! 2. On ne garde qu'une version des couples wmax,zmax et wmax_sec,zmax_sec
     14! l'idee etant que le choix se fasse a l'appel de thermcell_closure
     15! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
    616!-------------------------------------------------------------------------
    717      IMPLICIT NONE
    … …  
    919#include "iniprint.h"
    1020#include "thermcell.h"
    11       INTEGER ngrid,nlay
    12       INTEGER ig,k       
    13       REAL r_aspect,ptimestep
    14       integer lev_out                           ! niveau pour les print
     21INTEGER ngrid,nlay
     22INTEGER ig,k       
     23REAL r_aspect,ptimestep
     24integer lev_out                           ! niveau pour les print
    1525
    16       INTEGER lalim(ngrid)
    17       REAL alim_star(ngrid,nlay)
    18       REAL alim_star_tot(ngrid)
    19       REAL rho(ngrid,nlay)
    20       REAL zlev(ngrid,nlay)
    21       REAL zmax(ngrid),zmax_sec(ngrid)
    22       REAL wmax(ngrid),wmax_sec(ngrid)
    23       real zdenom
     26INTEGER lalim(ngrid)
     27REAL alim_star(ngrid,nlay)
     28REAL f_star(ngrid,nlay+1)
     29REAL rho(ngrid,nlay)
     30REAL zlev(ngrid,nlay)
     31REAL zmax(ngrid)
     32REAL wmax(ngrid)
     33REAL zdenom(ngrid)
     34REAL alim_star2(ngrid)
     35REAL f(ngrid)
    2436
    25       REAL alim_star2(ngrid)
     37REAL alim_star_tot(ngrid)
     38INTEGER llmax
    2639
    27       REAL f(ngrid)
     40!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     41print*,'THERMCELL CLOSURE 26E'
    2842
    29       do ig=1,ngrid
    30          alim_star2(ig)=0.
    31       enddo
    32       do ig=1,ngrid
    33          if (alim_star(ig,1).LT.1.e-10) then
    34             f(ig)=0.
    35          else   
    36              do k=1,lalim(ig)
    37                 alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
    38      &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
    39              enddo
    40              zdenom=max(500.,zmax(ig))*r_aspect*alim_star2(ig)
    41              if (zdenom<1.e-14) then
    42                 print*,'ig=',ig
    43                 print*,'alim_star2',alim_star2(ig)
    44                 print*,'zmax',zmax(ig)
    45                 print*,'r_aspect',r_aspect
    46                 print*,'zdenom',zdenom
    47                 print*,'alim_star',alim_star(ig,:)
    48                 print*,'zmax_sec',zmax_sec(ig)
    49                 print*,'wmax_sec',wmax_sec(ig)
    50                 stop
    51              endif
    52              if ((zmax_sec(ig).gt.1.e-10).and.(iflag_thermals_ed.eq.0)) then
    53              f(ig)=wmax_sec(ig)*alim_star_tot(ig)/(max(500.,zmax_sec(ig))*r_aspect  &
    54      &             *alim_star2(ig))
    55 !            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
    56 !    &                     zmax_sec(ig))*wmax_sec(ig))
    57              if(prt_level.GE.10) write(lunout,*)'closure dry',f(ig),wmax_sec(ig),alim_star_tot(ig),zmax_sec(ig)
    58              else
    59              f(ig)=wmax(ig)*alim_star_tot(ig)/zdenom
    60 !            f(ig)=f(ig)+(f0(ig)-f(ig))*exp((-ptimestep/  &
    61 !     &                     zmax(ig))*wmax(ig))
    62              if(prt_level.GE.10) print*,'closure moist',f(ig),wmax(ig),alim_star_tot(ig),zmax(ig)
    63              endif
    64          endif
    65 !         f0(ig)=f(ig)
    66       enddo
    67       if (prt_level.ge.1) print*,'apres fermeture'
     43alim_star2(:)=0.
     44alim_star_tot(:)=0.
     45f(:)=0.
    6846
    69 !
     47! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
     48llmax=1
     49do ig=1,ngrid
     50   if (lalim(ig)>llmax) llmax=lalim(ig)
     51enddo
     52
     53
     54! Calcul des integrales sur la verticale de alim_star et de
     55!   alim_star^2/(rho dz)
     56do k=1,llmax-1
     57   do ig=1,ngrid
     58      if (k<lalim(ig)) then
     59         alim_star2(ig)=alim_star2(ig)+alim_star(ig,k)**2  &
     60&                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
     61         alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,k)
     62      endif
     63   enddo
     64enddo
     65
     66
     67do ig=1,ngrid
     68   if (alim_star2(ig)>1.e-10) then
     69      f(ig)=wmax(ig)*alim_star_tot(ig)/  &
     70&     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
     71   endif
     72enddo
     73
     74
     75!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     76! TESTS POUR UNE NOUVELLE FERMETURE DANS LAQUELLE ALIM_STAR NE SERAIT
     77! PAS NORMALISE
     78!           f(ig)=f(ig)*f_star(ig,2)/(f_star(ig,lalim(ig)))
     79!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     80
    7081      return
    7182      end

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_dry.F90

    r938 r1294  
     1!
     2! $Header$
     3!
    14       SUBROUTINE thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
    25     &                            lalim,lmin,zmax,wmax,lev_out)
    … …  
    47!--------------------------------------------------------------------------
    58!thermcell_dry: calcul de zmax et wmax du thermique sec
     9! Calcul de la vitesse maximum et de la hauteur maximum pour un panache
     10! ascendant avec une fonction d'alimentation alim_star et sans changement
     11! de phase.
     12! Le calcul pourrait etre sans doute simplifier.
     13! La temperature potentielle virtuelle dans la panache ascendant est
     14! la temperature potentielle virtuelle pondérée par alim_star.
    615!--------------------------------------------------------------------------
     16
    717       IMPLICIT NONE
    818#include "YOMCST.h"       
    … …  
    4858!calcul de la vitesse a partir de la CAPE en melangeant thetav
    4959
    50 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    51 ! A eliminer
    52 ! Ce if complique etait fait pour reperer la premiere couche instable
    53 ! Ici, c'est lmin.
    54 !
    55 !       do l=1,nlay-2
    56 !         do ig=1,ngrid
    57 !            if (ztv(ig,l).gt.ztv(ig,l+1)  &
    58 !     &         .and.alim_star(ig,l).gt.1.e-10  &
    59 !     &         .and.zw2(ig,l).lt.1e-10) then
    60 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    61 
    6260
    6361! Calcul des F^*, integrale verticale de E^*
    … …  
    8482!  Premiere couche du panache thermique
    8583!------------------------------------------------------------------------
     84
    8685               zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
    8786     &                     *(zlev(ig,l+1)-zlev(ig,l))  &
    … …  
    9695! 3. la vitesse au carré en haut zw2(ig,l+1)
    9796!------------------------------------------------------------------------
    98 
    99 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    100 !  A eliminer : dans cette version, si zw2 est > 0 on a un therique.
    101 !  et donc, au dessus, f_star(ig,l+1) est forcement suffisamment
    102 !  grand puisque on n'a pas de detrainement.
    103 !  f_star est une fonction croissante.
    104 !  c'est donc vraiment sur zw2 uniquement qu'il faut faire le test.
    105 !           else if ((zw2(ig,l).ge.1e-10).and.  &
    106 !    &               (f_star(ig,l)+alim_star(ig,l).gt.1.e-10)) then
    107 !              f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)
    108 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    10997
    11098            else if (zw2(ig,l).ge.1e-10) then
    … …  
    145133       if (prt_level.ge.1) print*,'fin calcul zw2'
    146134!
    147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    148 ! A eliminer :
    149 ! Ce calcul de lmax est fait en meme temps que celui de linter, plus haut
    150 ! Calcul de la couche correspondant a la hauteur du thermique
    151 !      do ig=1,ngrid
    152 !         lmax(ig)=lalim(ig)
    153 !      enddo
    154 !      do ig=1,ngrid
    155 !         do l=nlay,lalim(ig)+1,-1
    156 !            if (zw2(ig,l).le.1.e-10) then
    157 !               lmax(ig)=l-1
    158 !            endif
    159 !         enddo
    160 !      enddo
    161 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    162 
    163 !   
    164135! Determination de zw2 max
    165136      do ig=1,ngrid
    … …  
    185156      do  ig=1,ngrid
    186157! calcul de zlevinter
    187 
    188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    189 ! FH A eliminer
    190 ! Simplification
    191 !          zlevinter(ig)=(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*  &
    192 !     &    linter(ig)+zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)  &
    193 !     &    -zlev(ig,lmax(ig)))
    194 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    195 
    196158          zlevinter(ig)=zlev(ig,lmax(ig)) + &
    197159     &    (linter(ig)-lmax(ig))*(zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))
    … …  
    199161      enddo
    200162
    201 ! Verification que lalim<=lmax
    202       do ig=1,ngrid
    203          if(lalim(ig)>lmax(ig)) then
    204            if ( prt_level > 1 ) THEN
    205             print*,'WARNING thermcell_dry ig=',ig,'  lalim=',lalim(ig),'  lmax(ig)=',lmax(ig)
    206            endif
    207            lmax(ig)=lalim(ig)
    208          endif
    209       enddo
    210      
    211163      RETURN
    212164      END

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_init.F90

    r1057 r1294  
     1!
     2! $Header$
     3!
    14      SUBROUTINE thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
    25     &                  lalim,lmin,alim_star,alim_star_tot,lev_out)
    … …  
    2629!def des alim_star tels que alim=f*alim_star     
    2730
    28       do l=1,nlay
    29          do ig=1,ngrid
    30             alim_star(ig,l)=0.
    31          enddo
    32       enddo
    33 ! determination de la longueur de la couche d entrainement
    34       do ig=1,ngrid
    35          lalim(ig)=1
    36       enddo
    3731
    38       if (iflag_thermals_ed.ge.1) then
    39 !si la première couche est instable, on declenche un thermique
     32      write(lunout,*)'THERM INIT V20C '
     33
     34      alim_star_tot(:)=0.
     35      alim_star(:,:)=0.
     36      lmin(:)=1
     37      lalim(:)=1
     38
     39      do l=1,nlay-1
    4040         do ig=1,ngrid
    41             if (ztv(ig,1).gt.ztv(ig,2)) then
    42                lmin(ig)=1
    43                lalim(ig)=2
    44                alim_star(ig,1)=1.
    45                alim_star_tot(ig)=alim_star(ig,1)
    46                if(prt_level.GE.10) print*,'init',alim_star(ig,1),alim_star_tot(ig)
    47             else
    48                 lmin(ig)=1
    49                 lalim(ig)=1
    50                 alim_star(ig,1)=0.
    51                 alim_star_tot(ig)=0.
    52             endif
    53          enddo
    54  
    55          else
    56 !else iflag_thermals_ed=0 ancienne def de l alim
    57 
    58 !on ne considere que les premieres couches instables
    59       do l=nlay-2,1,-1
    60          do ig=1,ngrid
    61             if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
    62      &          ztv(ig,l+1).le.ztv(ig,l+2)) then
    63                lalim(ig)=l+1
    64             endif
    65           enddo
    66       enddo
    67 
    68 ! determination du lmin: couche d ou provient le thermique
    69 
    70       do ig=1,ngrid
    71 ! FH initialisation de lmin a nlay plutot que 1.
    72 !        lmin(ig)=nlay
    73          lmin(ig)=1
    74       enddo
    75       do l=nlay,2,-1
    76          do ig=1,ngrid
    77             if (ztv(ig,l-1).gt.ztv(ig,l)) then
    78                lmin(ig)=l-1
     41            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
     42               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     43     &                       *sqrt(zlev(ig,l+1))
     44               lalim(:)=l+1
     45               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    7946            endif
    8047         enddo
    8148      enddo
    82 !
    83       zzalim(:)=0.
    84       do l=1,nlay-1
     49      do l=1,nlay
    8550         do ig=1,ngrid
    86              if (l<lalim(ig)) then
    87                 zzalim(ig)=zzalim(ig)+zlay(ig,l)*(ztv(ig,l)-ztv(ig,l+1))
    88              endif
    89           enddo
    90       enddo
    91       do ig=1,ngrid
    92           if (lalim(ig)>1) then
    93              zzalim(ig)=zlay(ig,1)+zzalim(ig)/(ztv(ig,1)-ztv(ig,lalim(ig)))
    94           else
    95              zzalim(ig)=zlay(ig,1)
    96           endif
    97       enddo
    98 
    99       if(prt_level.GE.10) print*,'ZZALIM LALIM ',zzalim,lalim,zlay(1,lalim(1))
    100 
    101 ! definition de l'entrainement des couches
    102       if (1.eq.1) then
    103       do l=1,nlay-1
    104          do ig=1,ngrid
    105             if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
    106      &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
    107 !def possibles pour alim_star: zdthetadz, dthetadz, zdtheta
    108              alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    109      &                       *sqrt(zlev(ig,l+1))
    110             endif
    111          enddo
    112       enddo
    113       else
    114       do l=1,nlay-1
    115          do ig=1,ngrid
    116             if (ztv(ig,l).gt.ztv(ig,l+1).and.  &
    117      &          l.ge.lmin(ig).and.l.lt.lalim(ig)) then
    118              alim_star(ig,l)=max(3.*zzalim(ig)-zlay(ig,l),0.) &
    119      &        *(zlev(ig,l+1)-zlev(ig,l))
    120             endif
    121          enddo
    122       enddo
    123       endif
    124      
    125 ! pas de thermique si couche 1 stable
    126       do ig=1,ngrid
    127 !CRnouveau test
    128         if (alim_star(ig,1).lt.1.e-10) then
    129             do l=1,nlay
    130                 alim_star(ig,l)=0.
    131             enddo
    132             lmin(ig)=1
    133          endif
    134       enddo
    135 ! calcul de l alimentation totale
    136       do ig=1,ngrid
    137          alim_star_tot(ig)=0.
    138       enddo
    139       do l=1,nlay
    140          do ig=1,ngrid
    141             alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    142          enddo
    143       enddo
    144 !
    145 ! Calcul entrainement normalise
    146       do l=1,nlay
    147          do ig=1,ngrid
    148             if (alim_star_tot(ig).gt.1.e-10) then
     51            if (alim_star_tot(ig) > 1.e-10 ) then
    14952               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    15053            endif
    15154         enddo
    15255      enddo
    153        
    154 !on remet alim_star_tot a 1
    155       do ig=1,ngrid
    156          alim_star_tot(ig)=1.
    157       enddo
     56      alim_star_tot(:)=1.
    15857
    159       endif
    160 !endif iflag_thermals_ed
    161       return
     58      return
    16259      end 

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_main.F90

    r1146 r1294  
    2222!   de "thermiques" explicitement representes avec processus nuageux
    2323!
    24 !   Réécriture à partir d'un listing papier à Habas, le 14/02/00
    25 !
    26 !   le thermique est supposé homogène et dissipé par mélange avec
    27 !   son environnement. la longueur l_mix contrôle l'efficacité du
    28 !   mélange
    29 !
    30 !   Le calcul du transport des différentes espèces se fait en prenant
     24!   Reecriture a partir d'un listing papier a Habas, le 14/02/00
     25!
     26!   le thermique est suppose homogene et dissipe par melange avec
     27!   son environnement. la longueur l_mix controle l'efficacite du
     28!   melange
     29!
     30!   Le calcul du transport des differentes especes se fait en prenant
    3131!   en compte:
    3232!     1. un flux de masse montant
    … …  
    8585      real linter(klon)
    8686      real zmix(klon)
    87       real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1)
     87      real zmax(klon),zw2(klon,klev+1),ztva(klon,klev),zw_est(klon,klev+1),ztva_est(klon,klev)
    8888!      real fraca(klon,klev)
    8989
    … …  
    115115! FH probleme de dimensionnement avec l'allocation dynamique
    116116!     common/comtherm/thetath2,wth2
     117      real wq(klon,klev)
     118      real wthl(klon,klev)
     119      real wthv(klon,klev)
    117120   
    118121      real ratqscth(klon,klev)
    … …  
    142145      real f_star(klon,klev+1),entr_star(klon,klev)
    143146      real detr_star(klon,klev)
    144       real alim_star_tot(klon),alim_star2(klon)
     147      real alim_star_tot(klon)
    145148      real alim_star(klon,klev)
     149      real alim_star_clos(klon,klev)
    146150      real f(klon), f0(klon)
    147151!FH/IM     save f0
    … …  
    149153      logical debut
    150154       real seuil
     155      real csc(klon,klev)
    151156
    152157!
    … …  
    220225      ENDIF
    221226!
    222 !Initialisation
    223 !
    224 !    IF (1.eq.0) THEN
    225 !     do ig=1,klon     
    226 !FH/IM 130308     if ((debut).or.((.not.debut).and.(f0(ig).lt.1.e-10))) then
    227 !     if ((.not.debut).and.(f0(ig).lt.1.e-10)) then
    228 !           f0(ig)=1.e-5
    229 !           zmax0(ig)=40.
    230 !v1d        therm=.false.
    231 !     endif
    232 !     enddo
    233 !    ENDIF !(1.eq.0) THEN
    234      if (prt_level.ge.10)write(lunout,*)                                &
    235     &     'WARNING thermcell_main f0=max(f0,1.e-2)'
     227     write(lunout,*)'WARNING thermcell_main f0=max(f0,1.e-2)'
    236228     do ig=1,klon
    237229      if (prt_level.ge.20) then
    … …  
    239231      endif
    240232         f0(ig)=max(f0(ig),1.e-2)
     233         zmax0(ig)=max(zmax0(ig),40.)
    241234!IMmarche pas ?!       if (f0(ig)<1.e-2) f0(ig)=1.e-2
    242235     enddo
    … …  
    364357
    365358!------------------------------------------------------------------
    366 !  1. alim_star est le profil vertical de l'alimentation à la base du
    367 !     panache thermique, calculé à partir de la flotabilité de l'air sec
     359!  1. alim_star est le profil vertical de l'alimentation a la base du
     360!     panache thermique, calcule a partir de la flotabilite de l'air sec
    368361!  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
    369362!------------------------------------------------------------------
    370363!
    371364      entr_star=0. ; detr_star=0. ; alim_star=0. ; alim_star_tot=0.
    372       CALL thermcell_init(ngrid,nlay,ztv,zlay,zlev,  &
    373      &                    lalim,lmin,alim_star,alim_star_tot,lev_out)
    374 
    375 call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lmin  ')
    376 call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_init lalim ')
    377 
    378 
    379       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_init'
    380       if (prt_level.ge.10) then
    381          write(lunout1,*) 'Dans thermcell_main 1'
    382          write(lunout1,*) 'lmin ',lmin(igout)
    383          write(lunout1,*) 'lalim ',lalim(igout)
    384          write(lunout1,*) ' ig l alim_star thetav'
    385          write(lunout1,'(i6,i4,2e15.5)') (igout,l,alim_star(igout,l) &
    386      &   ,ztv(igout,l),l=1,lalim(igout)+4)
    387       endif
    388 
    389 !v1d      do ig=1,klon
    390 !v1d     if (alim_star(ig,1).gt.1.e-10) then
    391 !v1d     therm=.true.
    392 !v1d     endif
    393 !v1d      enddo
     365      lmin=1
     366
    394367!-----------------------------------------------------------------------------
    395368!  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
    396369!     panache sec conservatif (e=d=0) alimente selon alim_star
    397370!     Il s'agit d'un calcul de type CAPE
    398 !     zmax_sec est utilisé pour déterminer la géométrie du thermique.
     371!     zmax_sec est utilise pour determiner la geometrie du thermique.
    399372!------------------------------------------------------------------------------
    400 !
     373!---------------------------------------------------------------------------------
     374!calcul du melange et des variables dans le thermique
     375!--------------------------------------------------------------------------------
     376!
     377      if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
     378!IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
     379      CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
     380     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     381     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     382     &           ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
     383     &            ,lev_out,lunout1,igout)
     384      if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
     385
     386      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
     387      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
     388
     389      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
     390      if (prt_level.ge.10) then
     391         write(lunout1,*) 'Dans thermcell_main 2'
     392         write(lunout1,*) 'lmin ',lmin(igout)
     393         write(lunout1,*) 'lalim ',lalim(igout)
     394         write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
     395         write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
     396     &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
     397      endif
     398
     399!-------------------------------------------------------------------------------
     400! Calcul des caracteristiques du thermique:zmax,zmix,wmax
     401!-------------------------------------------------------------------------------
     402!
     403      CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
     404     &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
     405
     406
     407
     408      call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
     409      call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
     410      call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
     411      call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
     412
     413      if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
     414
     415!-------------------------------------------------------------------------------
     416! Fermeture,determination de f
     417!-------------------------------------------------------------------------------
     418!
     419!
     420      write(lunout,*)'THERM NOUVEAU RIO2009 31B'
    401421      CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,  &
    402      &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
     422    &                      lalim,lmin,zmax_sec,wmax_sec,lev_out)
    403423
    404424call test_ltherm(ngrid,nlay,pplev,pplay,lmin,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
    … …  
    417437
    418438
    419 !---------------------------------------------------------------------------------
    420 !calcul du melange et des variables dans le thermique
    421 !--------------------------------------------------------------------------------
    422 !
    423       if (prt_level.ge.1) print*,'avant thermcell_plume ',lev_out
    424 !IM 140508   CALL thermcell_plume(ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    425       CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    426      &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
    427      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    428      &           ztla,zqla,zqta,zha,zw2,zw_est,zqsatth,lmix,lmix_bis,linter &
    429      &            ,lev_out,lunout1,igout)
    430       if (prt_level.ge.1) print*,'apres thermcell_plume ',lev_out
    431 
    432       call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
    433       call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
    434 
    435       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_plume'
    436       if (prt_level.ge.10) then
    437          write(lunout1,*) 'Dans thermcell_main 2'
    438          write(lunout1,*) 'lmin ',lmin(igout)
    439          write(lunout1,*) 'lalim ',lalim(igout)
    440          write(lunout1,*) ' ig l alim_star entr_star detr_star f_star '
    441          write(lunout1,'(i6,i4,4e15.5)') (igout,l,alim_star(igout,l),entr_star(igout,l),detr_star(igout,l) &
    442      &    ,f_star(igout,l+1),l=1,nint(linter(igout))+5)
    443       endif
    444 
    445 !-------------------------------------------------------------------------------
    446 ! Calcul des caracteristiques du thermique:zmax,zmix,wmax
    447 !-------------------------------------------------------------------------------
    448 !
    449       CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,  &
    450      &           zlev,lmax,zmax,zmax0,zmix,wmax,lev_out)
    451 
    452 
    453       call test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
    454       call test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
    455       call test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
    456       call test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
    457 
    458       if (prt_level.ge.1) print*,'thermcell_main apres thermcell_height'
    459 
    460 !-------------------------------------------------------------------------------
    461 ! Fermeture,determination de f
    462 !-------------------------------------------------------------------------------
    463 !
    464 !avant closure: on redéfinit lalim, alim_star_tot et alim_star
    465 !       do ig=1,klon
    466 !       do l=2,lalim(ig)
    467 !       alim_star(ig,l)=entr_star(ig,l)
    468 !       entr_star(ig,l)=0.
    469 !       enddo
    470 !       enddo
    471 
     439print*,'THERM 26JJJ'
     440
     441! Choix de la fonction d'alimentation utilisee pour la fermeture.
     442! Apparemment sans importance
     443      alim_star_clos(:,:)=alim_star(:,:)
     444      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
     445
     446! Appel avec la version seche
    472447      CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
    473      &   zlev,lalim,alim_star,alim_star_tot,zmax_sec,wmax_sec,zmax,wmax,f,lev_out)
     448     &   zlev,lalim,alim_star_clos,f_star,zmax_sec,wmax_sec,f,lev_out)
     449
     450!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     451! Appel avec les zmax et wmax tenant compte de la condensation
     452! Semble moins bien marcher
     453!     CALL thermcell_closure(ngrid,nlay,r_aspect,ptimestep,rho,  &
     454!    &   zlev,lalim,alim_star,f_star,zmax,wmax,f,lev_out)
     455!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    474456
    475457      if(prt_level.ge.1)print*,'thermcell_closure apres thermcell_closure'
    … …  
    484466! Test valable seulement en 1D mais pas genant
    485467      if (.not. (f0(1).ge.0.) ) then
    486            stop 'Dans thermcell_main'
     468           stop'Dans thermcell_main'
    487469      endif
    488470
    … …  
    511493         fm0=(1.-lambda)*fm+lambda*fm0
    512494         entr0=(1.-lambda)*entr+lambda*entr0
    513 !        detr0=(1.-lambda)*detr+lambda*detr0
     495         detr0=(1.-lambda)*detr+lambda*detr0
    514496      else
    515497         fm0=fm
    … …  
    560542     &    ,fraca,zmax  &
    561543     &    ,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
    562 !IM 050508    &    ,zu,zv,pduadj,pdvadj,zua,zva,igout,lev_out)
     544
    563545      else
    564546
    … …  
    636618!
    637619      if (prt_level.ge.10) print*,'14f OK convect8 ig,l,zha zh zpspsk ',ig,l,zha(ig,l),zh(ig,l),zpspsk(ig,l)
    638             thetath2(ig,l)=zf2*(zha(ig,l)-zh(ig,l)/zpspsk(ig,l))**2
     620            thetath2(ig,l)=zf2*(ztla(ig,l)-zthl(ig,l))**2
    639621            if(zw2(ig,l).gt.1.e-10) then
    640622             wth2(ig,l)=zf2*(zw2(ig,l))**2
    … …  
    651633         enddo
    652634      enddo
     635!calcul des flux: q, thetal et thetav
     636      do l=1,nlay
     637         do ig=1,ngrid
     638      wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
     639      wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
     640      wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
     641         enddo
     642      enddo
    653643!calcul de ale_bl et alp_bl
    654 !pour le calcul d'une valeur intégrée entre la surface et lmax
     644!pour le calcul d'une valeur integree entre la surface et lmax
    655645      do ig=1,ngrid
    656646      alp_int(ig)=0.
    … …  
    782772!     print*,'15 OK convect8'
    783773
     774#ifdef wrgrads_thermcell
    784775      if (prt_level.ge.1) print*,'thermcell_main sorties 3D'
    785 #ifdef wrgrads_thermcell
    786776#include "thermcell_out3d.h"
    787777#endif

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_out3d.h

    r1029 r1294  
    2727         call wrgradsfi(1,nlay,q2(igout,1:klev),'q2       ','q2       ')
    2828!
     29!
     30         call wrgradsfi(1,nlay,wthl(igout,1:klev),'wthl       ','wthl       ')
     31         call wrgradsfi(1,nlay,wthv(igout,1:klev),'wthv       ','wthv       ')
     32         call wrgradsfi(1,nlay,wq(igout,1:klev),'wq       ','wq       ')
     33         
    2934         call wrgradsfi(1,nlay,ztva(igout,1:klev),'ztva      ','ztva      ')
    3035         call wrgradsfi(1,nlay,ztv(igout,1:klev),'ztv       ','ztv       ')
    … …  
    5358      call wrgradsfi(1,1,f(igout),'f      ','f      ')
    5459      call wrgradsfi(1,1,alim_star_tot(igout),'a_s_t      ','a_s_t      ')
    55       call wrgradsfi(1,1,alim_star2(igout),'a_2      ','a_2      ')
    5660      call wrgradsfi(1,1,zmax(igout),'zmax      ','zmax      ')
    5761      call wrgradsfi(1,1,zmax_sec(igout),'z_sec      ','z_sec      ')

  • LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/thermcell_plume.F90

    r1057 r1294  
     1!
     2! $Header$
     3!
    14      SUBROUTINE thermcell_plume(itap,ngrid,klev,ptimestep,ztv,zthl,po,zl,rhobarz,  &
    2      &           zlev,pplev,pphi,zpspsk,l_mix,r_aspect,alim_star,alim_star_tot,  &
    3      &           lalim,zmax_sec,f0,detr_star,entr_star,f_star,ztva,  &
    4      &           ztla,zqla,zqta,zha,zw2,w_est,zqsatth,lmix,lmix_bis,linter &
     5     &           zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,  &
     6     &           lalim,f0,detr_star,entr_star,f_star,csc,ztva,  &
     7     &           ztla,zqla,zqta,zha,zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,linter &
    58     &           ,lev_out,lunout1,igout)
    69
    … …  
    3134      REAL zpspsk(ngrid,klev)
    3235      REAL alim_star(ngrid,klev)
    33       REAL zmax_sec(ngrid)
    3436      REAL f0(ngrid)
    35       REAL l_mix
    36       REAL r_aspect
    3737      INTEGER lalim(ngrid)
    3838      integer lev_out                           ! niveau pour les print
    … …  
    4444      REAL ztla(ngrid,klev)
    4545      REAL zqla(ngrid,klev)
    46       REAL zqla0(ngrid,klev)
    4746      REAL zqta(ngrid,klev)
    4847      REAL zha(ngrid,klev)
    … …  
    5049      REAL detr_star(ngrid,klev)
    5150      REAL coefc
    52       REAL detr_stara(ngrid,klev)
    53       REAL detr_starb(ngrid,klev)
    54       REAL detr_starc(ngrid,klev)
    55       REAL detr_star0(ngrid,klev)
    56       REAL detr_star1(ngrid,klev)
    57       REAL detr_star2(ngrid,klev)
    58 
    5951      REAL entr_star(ngrid,klev)
    60       REAL entr_star1(ngrid,klev)
    61       REAL entr_star2(ngrid,klev)
    6252      REAL detr(ngrid,klev)
    6353      REAL entr(ngrid,klev)
     54
     55      REAL csc(ngrid,klev)
    6456
    6557      REAL zw2(ngrid,klev+1)
    … …  
    7264      REAL zqsatth(ngrid,klev)
    7365      REAL zta_est(ngrid,klev)
     66      REAL zdw2
     67      REAL zw2modif
     68      REAL zeps
    7469
    7570      REAL linter(ngrid)
    … …  
    8075      INTEGER ig,l,k
    8176
     77      real zdz,zfact,zbuoy,zalpha
    8278      real zcor,zdelta,zcvm5,qlbef
    8379      real Tbef,qsatbef
    … …  
    8682      PARAMETER (DDT0=.01)
    8783      logical Zsat
    88       REAL fact_gamma,fact_epsilon
     84      LOGICAL active(ngrid),activetmp(ngrid)
     85      REAL fact_gamma,fact_epsilon,fact_gamma2
    8986      REAL c2(ngrid,klev)
    9087
     88      REAL zw2fact,expa
    9189      Zsat=.false.
    9290! Initialisation
    … …  
    9795         fact_epsilon=1.
    9896      else if (iflag_thermals_ed==1)  then
     97! Valeurs au moment de la premiere soumission des papiers
    9998         fact_gamma=1.
    100          fact_epsilon=1.
     99         fact_epsilon=0.002
     100         fact_gamma2=0.6
     101! Suggestions des Fleurs, Septembre 2009
     102         fact_epsilon=0.015
     103!test cr
     104!         fact_epsilon=0.002
     105         fact_gamma=0.9
     106         fact_gamma2=0.7
     107
    101108      else if (iflag_thermals_ed==2)  then
    102109         fact_gamma=1.
    103110         fact_epsilon=2.
     111      else if (iflag_thermals_ed==3)  then
     112         fact_gamma=3./4.
     113         fact_epsilon=3.
    104114      endif
    105115
    106       do l=1,klev
     116      write(lunout,*)'THERM 31H '
     117
     118! Initialisations des variables reeles
     119if (1==0) then
     120      ztva(:,:)=ztv(:,:)
     121      ztva_est(:,:)=ztva(:,:)
     122      ztla(:,:)=zthl(:,:)
     123      zqta(:,:)=po(:,:)
     124      zha(:,:) = ztva(:,:)
     125else
     126      ztva(:,:)=0.
     127      ztva_est(:,:)=0.
     128      ztla(:,:)=0.
     129      zqta(:,:)=0.
     130      zha(:,:) =0.
     131endif
     132
     133      zqla_est(:,:)=0.
     134      zqsatth(:,:)=0.
     135      zqla(:,:)=0.
     136      detr_star(:,:)=0.
     137      entr_star(:,:)=0.
     138      alim_star(:,:)=0.
     139      alim_star_tot(:)=0.
     140      csc(:,:)=0.
     141      detr(:,:)=0.
     142      entr(:,:)=0.
     143      zw2(:,:)=0.
     144      w_est(:,:)=0.
     145      f_star(:,:)=0.
     146      wa_moy(:,:)=0.
     147      linter(:)=1.
     148      linter(:)=1.
     149
     150! Initialisation des variables entieres
     151      lmix(:)=1
     152      lmix_bis(:)=2
     153      wmaxa(:)=0.
     154      lalim(:)=1
     155
     156!-------------------------------------------------------------------------
     157! On ne considere comme actif que les colonnes dont les deux premieres
     158! couches sont instables.
     159!-------------------------------------------------------------------------
     160      active(:)=ztv(:,1)>ztv(:,2)
     161
     162!-------------------------------------------------------------------------
     163! Definition de l'alimentation a l'origine dans thermcell_init
     164!-------------------------------------------------------------------------
     165      do l=1,klev-1
    107166         do ig=1,ngrid
    108             zqla_est(ig,l)=0.
    109             ztva_est(ig,l)=ztva(ig,l)
    110             zqsatth(ig,l)=0.
     167            if (ztv(ig,l)> ztv(ig,l+1) .and. ztv(ig,1)>=ztv(ig,l) ) then
     168               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
     169     &                       *sqrt(zlev(ig,l+1))
     170               lalim(:)=l+1
     171               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     172            endif
    111173         enddo
    112174      enddo
    113 
    114 !CR: attention test couche alim
    115 !     do l=2,klev
    116 !     do ig=1,ngrid
    117 !        alim_star(ig,l)=0.
    118 !     enddo
    119 !     enddo
    120 !AM:initialisations du thermique
    121       do k=1,klev
    122          do ig=1,ngrid
    123             ztva(ig,k)=ztv(ig,k)
    124             ztla(ig,k)=zthl(ig,k)
    125             zqla(ig,k)=0.
    126             zqta(ig,k)=po(ig,k)
    127 !
    128             ztva(ig,k) = ztla(ig,k)*zpspsk(ig,k)+RLvCp*zqla(ig,k)
    129             ztva(ig,k) = ztva(ig,k)/zpspsk(ig,k)
    130             zha(ig,k) = ztva(ig,k)
    131 !
    132          enddo
    133       enddo
    134       do k=1,klev
    135         do ig=1,ngrid
    136            detr_star(ig,k)=0.
    137            entr_star(ig,k)=0.
    138 
    139            detr_stara(ig,k)=0.
    140            detr_starb(ig,k)=0.
    141            detr_starc(ig,k)=0.
    142            detr_star0(ig,k)=0.
    143            zqla0(ig,k)=0.
    144            detr_star1(ig,k)=0.
    145            detr_star2(ig,k)=0.
    146            entr_star1(ig,k)=0.
    147            entr_star2(ig,k)=0.
    148 
    149            detr(ig,k)=0.
    150            entr(ig,k)=0.
    151         enddo
    152       enddo
    153       if (prt_level.ge.1) print*,'7 OK convect8'
    154       do k=1,klev+1
    155          do ig=1,ngrid
    156             zw2(ig,k)=0.
    157             w_est(ig,k)=0.
    158             f_star(ig,k)=0.
    159             wa_moy(ig,k)=0.
     175      do l=1,klev
     176         do ig=1,ngrid
     177            if (alim_star_tot(ig) > 1.e-10 ) then
     178               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
     179            endif
    160180         enddo
    161181      enddo
    162 
    163       if (prt_level.ge.1) print*,'8 OK convect8'
    164       do ig=1,ngrid
    165          linter(ig)=1.
    166          lmix(ig)=1
    167          lmix_bis(ig)=2
    168          wmaxa(ig)=0.
    169       enddo
    170 
    171 !-----------------------------------------------------------------------------------
    172 !boucle de calcul de la vitesse verticale dans le thermique
    173 !-----------------------------------------------------------------------------------
    174       do l=1,klev-1
    175          do ig=1,ngrid
    176 
    177 
    178 
    179 ! Calcul dans la premiere couche active du thermique (ce qu'on teste
    180 ! en disant que la couche est instable et que w2 en bas de la couche
    181 ! est nulle.
    182 
    183             if (ztv(ig,l).gt.ztv(ig,l+1)  &
    184      &         .and.alim_star(ig,l).gt.1.e-10  &
    185      &         .and.zw2(ig,l).lt.1e-10) then
    186 
    187 
     182      alim_star_tot(:)=1.
     183
     184
     185!------------------------------------------------------------------------------
     186! Calcul dans la premiere couche
     187! On decide dans cette version que le thermique n'est actif que si la premiere
     188! couche est instable.
     189! Pourrait etre change si on veut que le thermiques puisse se déclencher
     190! dans une couche l>1
     191!------------------------------------------------------------------------------
     192do ig=1,ngrid
    188193! Le panache va prendre au debut les caracteristiques de l'air contenu
    189194! dans cette couche.
    190                ztla(ig,l)=zthl(ig,l)
    191                zqta(ig,l)=po(ig,l)
    192                zqla(ig,l)=zl(ig,l)
    193                f_star(ig,l+1)=alim_star(ig,l)
    194 
    195                zw2(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
    196      &                     *(zlev(ig,l+1)-zlev(ig,l))  &
    197      &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
    198                w_est(ig,l+1)=zw2(ig,l+1)
     195    if (active(ig)) then
     196    ztla(ig,1)=zthl(ig,1)
     197    zqta(ig,1)=po(ig,1)
     198    zqla(ig,1)=zl(ig,1)
     199!cr: attention, prise en compte de f*(1)=1
     200    f_star(ig,2)=alim_star(ig,1)
     201    zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     202&                     *(zlev(ig,2)-zlev(ig,1))  &
     203&                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
     204    w_est(ig,2)=zw2(ig,2)
     205    endif
     206enddo
    199207!
    200208
    201 
    202             else if ((zw2(ig,l).ge.1e-10).and.  &
    203      &         (f_star(ig,l)+alim_star(ig,l)).gt.1.e-10) then
    204 !estimation du detrainement a partir de la geometrie du pas precedent
    205 !tests sur la definition du detr
    206 !calcul de detr_star et entr_star
    207 
    208 
    209 
    210 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    211 ! FH le test miraculeux de Catherine ? Le bout du tunel ?
    212 !               w_est(ig,3)=zw2(ig,2)*  &
    213 !    &                   ((f_star(ig,2))**2)  &
    214 !    &                   /(f_star(ig,2)+alim_star(ig,2))**2+  &
    215 !    &                   2.*RG*(ztva(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    216 !    &                   *(zlev(ig,3)-zlev(ig,2))
    217 !     if (l.gt.2) then
    218 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     209!==============================================================================
     210!boucle de calcul de la vitesse verticale dans le thermique
     211!==============================================================================
     212do l=2,klev-1
     213!==============================================================================
     214
     215
     216! On decide si le thermique est encore actif ou non
     217! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     218    do ig=1,ngrid
     219       active(ig)=active(ig) &
     220&                 .and. zw2(ig,l)>1.e-10 &
     221&                 .and. f_star(ig,l)+alim_star(ig,l)>1.e-10
     222    enddo
    219223
    220224
    … …  
    222226! Premier calcul de la vitesse verticale a partir de la temperature
    223227! potentielle virtuelle
    224 
    225 ! FH CESTQUOI CA ????
    226 #define int1d2
    227 !#undef int1d2
    228 #ifdef int1d2
    229       if (l.ge.2) then
    230 #else
    231       if (l.gt.2) then
    232 #endif
    233 
    234       if (1.eq.1) then
    235           w_est(ig,3)=zw2(ig,2)* &
    236      &      ((f_star(ig,2))**2) &
    237      &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
    238      &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
    239 !     &      *1./3. &
    240      &      *(zlev(ig,3)-zlev(ig,2))
    241        endif
     228!     if (1.eq.1) then
     229!         w_est(ig,3)=zw2(ig,2)* &
     230!    &      ((f_star(ig,2))**2) &
     231!    &      /(f_star(ig,2)+alim_star(ig,2))**2+ &
     232!    &      2.*RG*(ztva(ig,2)-ztv(ig,2))/ztv(ig,2) &
     233!    &      *(zlev(ig,3)-zlev(ig,2))
     234!      endif
    242235
    243236
    244237!---------------------------------------------------------------------------
    245 !calcul de l entrainement et du detrainement lateral
     238! calcul des proprietes thermodynamiques et de la vitesse de la couche l
     239! sans tenir compte du detrainement et de l'entrainement dans cette
     240! couche
     241! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
     242! avant) a l'alimentation pour avoir un calcul plus propre
    246243!---------------------------------------------------------------------------
    247 !
    248 !test:estimation de ztva_new_est sans entrainement
    249 
    250                Tbef=ztla(ig,l-1)*zpspsk(ig,l)
    251                zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    252                qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
    253                qsatbef=MIN(0.5,qsatbef)
    254                zcor=1./(1.-retv*qsatbef)
    255                qsatbef=qsatbef*zcor
    256                Zsat = (max(0.,zqta(ig,l-1)-qsatbef) .gt. 1.e-10)
    257                if (Zsat) then
    258                qlbef=max(0.,zqta(ig,l-1)-qsatbef)
    259                DT = 0.5*RLvCp*qlbef
    260                do while (abs(DT).gt.DDT0)
    261                  Tbef=Tbef+DT
    262                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    263                  qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
    264                  qsatbef=MIN(0.5,qsatbef)
    265                  zcor=1./(1.-retv*qsatbef)
    266                  qsatbef=qsatbef*zcor
    267                  qlbef=zqta(ig,l-1)-qsatbef
    268 
    269                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    270                  zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    271                  zcor=1./(1.-retv*qsatbef)
    272                  dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
    273                  num=-Tbef+ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*qlbef
    274                  denom=1.+RLvCp*dqsat_dT
    275                  DT=num/denom
    276                enddo
    277                  zqla_est(ig,l) = max(0.,zqta(ig,l-1)-qsatbef)
    278                endif
     244
     245     call thermcell_condens(ngrid,active, &
     246&          zpspsk(:,l),pplev(:,l),ztla(:,l-1),zqta(:,l-1),zqla_est(:,l))
     247
     248    do ig=1,ngrid
     249        if(active(ig)) then
    279250        ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)
     251        zta_est(ig,l)=ztva_est(ig,l)
    280252        ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)
    281         zta_est(ig,l)=ztva_est(ig,l)
    282253        ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)  &
    283254     &      -zqla_est(ig,l))-zqla_est(ig,l))
    … …  
    287258     &                   /(f_star(ig,l)+alim_star(ig,l))**2+  &
    288259     &                   2.*RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    289 !     &                   *1./3. &
    290260     &                   *(zlev(ig,l+1)-zlev(ig,l))
    291261             if (w_est(ig,l+1).lt.0.) then
    292262                w_est(ig,l+1)=zw2(ig,l)
    293263             endif
    294 !
    295 !calcul du detrainement
    296 !=======================
    297 
    298 !CR:on vire les modifs
    299          if (iflag_thermals_ed==0) then
    300 
    301 ! Modifications du calcul du detrainement.
    302 ! Dans la version de la these de Catherine, on passe brusquement
    303 ! de la version seche a la version nuageuse pour le detrainement
    304 ! ce qui peut occasioner des oscillations.
    305 ! dans la nouvelle version, on commence par calculer un detrainement sec.
    306 ! Puis un autre en cas de nuages.
    307 ! Puis on combine les deux lineairement en fonction de la quantite d'eau.
    308 
    309 #define int1d3
    310 !#undef int1d3
    311 #define RIO_TH
    312 #ifdef RIO_TH
    313 !1. Cas non nuageux
    314 ! 1.1 on est sous le zmax_sec et w croit
    315           if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
    316      &       (zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
    317 #ifdef int1d3
    318      &       (zqla_est(ig,l).lt.1.e-10)) then
    319 #else
    320      &       (zqla(ig,l-1).lt.1.e-10)) then
    321 #endif
    322              detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
    323      &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
    324      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
    325      &       /(r_aspect*zmax_sec(ig)))
    326              detr_stara(ig,l)=detr_star(ig,l)
    327 
    328        if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l',ig,l
    329 
    330 ! 1.2 on est sous le zmax_sec et w decroit
    331           else if ((zlev(ig,l+1).lt.zmax_sec(ig)).and.  &
    332 #ifdef int1d3
    333      &            (zqla_est(ig,l).lt.1.e-10)) then
    334 #else
    335      &            (zqla(ig,l-1).lt.1.e-10)) then
    336 #endif
    337              detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
    338      &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
    339      &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
    340      &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
    341      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
    342      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
    343      &       *((zmax_sec(ig)-zlev(ig,l))/  &
    344      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
    345              detr_starb(ig,l)=detr_star(ig,l)
    346 
    347         if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l',ig,l
    348 
     264       endif
     265    enddo
     266
     267!-------------------------------------------------
     268!calcul des taux d'entrainement et de detrainement
     269!-------------------------------------------------
     270
     271     do ig=1,ngrid
     272        if (active(ig)) then
     273          zdz=zlev(ig,l+1)-zlev(ig,l)
     274          zbuoy=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     275          zfact=fact_gamma/(1.+fact_gamma)
     276 
     277! estimation de la fraction couverte par les thermiques
     278          zalpha=f0(ig)*f_star(ig,l)/sqrt(w_est(ig,l))/rhobarz(ig,l)
     279
     280!calcul de la soumission papier
     281          if (1.eq.1) then
     282             fact_epsilon=0.0007
     283!             zfact=0.9/(1.+0.9)
     284             zfact=0.3
     285             fact_gamma=0.7
     286             fact_gamma2=0.6
     287             expa=0.25
     288! Calcul  du taux d'entrainement entr_star (epsilon)
     289           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
     290     &     zbuoy/w_est(ig,l+1) )&
     291!- fact_epsilon/zalpha**0.25  ) &
     292     &     +0.000 )
     293
     294!           entr_star(ig,l)=f_star(ig,l)*zdz * (  1./3 * MAX(0.,  &     
     295!     &     zbuoy/w_est(ig,l+1) - 1./zalpha**0.25  ) &
     296!     &     +0.000 )
     297! Calcul du taux de detrainement detr_star (delta)
     298!           if (zqla_est(ig,l).lt.1.e-10) then
     299!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     300!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
     301!     &     +0.0006 )
     302!           else
     303!           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     304!     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
     305!     &     +0.002 )
     306!           endif
     307           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     308     &     fact_gamma2 * MAX(0., &
     309!fact_epsilon/zalpha**0.25
     310     &      -zbuoy/w_est(ig,l+1) )       &
     311!     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
     312!test
     313     &     +  0.006*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &   
     314     &     +0.0000 )
    349315          else
    350316
    351 ! 1.3 dans les autres cas
    352              detr_star(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
    353      &                      *(zlev(ig,l+1)-zlev(ig,l))
    354              detr_starc(ig,l)=detr_star(ig,l)
    355 
    356         if (prt_level.ge.20) print*,'coucou calcul detr 3 n: ig, l',ig, l
    357              
     317! Calcul  du taux d'entrainement entr_star (epsilon)
     318           entr_star(ig,l)=f_star(ig,l)*zdz * (  zfact * MAX(0.,  &     
     319     &     zbuoy/w_est(ig,l+1) - fact_epsilon  ) &
     320     &     +0.0000 )
     321
     322! Calcul du taux de detrainement detr_star (delta)
     323           detr_star(ig,l)=f_star(ig,l)*zdz * (                           &
     324     &     fact_gamma2 * MAX(0.,fact_epsilon-zbuoy/w_est(ig,l+1) )       &
     325     &     + 0.002*(max(zqta(ig,l-1)-po(ig,l),0.)/po(ig,l))**0.6  &
     326     &     +0.0000 )
     327
    358328          endif
    359 
    360 #else
    361 
    362 ! 1.1 on est sous le zmax_sec et w croit
    363           if ((w_est(ig,l+1).gt.w_est(ig,l)).and.  &
    364      &       (zlev(ig,l+1).lt.zmax_sec(ig)) ) then
    365              detr_star(ig,l)=MAX(0.,(rhobarz(ig,l+1)  &
    366      &       *sqrt(w_est(ig,l+1))*sqrt(l_mix*zlev(ig,l+1))  &
    367      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))*sqrt(l_mix*zlev(ig,l)))  &
    368      &       /(r_aspect*zmax_sec(ig)))
    369 
    370        if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
    371 
    372 ! 1.2 on est sous le zmax_sec et w decroit
    373           else if ((zlev(ig,l+1).lt.zmax_sec(ig)) ) then
    374              detr_star(ig,l)=-f0(ig)*f_star(ig,lmix(ig))  &
    375      &       /(rhobarz(ig,lmix(ig))*wmaxa(ig))*  &
    376      &       (rhobarz(ig,l+1)*sqrt(w_est(ig,l+1))  &
    377      &       *((zmax_sec(ig)-zlev(ig,l+1))/  &
    378      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.  &
    379      &       -rhobarz(ig,l)*sqrt(w_est(ig,l))  &
    380      &       *((zmax_sec(ig)-zlev(ig,l))/  &
    381      &       ((zmax_sec(ig)-zlev(ig,lmix(ig)))))**2.)
    382        if (prt_level.ge.20) print*,'coucou calcul detr 1: ig, l', ig, l
    383 
    384           else
    385              detr_star=0.
    386           endif
    387 
    388 ! 1.3 dans les autres cas
    389           detr_starc(ig,l)=0.002*f0(ig)*f_star(ig,l)  &
    390      &                      *(zlev(ig,l+1)-zlev(ig,l))
    391 
    392           coefc=min(zqla(ig,l-1)/1.e-3,1.)
    393           if (zlev(ig,l+1).ge.zmax_sec(ig)) coefc=1.
    394           coefc=1.
    395 ! il semble qu'il soit important de baser le calcul sur
    396 ! zqla_est(ig,l-1) plutot que sur zqla_est(ig,l)
    397           detr_star(ig,l)=detr_starc(ig,l)*coefc+detr_star(ig,l)*(1.-coefc)
    398 
    399         if (prt_level.ge.20) print*,'coucou calcul detr 2: ig, l', ig, l
    400 
    401 #endif
    402 
    403 
    404         if (prt_level.ge.20) print*,'coucou calcul detr 444: ig, l', ig, l
    405 !IM 730508 beg
    406 !        if(itap.GE.7200) THEN
    407 !         print*,'th_plume ig,l,itap,zqla_est=',ig,l,itap,zqla_est(ig,l)
    408 !        endif
    409 !IM 730508 end
    410          
    411          zqla0(ig,l)=zqla_est(ig,l)
    412          detr_star0(ig,l)=detr_star(ig,l)
    413 !IM 060508 beg
    414 !         if(detr_star(ig,l).GT.1.) THEN
    415 !          print*,'th_plumeBEF ig l detr_star detr_starc coefc',ig,l,detr_star(ig,l) &
    416 !   &      ,detr_starc(ig,l),coefc
    417 !         endif
    418 !IM 060508 end
    419 !IM 160508 beg
    420 !IM 160508       IF (f0(ig).NE.0.) THEN
    421            detr_star(ig,l)=detr_star(ig,l)/f0(ig)
    422 !IM 160508       ELSE IF(detr_star(ig,l).EQ.0.) THEN
    423 !IM 160508        print*,'WARNING1  : th_plume f0=0, detr_star=0: ig, l, itap',ig,l,itap
    424 !IM 160508       ELSE
    425 !IM 160508        print*,'WARNING2  : th_plume f0=0, ig, l, itap, detr_star',ig,l,itap,detr_star(ig,l)
    426 !IM 160508       ENDIF
    427 !IM 160508 end
    428 !IM 060508 beg
    429 !        if(detr_star(ig,l).GT.1.) THEN
    430 !         print*,'th_plumeAFT ig l detr_star f0 1/f0',ig,l,detr_star(ig,l),f0(ig), &
    431 !   &     float(1)/f0(ig)
    432 !        endif
    433 !IM 060508 end
    434         if (prt_level.ge.20) print*,'coucou calcul detr 445: ig, l', ig, l
    435 !
    436 !calcul de entr_star
    437 
    438 ! #undef test2
    439 ! #ifdef test2
    440 ! La version test2 destabilise beaucoup le modele.
    441 ! Il semble donc que ca aide d'avoir un entrainement important sous
    442 ! le nuage.
    443 !         if (zqla_est(ig,l-1).ge.1.e-10.and.l.gt.lalim(ig)) then
    444 !          entr_star(ig,l)=0.4*detr_star(ig,l)
    445 !         else
    446 !          entr_star(ig,l)=0.
    447 !         endif
    448 ! #else
    449 !
    450 ! Deplacement du calcul de entr_star pour eviter d'avoir aussi
    451 ! entr_star > fstar.
    452 ! Redeplacer suite a la transformation du cas detr>f
    453 ! FH
    454 
    455         if (prt_level.ge.20) print*,'coucou calcul detr 446: ig, l', ig, l
    456 #define int1d
    457 !FH 070508 #define int1d4
    458 !#undef int1d4
    459 ! L'option int1d4 correspond au choix dans le cas ou le detrainement
    460 ! devient trop grand.
    461 
    462 #ifdef int1d
    463 
    464 #ifdef int1d4
    465 #else
    466        detr_star(ig,l)=min(detr_star(ig,l),f_star(ig,l))
    467 !FH 070508 plus
    468        detr_star(ig,l)=min(detr_star(ig,l),1.)
    469 #endif
    470 
    471        entr_star(ig,l)=max(0.4*detr_star(ig,l)-alim_star(ig,l),0.)
    472 
    473         if (prt_level.ge.20) print*,'coucou calcul detr 447: ig, l', ig, l
    474 #ifdef int1d4
    475 ! Si le detrainement excede le flux en bas + l'entrainement, le thermique
    476 ! doit disparaitre.
    477        if (detr_star(ig,l)>f_star(ig,l)+entr_star(ig,l)) then
    478           detr_star(ig,l)=f_star(ig,l)+entr_star(ig,l)
    479           f_star(ig,l+1)=0.
    480           linter(ig)=l+1
    481           zw2(ig,l+1)=-1.e-10
    482        endif
    483 #endif
    484 
    485 
    486 #else
    487 
    488         if (prt_level.ge.20) print*,'coucou calcul detr 448: ig, l', ig, l
    489         if(l.gt.lalim(ig)) then
    490          entr_star(ig,l)=0.4*detr_star(ig,l)
    491         else
    492 
    493 ! FH :
    494 ! Cette ligne doit permettre de garantir qu'on a toujours un flux = 1
    495 ! en haut de la couche d'alimentation.
    496 ! A remettre en questoin a la premiere occasion mais ca peut aider a
    497 ! ecrire un code robuste.
    498 ! Que ce soit avec ca ou avec l'ancienne facon de faire (e* = 0 mais
    499 ! d* non nul) on a une discontinuité de e* ou d* en haut de la couche
    500 ! d'alimentation, ce qui n'est pas forcement heureux.
    501 
    502         if (prt_level.ge.20) print*,'coucou calcul detr 449: ig, l', ig, l
    503 #undef pre_int1c
    504 #ifdef pre_int1c
    505          entr_star(ig,l)=max(detr_star(ig,l)-alim_star(ig,l),0.)
    506          detr_star(ig,l)=entr_star(ig,l)
    507 #else
    508          entr_star(ig,l)=0.
    509 #endif
    510 
     329!endif choix du calcul de E* et D*
     330
     331!cr test
     332!           entr_star(ig,l)=entr_star(ig,l)+0.2*detr_star(ig,l)     
     333
     334! Prise en compte de la fraction
     335!      detr_star(ig,l)=detr_star(ig,l)*sqrt(0.01/max(zalpha,1.e-5))
     336
     337! En dessous de lalim, on prend le max de alim_star et entr_star pour
     338! alim_star et 0 sinon
     339        if (l.lt.lalim(ig)) then
     340          alim_star(ig,l)=max(alim_star(ig,l),entr_star(ig,l))
     341          entr_star(ig,l)=0.
    511342        endif
    512343
    513 #endif
    514 
    515         if (prt_level.ge.20) print*,'coucou calcul detr 440: ig, l', ig, l
    516         entr_star1(ig,l)=entr_star(ig,l)
    517         detr_star1(ig,l)=detr_star(ig,l)
    518 !
    519 
    520 #ifdef int1d
    521 #else
    522         if (detr_star(ig,l).gt.f_star(ig,l)) then
    523 
    524 !  Ce test est là entre autres parce qu'on passe par des valeurs
    525 !  delirantes de detr_star.
    526 !  ca vaut sans doute le coup de verifier pourquoi.
    527 
    528            detr_star(ig,l)=f_star(ig,l)
    529 #ifdef pre_int1c
    530            if (l.gt.lalim(ig)+1) then
    531                entr_star(ig,l)=0.
    532                alim_star(ig,l)=0.
    533 ! FH ajout pour forcer a stoper le thermique juste sous le sommet
    534 ! de la couche (voir calcul de finter)
    535                zw2(ig,l+1)=-1.e-10
    536                linter(ig)=l+1
    537             else
    538                entr_star(ig,l)=0.4*detr_star(ig,l)
    539             endif
    540 #else
    541            entr_star(ig,l)=0.4*detr_star(ig,l)
    542 #endif
    543         endif
    544 #endif
    545 
    546       else !l > 2
    547          detr_star(ig,l)=0.
    548          entr_star(ig,l)=0.
    549       endif
    550 
    551         entr_star2(ig,l)=entr_star(ig,l)
    552         detr_star2(ig,l)=detr_star(ig,l)
    553         if (prt_level.ge.20) print*,'coucou calcul detr 450: ig, l', ig, l
    554 
    555        endif  ! iflag_thermals_ed==0
    556 
    557 !CR:nvlle def de entr_star et detr_star
    558       if (iflag_thermals_ed>=1) then
    559 !      if (l.lt.lalim(ig)) then
    560 !      if (l.lt.2) then
    561 !        entr_star(ig,l)=0.
    562 !        detr_star(ig,l)=0.
    563 !      else
    564 !      if (0.001.gt.(RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))/(2.*w_est(ig,l+1)))) then
    565 !         entr_star(ig,l)=0.001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    566 !      else
    567 !         entr_star(ig,l)=  &
    568 !     &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
    569 !     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
    570 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    571 
    572  
    573          entr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &         
    574      &                f_star(ig,l)/(2.*w_est(ig,l+1))        &
    575      &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)   &
    576      &                *(zlev(ig,l+1)-zlev(ig,l))) &
    577      &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    578 
    579         if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
    580             alim_star_tot(ig)=alim_star_tot(ig)+entr_star(ig,l)
    581             lalim(ig)=lmix_bis(ig)
    582             if(prt_level.GE.10) print*,'alim_star_tot',alim_star_tot(ig),entr_star(ig,l)
    583         endif
    584 
    585         if (((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10).and.(l.le.lmix_bis(ig))) then
    586 !        c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
    587          c2(ig,l)=0.001
    588          detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
    589      &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
    590      &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
    591      &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
    592      &                *(zlev(ig,l+1)-zlev(ig,l)))                    &
    593      &                +0.0001*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    594 
    595        else
    596 !         c2(ig,l)=2500000.*zqla_est(ig,l)/(1004.*zta_est(ig,l))
    597           c2(ig,l)=0.003
    598 
    599          detr_star(ig,l)=MAX(0.*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)),  &
    600      &                c2(ig,l)*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l)) &
    601      &                -f_star(ig,l)/(2.*w_est(ig,l+1))       &
    602      &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)       &
    603      &                *(zlev(ig,l+1)-zlev(ig,l))) &
    604      &                +0.0002*f_star(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    605        endif
    606          
    607            
    608 !        detr_star(ig,l)=detr_star(ig,l)*3.
    609 !        if (l.lt.lalim(ig)) then
    610 !          entr_star(ig,l)=0.
    611 !        endif
    612 !        if (l.lt.2) then
    613 !          entr_star(ig,l)=0.
    614 !          detr_star(ig,l)=0.
    615 !        endif
    616 
    617 
    618 !      endif
    619 !      else if ((ztva_est(ig,l)-ztv(ig,l)).gt.1.e-10) then
    620 !      entr_star(ig,l)=MAX(0.,0.8*f_star(ig,l)/(2.*w_est(ig,l+1))        &
    621 !     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))   &
    622 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    623 !      detr_star(ig,l)=0.002*f_star(ig,l)                         &
    624 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    625 !      else
    626 !      entr_star(ig,l)=0.001*f_star(ig,l)                         &
    627 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    628 !      detr_star(ig,l)=MAX(0.,-0.2*f_star(ig,l)/(2.*w_est(ig,l+1))       &
    629 !     &                *RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l))       &
    630 !     &                *(zlev(ig,l+1)-zlev(ig,l))                      &
    631 !     &                +0.002*f_star(ig,l)                             &
    632 !     &                *(zlev(ig,l+1)-zlev(ig,l))
    633 !      endif
    634 
    635       endif   ! iflag_thermals_ed==1
    636 
    637 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    638 ! FH inutile si on conserve comme on l'a fait plus haut entr=detr
    639 ! dans la couche d'alimentation
    640 !pas d entrainement dans la couche alim
    641 !      if ((l.le.lalim(ig))) then
    642 !           entr_star(ig,l)=0.
    643 !      endif
    644 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    645 !
    646 !prise en compte du detrainement et de l entrainement dans le calcul du flux
    647 
     344! Calcul du flux montant normalise
    648345      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    649346     &              -detr_star(ig,l)
    650347
    651 !test sur le signe de f_star
    652         if (prt_level.ge.20) print*,'coucou calcul detr 451: ig, l', ig, l
    653        if (f_star(ig,l+1).gt.1.e-10) then
     348      endif
     349   enddo
     350
    654351!----------------------------------------------------------------------------
    655352!calcul de la vitesse verticale en melangeant Tl et qt du thermique
    656353!---------------------------------------------------------------------------
    657 !
    658        Zsat=.false.
    659        ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
     354   activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
     355   do ig=1,ngrid
     356       if (activetmp(ig)) then
     357           Zsat=.false.
     358           ztla(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    660359     &            (alim_star(ig,l)+entr_star(ig,l))*zthl(ig,l))  &
    661360     &            /(f_star(ig,l+1)+detr_star(ig,l))
    662 !
    663        zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
     361           zqta(ig,l)=(f_star(ig,l)*zqta(ig,l-1)+  &
    664362     &            (alim_star(ig,l)+entr_star(ig,l))*po(ig,l))  &
    665363     &            /(f_star(ig,l+1)+detr_star(ig,l))
    666 ! 
    667                Tbef=ztla(ig,l)*zpspsk(ig,l)
    668                zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    669                qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)               
    670                qsatbef=MIN(0.5,qsatbef)
    671                zcor=1./(1.-retv*qsatbef)
    672                qsatbef=qsatbef*zcor
    673                Zsat = (max(0.,zqta(ig,l)-qsatbef) .gt. 1.e-10)
    674                if (Zsat) then
    675                qlbef=max(0.,zqta(ig,l)-qsatbef)
    676                DT = 0.5*RLvCp*qlbef
    677                do while (abs(DT).gt.DDT0)
    678                  Tbef=Tbef+DT
    679                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    680                  qsatbef= R2ES * FOEEW(Tbef,zdelta)/pplev(ig,l)
    681                  qsatbef=MIN(0.5,qsatbef)
    682                  zcor=1./(1.-retv*qsatbef)
    683                  qsatbef=qsatbef*zcor
    684                  qlbef=zqta(ig,l)-qsatbef
    685 
    686                  zdelta=MAX(0.,SIGN(1.,RTT-Tbef))
    687                  zcvm5=R5LES*(1.-zdelta) + R5IES*zdelta
    688                  zcor=1./(1.-retv*qsatbef)
    689                  dqsat_dT=FOEDE(Tbef,zdelta,zcvm5,qsatbef,zcor)
    690                  num=-Tbef+ztla(ig,l)*zpspsk(ig,l)+RLvCp*qlbef
    691                  denom=1.+RLvCp*dqsat_dT
    692                  DT=num/denom
    693               enddo
    694                  zqla(ig,l) = max(0.,qlbef)
    695               endif
    696 !   
     364
     365        endif
     366    enddo
     367
     368   call thermcell_condens(ngrid,activetmp,zpspsk(:,l),pplev(:,l),ztla(:,l),zqta(:,l),zqla(:,l))
     369
     370
     371   do ig=1,ngrid
     372      if (activetmp(ig)) then
    697373        if (prt_level.ge.20) print*,'coucou calcul detr 4512: ig, l', ig, l
    698374! on ecrit de maniere conservative (sat ou non)
    … …  
    707383!on ecrit zqsat
    708384           zqsatth(ig,l)=qsatbef 
    709 !calcul de vitesse
    710            zw2(ig,l+1)=zw2(ig,l)*  &
    711      &                 ((f_star(ig,l))**2)  &
    712 !  Tests de Catherine
    713 !     &                 /(f_star(ig,l+1)+detr_star(ig,l))**2+             &
    714      &      /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-fact_epsilon))**2+ &
    715      &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    716      &                 *fact_gamma &
    717      &                 *(zlev(ig,l+1)-zlev(ig,l))
    718 !prise en compte des forces de pression que qd flottabilité<0
    719 !              zw2(ig,l+1)=zw2(ig,l)*  &
    720 !     &            1./(1.+2.*entr_star(ig,l)/f_star(ig,l)) + &       
    721 !     &                 (f_star(ig,l))**2 &
    722 !     &                 /(f_star(ig,l)+entr_star(ig,l))**2+ &
    723 !     &                 (f_star(ig,l)-2.*entr_star(ig,l))**2/(f_star(ig,l)+2.*entr_star(ig,l))**2+  &       
    724 !     &                 /(f_star(ig,l+1)+detr_star(ig,l)-entr_star(ig,l)*(1.-2.))**2+ &
    725 !     &                 /(f_star(ig,l)**2+2.*2.*detr_star(ig,l)*f_star(ig,l)+2.*entr_star(ig,l)*f_star(ig,l))+ &
    726 !     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    727 !     &                 *1./3. &
     385
     386!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     387!          zw2(ig,l+1)=&
     388!     &                 zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*(zlev(ig,l+1)-zlev(ig,l))) &
     389!     &                 +2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
     390!     &                 *1./(1.+fact_gamma) &
    728391!     &                 *(zlev(ig,l+1)-zlev(ig,l))
    729          
    730 !        write(30,*),l+1,zw2(ig,l+1)-zw2(ig,l), &
    731 !     &              -2.*entr_star(ig,l)/f_star(ig,l)*zw2(ig,l), &
    732 !     &               2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    733 
    734  
    735 !             zw2(ig,l+1)=zw2(ig,l)*  &
    736 !     &                 (2.-2.*entr_star(ig,l)/f_star(ig,l)) & 
    737 !     &                 -zw2(ig,l-1)+  &       
    738 !     &                 2.*RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    739 !     &                 *1./3. &
    740 !     &                 *(zlev(ig,l+1)-zlev(ig,l))             
    741 
    742             endif
    743         endif
     392!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     393! La meme en plus modulaire :
     394           zbuoy=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     395           zdz=zlev(ig,l+1)-zlev(ig,l)
     396
     397
     398           zeps=(entr_star(ig,l)+alim_star(ig,l))/(f_star(ig,l)*zdz)
     399
     400if (1==0) then
     401           zw2modif=zw2(ig,l)*(1-fact_epsilon/(1.+fact_gamma)*2.*zdz)
     402           zdw2=2.*zbuoy/(1.+fact_gamma)*zdz
     403           zw2(ig,l+1)=zw2modif+zdw2
     404else
     405! Tentative de reecriture de l'equation de w2. A reprendre ...
     406!           zdw2=2*zdz*zbuoy
     407!           zw2modif=zw2(ig,l)*(1.-2.*zdz*(zeps+fact_epsilon))
     408!!!!!      zw2(ig,l+1)=(zw2(ig,l)+zdw2)/(1.+2.*zfactw2(ig,l))
     409!formulation Arnaud
     410!           zw2fact=zbuoy*zalpha**expa/fact_epsilon
     411!           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa*(1+fact_gamma))*zdz) &
     412!      &    +zw2fact
     413           if (zbuoy.gt.1.e-10) then
     414           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma-zfact)
     415           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
     416      &    +zw2fact
     417           else
     418           zw2fact=zbuoy*zalpha**expa/fact_epsilon*(fact_gamma)
     419           zw2(ig,l+1)=(zw2(ig,l)-zw2fact)*exp(-2.*fact_epsilon/(zalpha**expa)*zdz) &
     420      &    +zw2fact
     421
     422           endif
     423
     424endif
     425!           zw2(ig,l+1)=zw2modif+zdw2
     426      endif
     427   enddo
     428
    744429        if (prt_level.ge.20) print*,'coucou calcul detr 460: ig, l',ig, l
    745430!
     431!---------------------------------------------------------------------------
    746432!initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
    747 
     433!---------------------------------------------------------------------------
     434
     435   do ig=1,ngrid
    748436            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    749437!               stop'On tombe sur le cas particulier de thermcell_dry'
    … …  
    753441            endif
    754442
    755 !        if ((zw2(ig,l).gt.0.).and. (zw2(ig,l+1).le.0.)) then
    756443        if (zw2(ig,l+1).lt.0.) then
    757444           linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    … …  
    771458            wmaxa(ig)=wa_moy(ig,l+1)
    772459        endif
    773         enddo
     460   enddo
     461
     462!=========================================================================
     463! FIN DE LA BOUCLE VERTICALE
    774464      enddo
    775 
    776 !on remplace a* par e* ds premiere couche
    777 !      if (iflag_thermals_ed.ge.1) then
    778 !       do ig=1,ngrid
    779 !       do l=2,klev
    780 !          if (l.lt.lalim(ig)) then
    781 !             alim_star(ig,l)=entr_star(ig,l)
    782 !          endif
    783 !       enddo
    784 !       enddo
    785 !       do ig=1,ngrid
    786 !          lalim(ig)=lmix_bis(ig)
    787 !       enddo
    788 !      endif
    789        if (iflag_thermals_ed.ge.1) then
    790           do ig=1,ngrid
    791              do l=2,lalim(ig)
    792                 alim_star(ig,l)=entr_star(ig,l)
    793                 entr_star(ig,l)=0.
    794              enddo
    795            enddo
    796        endif
     465!=========================================================================
     466
     467!on recalcule alim_star_tot
     468       do ig=1,ngrid
     469          alim_star_tot(ig)=0.
     470       enddo
     471       do ig=1,ngrid
     472          do l=1,lalim(ig)-1
     473          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
     474          enddo
     475       enddo
     476       
     477
    797478        if (prt_level.ge.20) print*,'coucou calcul detr 470: ig, l', ig, l
    798479
Note: See TracChangeset for help on using the changeset viewer.