Ignore:
Timestamp:
Apr 29, 2019, 10:07:47 AM (6 years ago)
Author:
aboissinot
Message:

Fix in convadj.F
New version of the thermal plume model (new entrainment and detrainment formulae, alimentation is removed)

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/convadj.F

    r2107 r2127  
    197197            l2 = l2 + 1
    198198            IF (l2 .GT. nlay) EXIT
    199             IF ((zhc(i, l2).LT.zhc(i, l2-1)).and.(l.GT.lmax(i))) THEN
     199            IF ((zhc(i, l2).LT.zhc(i, l2-1)).and.(l2.GT.lmax(i))) THEN
    200200 
    201201!     l2 is the highest level of the unstable column
    … …  
    226226                down = .false.
    227227                IF (l1 .ne. 1) then    !--  and then
    228                   IF ((zhmc.LT.zhc(i, l1-1)).and.(l.GT.lmax(i))) then
     228                  IF ((zhmc.LT.zhc(i, l1-1)).and.(l1.GT.lmax(i))) then
    229229                    down = .true.
    230230                  END IF

  • trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

    r2121 r2127  
    11901190     
    11911191      IF (calltherm) THEN
    1192 ! AB : We need to evaporate ice before calling thermcell_main
     1192         
     1193! AB: We need to evaporate ice before calling thermcell_main.
    11931194         IF (water) THEN
    11941195            CALL evap(ngrid,nlayer,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,zqtherm,zttherm)
    … …  
    11981199         ENDIF
    11991200         
    1200 ! AB:  If a plume stops, the parametrization never look above if somewhere the atmosphere is still unstable!
    1201 !      Maybe it's a good idea to call convective adjustment too.
    1202          CALL thermcell_main(icount, ngrid, nlayer, ptimestep,                   &
    1203                              pplay, pplev, pphi, firstcall,                      &
    1204                              pu, pv, zttherm, zqtherm(:,:,igcm_h2o_vap),         &
    1205                              zdutherm, zdvtherm, zdttherm, zdotherm,             &
    1206                              f0, fm0, entr0, detr0,                              &
     1201! AB: If a plume stops, the parametrization never look above if somewhere the atmosphere is still unstable!
     1202!     Maybe it's a good idea to call convective adjustment too.
     1203         CALL thermcell_main(ngrid, nlayer, nq, ptimestep, firstcall,            &
     1204                             pplay, pplev, pphi, zpopsk,                         &
     1205                             pu, pv, zttherm, zqtherm,                           &
     1206                             zdutherm, zdvtherm, zdttherm, zdqtherm,             &
     1207                             f0, fm0, entr0, detr0, zw2, fraca,                  &
    12071208                             zqta, zqla, ztv, ztva, ztla, zthl, zqsatth,         &
    1208                              zw2, fraca,                                         &
    1209                              lmin, lmix, lalim, lmax,                            &
    1210                              zpopsk, ratqscth, ratqsdiff,                        &
    1211 ! AB : Next variables are only used for diagnoses
    1212                              Ale_bl,Alp_bl,lalim_conv,wght_th,                   &
    1213                              pbl_tke,pctsrf,omega_therm,airephy,                 &
    1214                              zlcl,fraca0,w0,w_conv,                              &
    1215                              therm_tke_max0,env_tke_max0,                        &
    1216                              n2,s2,ale_bl_stat,                                  &
    1217                              therm_tke_max,env_tke_max,                          &
    1218                              alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,         &
    1219                              alp_bl_conv,alp_bl_stat)
     1209                             lmin, lmix, lalim, lmax)
    12201210         
    12211211         DO ig=1,ngrid
    … …  
    12281218         pdt(1:ngrid,1:nlayer) = pdt(1:ngrid,1:nlayer) + zdttherm(1:ngrid,1:nlayer)
    12291219         
    1230          IF(tracer) THEN
    1231             pdq(1:ngrid,1:nlayer,igcm_h2o_vap) = pdq(1:ngrid,1:nlayer,igcm_h2o_vap) + zdotherm(1:ngrid,1:nlayer)
     1220         IF (tracer) THEN
     1221            pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqtherm(1:ngrid,1:nlayer,1:nq)
    12321222         ENDIF
    12331223         
    12341224      ELSE
    12351225         
    1236          lmax(1:ngrid) = 0
     1226         lmax(1:ngrid) = 1
    12371227         
    12381228      ENDIF ! end of 'calltherm'
    … …  
    23582348         CALL send_xios_field('entr',entr0)
    23592349         CALL send_xios_field('detr',detr0)
    2360          CALL send_xios_field('dt_plm',zdttherm)
    2361          CALL send_xios_field('pmax',pmax)
    23622350      ENDIF
     2351     
    23632352      IF (water) THEN
    23642353         CALL send_xios_field('h2o_vap',zq(:,:,igcm_h2o_vap))

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_closure.F90

    r2064 r2127  
    22!
    33!
    4       SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,             &
    5       &                            lalim,alim_star,f_star,                    &
    6       &                            zmax,wmax,f,lev_out)
     4SUBROUTINE thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
     5                             lmax,alim_star,zmax,wmax,f)
    76     
    8 !==============================================================================
    9 !thermcell_closure: fermeture, determination de f
     7!===============================================================================
     8!  Purpose: fermeture, determination de f
    109!
    1110! Modification 7 septembre 2009
    … …  
    1514! l'idee etant que le choix se fasse a l'appel de thermcell_closure
    1615! 3. Vectorisation en mettant les boucles en l l'exterieur avec des if
    17 !==============================================================================
     16!===============================================================================
    1817     
    1918      USE thermcell_mod
    … …  
    2221     
    2322     
    24 !=============================================================================
     23!===============================================================================
    2524! Declaration
    26 !=============================================================================
     25!===============================================================================
    2726     
    28 !      inputs:
    29 !      -------
     27!     Inputs:
     28!     -------
    3029     
    31       INTEGER ngrid,nlay
    32       INTEGER lalim(ngrid)
    33       INTEGER lev_out                           ! niveau pour les print
     30      INTEGER ngrid, nlay
     31      INTEGER lmax(ngrid)
    3432     
    35       REAL alim_star(ngrid,nlay)
    36       REAL f_star(ngrid,nlay+1)
     33      REAL ptimestep
    3734      REAL rho(ngrid,nlay)
    3835      REAL zlev(ngrid,nlay)
     36      REAL alim_star(ngrid,nlay)
    3937      REAL zmax(ngrid)
    4038      REAL wmax(ngrid)
    41       REAL ptimestep
    4239     
    43 !      outputs:
    44 !      --------
     40!     Outputs:
     41!     --------
    4542     
    4643      REAL f(ngrid)
    4744     
    48 !      local:
    49 !      ------
     45!     Local:
     46!     ------
    5047     
    51       INTEGER ig,k
     48      INTEGER ig, l
    5249      INTEGER llmax
    5350     
    … …  
    5653      REAL zdenom(ngrid)
    5754     
    58 !==============================================================================
     55!===============================================================================
    5956! Initialization
    60 !==============================================================================
     57!===============================================================================
    6158     
    6259      alim_star2(:) = 0.
    … …  
    6764      llmax = 1
    6865     
    69 !==============================================================================
     66!===============================================================================
    7067! Closure
    71 !==============================================================================
     68!===============================================================================
    7269     
    73 !------------------------------------------------------------------------------
    74 ! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
    75 !------------------------------------------------------------------------------
     70!-------------------------------------------------------------------------------
     71! Indice vertical max atteint par les thermiques sur le domaine
     72!-------------------------------------------------------------------------------
    7673     
    7774      DO ig=1,ngrid
    78          IF (lalim(ig)>llmax) THEN
    79             llmax = lalim(ig)
     75         IF (lmax(ig).GT.llmax) THEN
     76            llmax = lmax(ig)
    8077         ENDIF
    8178      ENDDO
    8279     
    83 !------------------------------------------------------------------------------
     80!-------------------------------------------------------------------------------
    8481! Calcul des integrales sur la verticale de alim_star et de alim_star^2/(rho dz)
    85 !------------------------------------------------------------------------------
     82!-------------------------------------------------------------------------------
    8683     
    87       DO k=1,llmax-1
     84      DO l=1,llmax-1
    8885         DO ig=1,ngrid
    89             IF (k<lalim(ig)) THEN
    90                alim_star2(ig) = alim_star2(ig) + alim_star(ig,k)**2           &
    91                &              / (rho(ig,k) * (zlev(ig,k+1) - zlev(ig,k)))
    92                alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,k)
     86            IF (l < lmax(ig)) THEN
     87               alim_star2(ig) = alim_star2(ig) + alim_star(ig,l)**2           &
     88               &              / (rho(ig,l) * (zlev(ig,l+1) - zlev(ig,l)))
     89               alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
    9390            ENDIF
    9491         ENDDO
    … …  
    9693     
    9794      DO ig=1,ngrid
    98          IF (alim_star2(ig)>1.e-10) THEN
     95         IF (alim_star2(ig).GT.0.) THEN
    9996            f(ig) = wmax(ig) * alim_star_tot(ig)                              &
    100             &     / (max(500.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
     97            &     / (max(1.,zmax(ig)) * r_aspect_thermals * alim_star2(ig))
     98         ELSE
     99            f(ig) = 0.
    101100         ENDIF
    102101      ENDDO
    103102     
     103     
    104104RETURN
    105105END

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_dq.F90

    r2102 r2127  
    22!
    33!
    4 SUBROUTINE thermcell_dq(ngrid,nlay,impl,ptimestep,fm,entr,detr,masse,         &
    5                         q,dq,qa,lmin,lev_out)
     4SUBROUTINE thermcell_dq(ngrid,nlay,ptimestep,fm,entr,detr,masse,              &
     5                        q,dq,qa,lmin)
    66     
    77     
    8 !==============================================================================
    9 !  Calcul du transport verticale dans la couche limite en presence
    10 !  de "thermiques" explicitement representes
    11 !  calcul du dq/dt une fois qu'on connait les ascendances
    12 !
     8!===============================================================================
     9!  Purpose: Calcul du transport verticale dans la couche limite en presence de
     10!          "thermiques" explicitement representes
     11!           Calcul du dq/dt une fois qu'on connait les ascendances
     12! 
    1313!  Modif 2013/01/04 (FH hourdin@lmd.jussieu.fr)
    14 !  Introduction of an implicit computation of vertical advection in
    15 !  the environment of thermal plumes in thermcell_dq
    16 !  impl =     0 : explicit, 1 : implicit, -1 : old version
    17 !
    18 !  Modif 2018/11/05 (AB alexandre.boissinot@lmd.jussieu.fr)
     14!  Introduction of an implicit computation of vertical advection in the environ-
     15!     ment of thermal plumes in thermcell_dq
    1916! 
    20 !
    21 !==============================================================================
     17!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
     18!     dqimpl = 1 : implicit scheme
     19!     dqimpl = 0 : explicit scheme
     20! 
     21!===============================================================================
    2222     
    2323      USE print_control_mod, ONLY: prt_level
     24      USe thermcell_mod, ONLY: dqimpl
    2425     
    2526      IMPLICIT NONE
    2627     
    2728     
    28 !==============================================================================
     29!===============================================================================
    2930! Declaration
    30 !==============================================================================
     31!===============================================================================
    3132     
    32 !      inputs:
    33 !      -------
     33!     Inputs:
     34!     -------
    3435     
    3536      INTEGER ngrid, nlay
    36       INTEGER impl
    3737      INTEGER lmin(ngrid)
    38       INTEGER lev_out                           ! niveau pour les print
    3938     
    4039      REAL ptimestep
    … …  
    4443      REAL detr(ngrid,nlay)
    4544      REAL q(ngrid,nlay)
     45     
     46!     Outputs:
     47!     --------
     48     
     49      REAL dq(ngrid,nlay)
    4650      REAL qa(ngrid,nlay)
    4751     
    48 !      outputs:
    49 !      --------
    50      
    51       REAL dq(ngrid,nlay)
    52      
    53 !      local:
    54 !      ------
     52!     Local:
     53!     ------
    5554     
    5655      INTEGER ig, l
    … …  
    6362      REAL zzm
    6463     
    65       CHARACTER (LEN=20) :: modname
    66       CHARACTER (LEN=80) :: abort_message
     64!===============================================================================
     65! Initialization
     66!===============================================================================
    6767     
    68 !==============================================================================
    69 ! Initialization
    70 !==============================================================================
     68      qold(:,:) = q(:,:)
    7169     
    72       qold = q
     70!===============================================================================
     71! Tracer variation computation
     72!===============================================================================
    7373     
    74       modname = 'thermcell_dq'
    75      
    76       IF (impl.lt.0) THEN
    77          print *, 'OLD EXPLICIT SCHEME IS USED'
    78          CALL thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,               &
    79          &                   masse,q,dq,qa,lev_out)
    80          RETURN
    81       ENDIF
    82      
    83 !------------------------------------------------------------------------------
     74!-------------------------------------------------------------------------------
    8475! CFL criterion computation for advection in downdraft
    85 !------------------------------------------------------------------------------
     76!-------------------------------------------------------------------------------
    8677     
    8778      cfl = 0.
    … …  
    10394               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
    10495               print *, 'fm-', fm(ig,l-1)
    105                abort_message = 'entr dt > m, 1st'
    106                CALL abort_physic(modname,abort_message,1)
     96               CALL abort
    10797            ENDIF
    10898         ENDDO
    10999      ENDDO
    110100     
    111 !------------------------------------------------------------------------------
     101!-------------------------------------------------------------------------------
    112102! Computation of tracer concentrations in the ascending plume
    113 !------------------------------------------------------------------------------
     103!-------------------------------------------------------------------------------
    114104     
    115105      DO ig=1,ngrid
    … …  
    121111      DO ig=1,ngrid
    122112         DO l=lmin(ig)+1,nlay
    123             if ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-5*masse(ig,l)) then
     113            IF ((fm(ig,l+1)+detr(ig,l))*ptimestep.gt.1.e-6*masse(ig,l)) THEN
    124114               qa(ig,l) = (fm(ig,l) * qa(ig,l-1) + entr(ig,l) * q(ig,l))      &
    125115               &        / (fm(ig,l+1) + detr(ig,l))
    126             else
     116            ELSE
    127117               qa(ig,l) = q(ig,l)
    128             endif
     118            ENDIF
    129119           
    130120!            IF (qa(ig,l).lt.0.) THEN
    … …  
    140130      ENDDO
    141131     
    142 !------------------------------------------------------------------------------
    143 ! Plume vertical flux of q
    144 !------------------------------------------------------------------------------
     132!-------------------------------------------------------------------------------
     133! Plume vertical flux of tracer
     134!-------------------------------------------------------------------------------
    145135     
    146136      DO l=2,nlay-1
    … …  
    151141      fqa(:,nlay) = 0.
    152142     
    153 !------------------------------------------------------------------------------
     143!-------------------------------------------------------------------------------
    154144! Trace species evolution
    155 !------------------------------------------------------------------------------
     145!-------------------------------------------------------------------------------
    156146     
    157       IF (impl==0) THEN
    158          do l=1,nlay-1
     147      IF (dqimpl==0) THEN
     148         DO l=1,nlay-1
    159149            q(:,l) = q(:,l) + (fqa(:,l) - fqa(:,l+1) - fm(:,l) * q(:,l)       &
    160150            &      + fm(:,l+1) * q(:,l+1)) * ptimestep / masse(:,l)
    161          enddo
    162       ELSE
    163          do l=nlay-1,1,-1
     151         ENDDO
     152      ELSEIF (dqimpl==1) THEN
     153         DO l=nlay-1,1,-1
    164154            q(:,l) = ( q(:,l) + ptimestep / masse(:,l)                        &
    165155            &      * ( fqa(:,l) - fqa(:,l+1) + fm(:,l+1) * q(:,l+1) ) )       &
    166156            &      / ( 1. + fm(:,l) * ptimestep / masse(:,l) )
    167157         ENDDO
     158      ELSE
     159         print *, 'ERROR: No corresponding scheme for mixing computations!'
     160         print *, '       dqimpl must be equal to 1, 0 or -1 but'
     161         print *, 'dqimpl =', dqimpl
     162         call abort
    168163      ENDIF
    169164     
    170 !==============================================================================
     165!===============================================================================
    171166! Tendencies
    172 !==============================================================================
     167!===============================================================================
    173168     
    174169      DO l=1,nlay
    … …  
    182177RETURN
    183178END
    184 
    185 
    186 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    187 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    188 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    189 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    190 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    192 
    193 
    194 Subroutine thermcell_dq_o(ngrid,nlay,impl,ptimestep,fm,entr,masse,            &
    195                           q,dq,qa,lev_out)
    196      
    197      
    198 !==============================================================================
    199 !
    200 !   Calcul du transport verticale dans la couche limite en presence
    201 !   de "thermiques" explicitement representes
    202 !   calcul du dq/dt une fois qu'on connait les ascendances
    203 !
    204 !==============================================================================
    205      
    206       USE print_control_mod, ONLY: prt_level
    207      
    208       implicit none
    209      
    210      
    211 !==============================================================================
    212 ! Declaration
    213 !==============================================================================
    214      
    215       integer ngrid,nlay,impl
    216 
    217       real ptimestep
    218       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    219       real entr(ngrid,nlay)
    220       real q(ngrid,nlay)
    221       real dq(ngrid,nlay)
    222       integer lev_out                           ! niveau pour les print
    223 
    224       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
    225 
    226       real zzm
    227 
    228       integer ig,l
    229       real cfl
    230 
    231       real qold(ngrid,nlay)
    232       real ztimestep
    233       integer niter,iter
    234       CHARACTER (LEN=20) :: modname='thermcell_dq'
    235       CHARACTER (LEN=80) :: abort_message
    236 
    237 
    238 
    239 ! Calcul du critere CFL pour l'advection dans la subsidence
    240       cfl = 0.
    241       do l=1,nlay
    242          do ig=1,ngrid
    243             zzm=masse(ig,l)/ptimestep
    244             cfl=max(cfl,fm(ig,l)/zzm)
    245             if (entr(ig,l).gt.zzm) then
    246                print*,'entr*dt>m,2',l,entr(ig,l)*ptimestep,masse(ig,l)
    247                abort_message = 'entr dt > m, 2nd'
    248                CALL abort_physic (modname,abort_message,1)
    249             endif
    250          enddo
    251       enddo
    252      
    253 !IM 090508     print*,'CFL CFL CFL CFL ',cfl
    254      
    255 #undef CFL
    256 #ifdef CFL
    257 ! On subdivise le calcul en niter pas de temps.
    258       niter=int(cfl)+1
    259 #else
    260       niter=1
    261 #endif
    262      
    263       ztimestep=ptimestep/niter
    264       qold=q
    265      
    266       do iter=1,niter
    267          if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    268          
    269 !   calcul du detrainement
    270          do l=1,nlay
    271             do ig=1,ngrid
    272                detr(ig,l)=fm(ig,l)-fm(ig,l+1)+entr(ig,l)
    273 !               print*,'Q2 DQ ',detr(ig,l),fm(ig,l),entr(ig,l)
    274 !test
    275                if (detr(ig,l).lt.0.) then
    276                   entr(ig,l)=entr(ig,l)-detr(ig,l)
    277                   detr(ig,l)=0.
    278 !                  print*,'detr2<0!!!','ig=',ig,'l=',l,'f=',fm(ig,l),
    279 !     s            'f+1=',fm(ig,l+1),'e=',entr(ig,l),'d=',detr(ig,l)
    280                endif
    281                
    282 !               if (fm(ig,l+1).lt.0.) then
    283 !                  print*,'fm2<0!!!'
    284 !               endif
    285                
    286 !               if (entr(ig,l).lt.0.) then
    287 !                  print*,'entr2<0!!!'
    288 !               endif
    289             enddo
    290          enddo
    291          
    292 ! calcul de la valeur dans les ascendances
    293         do ig=1,ngrid
    294             qa(ig,1)=q(ig,1)
    295          enddo
    296          
    297          do l=2,nlay
    298             do ig=1,ngrid
    299                if ((fm(ig,l+1)+detr(ig,l))*ztimestep.gt.1.e-5*masse(ig,l)) then
    300                   qa(ig,l) = (fm(ig,l)*qa(ig,l-1)+entr(ig,l)*q(ig,l))         &
    301                   &        / (fm(ig,l+1)+detr(ig,l))
    302                else
    303                   qa(ig,l)=q(ig,l)
    304                endif
    305                
    306                if (qa(ig,l).lt.0.) then
    307 !                  print*,'qa<0!!!'
    308                endif
    309                
    310                if (q(ig,l).lt.0.) then
    311 !                  print*,'q<0!!!'
    312                endif
    313             enddo
    314          enddo
    315          
    316 ! Calcul du flux subsident
    317          
    318          do l=2,nlay
    319             do ig=1,ngrid
    320 #undef centre
    321 #ifdef centre
    322                 wqd(ig,l)=fm(ig,l)*0.5*(q(ig,l-1)+q(ig,l))
    323 #else
    324 
    325 #define plusqueun
    326 #ifdef plusqueun
    327 ! Schema avec advection sur plus qu'une maille.
    328                zzm=masse(ig,l)/ztimestep
    329                if (fm(ig,l)>zzm) then
    330                  wqd(ig,l)=zzm*q(ig,l)+(fm(ig,l)-zzm)*q(ig,l+1)
    331                else
    332                   wqd(ig,l)=fm(ig,l)*q(ig,l)
    333                endif
    334 #else
    335                wqd(ig,l)=fm(ig,l)*q(ig,l)
    336 #endif
    337 #endif
    338 
    339                if (wqd(ig,l).lt.0.) then
    340 !                  print*,'wqd<0!!!'
    341                endif
    342             enddo
    343          enddo
    344          
    345          do ig=1,ngrid
    346             wqd(ig,1)=0.
    347             wqd(ig,nlay+1)=0.
    348          enddo
    349          
    350 ! Calcul des tendances
    351          do l=1,nlay
    352             do ig=1,ngrid
    353                q(ig,l)=q(ig,l)+(detr(ig,l)*qa(ig,l)-entr(ig,l)*q(ig,l)           &
    354                &        -wqd(ig,l)+wqd(ig,l+1))                                  &
    355                &        *ztimestep/masse(ig,l)
    356 !               if (dq(ig,l).lt.0.) then
    357 !                  print*,'dq<0!!!'
    358 !               endif
    359            enddo
    360          enddo
    361       enddo
    362      
    363 ! Calcul des tendances
    364       do l=1,nlay
    365          do ig=1,ngrid
    366             dq(ig,l) = (q(ig,l) - qold(ig,l)) / ptimestep
    367             q(ig,l) = qold(ig,l)
    368          enddo
    369       enddo
    370      
    371      
    372 return
    373 end

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_env.F90

    r2071 r2127  
    22!
    33!
    4       SUBROUTINE thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,  &
    5       &           pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,pqsat,lev_out)
    6 
    7 !==============================================================================
    8 ! thermcell_env: calcule les caracteristiques de l environnement
    9 !                necessaires au calcul des proprietes dans le thermique
    10 !==============================================================================
    11 
     4SUBROUTINE thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
     5                         zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
     6     
     7     
     8!===============================================================================
     9!  Purpose: calcul des caracteristiques de l'environnement necessaires au calcul
     10!           des proprietes dans le thermique.
     11! 
     12!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
     13! 
     14!  Nota Bene
     15!     ql   means "non-gaseous water mass mixing ratio" (liquid and solid)
     16!     qv   means "vapor mass mixing ratio"
     17!     qt   means "total water mass mixing ratio"
     18!     TP   means "potential temperature"
     19!     TRPV means "virtual potential temperature with latent heat release" 
     20!     TPV  means "virtual potential temperature"
     21!     TR   means "temperature with latent heat release"
     22!===============================================================================
     23     
    1224      USE print_control_mod, ONLY: prt_level
    1325      USE thermcell_mod, ONLY: RKAPPA
    1426      USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
    15       USE planete_mod, ONLY: preff
     27      USE tracer_h, ONLY: igcm_h2o_vap, igcm_h2o_ice
     28      USE callkeys_mod, ONLY: water
    1629     
    17       IMPLICIT NONE   
     30      IMPLICIT NONE
    1831     
    19 !==============================================================================
     32     
     33!===============================================================================
    2034! Declaration
    21 !==============================================================================
     35!===============================================================================
    2236     
    23 !      inputs:
    24 !      -------
     37!     Inputs:
     38!     -------
    2539     
    26       INTEGER ngrid
    27       INTEGER nlay
    28       INTEGER lev_out                           ! niveau pour les print
     40      INTEGER ngrid, nlay, nq
    2941     
    30       REAL po(ngrid,nlay)                       ! large scale water
    31       REAL pt(ngrid,nlay)                       ! large scale temperature
    32       REAL pu(ngrid,nlay)                       ! large scale zonal wind
    33       REAL pv(ngrid,nlay)                       ! large scale meridional wind
    34       REAL pplay(ngrid,nlay)                    ! layers pressure
    35       REAL pplev(ngrid,nlay+1)                  ! levels pressure
     42      REAL pq(ngrid,nlay,nq)                    ! Large scale water
     43      REAL pt(ngrid,nlay)                       ! Large scale temperature
     44      REAL pu(ngrid,nlay)                       ! Large scale zonal wind
     45      REAL pv(ngrid,nlay)                       ! Large scale meridional wind
     46      REAL pplay(ngrid,nlay)                    ! Layers pressure
     47      REAL pplev(ngrid,nlay+1)                  ! Levels pressure
     48      REAL zpopsk(ngrid,nlay)                   ! Exner function
    3649     
    37 !      outputs:
    38 !      --------
     50!     Outputs:
     51!     --------
    3952     
    40       REAL zo(ngrid,nlay)                       ! qt  environment
    41       REAL zl(ngrid,nlay)                       ! ql  environment
    42       REAL zh(ngrid,nlay)                       ! T   environment
    43       REAL ztv(ngrid,nlay)                      ! TV  environment
    44       REAL zthl(ngrid,nlay)                     ! TPV environment
    45       REAL zu(ngrid,nlay)                       ! u   environment
    46       REAL zv(ngrid,nlay)                       ! v   environment
     53      REAL zqt(ngrid,nlay)                      ! q    environment
     54      REAL zql(ngrid,nlay)                      ! q    environment
     55      REAL zt(ngrid,nlay)                       ! T    environment
     56      REAL ztv(ngrid,nlay)                      ! TRPV environment
     57      REAL zhl(ngrid,nlay)                      ! TP   environment
     58      REAL zu(ngrid,nlay)                       ! u    environment
     59      REAL zv(ngrid,nlay)                       ! v    environment
     60      REAL zqs(ngrid,nlay)                      ! qsat environment
    4761     
    48       REAL zpspsk(ngrid,nlay)                   ! Exner function
    49       REAL pqsat(ngrid,nlay)                    ! saturation vapor pressure
     62!     Local:
     63!     ------
    5064     
    51 !      local:
    52 !      ------
     65      INTEGER ig, l
    5366     
    54       INTEGER ig, ll
     67      REAL psat                                 ! Dummy argument for Psat_water()
    5568     
    56       REAL dummy
     69!===============================================================================
     70! Initialization
     71!===============================================================================
    5772     
    58 !==============================================================================
    59 ! Initialization
    60 !==============================================================================
     73      zu(:,:) = pu(:,:)
     74      zv(:,:) = pv(:,:)
    6175     
    62 !------------------------------------------------------------------------------
    63 ! Calcul des caracteristiques de l'environnement
    64 !------------------------------------------------------------------------------
    65       DO ll=1,nlay
    66          DO ig=1,ngrid
    67             zo(ig,ll) = po(ig,ll)
    68             zl(ig,ll) = 0.
    69             zh(ig,ll) = pt(ig,ll)
     76      zhl(:,:) = pt(:,:) / zpopsk(:,:)
     77     
     78      zqt(:,:) = pq(:,:,igcm_h2o_vap)
     79     
     80!===============================================================================
     81! Condensation and latent heat release
     82!===============================================================================
     83     
     84      IF (water) THEN
     85         
     86         DO l=1,nlay
     87            DO ig=1,ngrid
     88               CALL Psat_water(pt(ig,l), pplev(ig,l), psat, zqs(ig,l))
     89            ENDDO
    7090         ENDDO
    71       ENDDO
    72      
    73 !==============================================================================
    74 ! Computation of qsat and condensation
    75 !==============================================================================
    76 
    77       DO ll=1,nlay
    78          DO ig=1,ngrid
    79             CALL Psat_water(pt(ig,ll), pplev(ig,ll), dummy, pqsat(ig,ll))
     91         
     92         DO l=1,nlay
     93            DO ig=1,ngrid
     94               zql(ig,l) = max(0.,pq(ig,l,igcm_h2o_vap) - zqs(ig,l))
     95               zt(ig,l) = pt(ig,l) + RLvCp * zql(ig,l)
     96               ztv(ig,l) = zt(ig,l) / zpopsk(ig,l)                            &
     97               &         * (1. + RETV * (zqt(ig,l)-zql(ig,l)) - zql(ig,l))
     98            ENDDO
    8099         ENDDO
    81       ENDDO
    82      
    83       DO ll=1,nlay
    84          DO ig=1,ngrid
    85             zl(ig,ll) = max(0.,po(ig,ll) - pqsat(ig,ll))                      ! zl   set to ql   env
    86             zh(ig,ll) = pt(ig,ll) + RLvCp * zl(ig,ll)                         ! zh   set to TR   env
    87             zo(ig,ll) = po(ig,ll) - zl(ig,ll)                                 ! zo   set to qv   env
    88          ENDDO
    89       ENDDO
    90      
    91       DO ll=1,nlay
    92          DO ig=1,ngrid
    93              zpspsk(ig,ll) = (pplay(ig,ll)/preff)**RKAPPA
    94              zu(ig,ll) = pu(ig,ll)
    95              zv(ig,ll) = pv(ig,ll)
    96              
    97 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    98 ! AB : WARNING, zh is no longer the potentiel temperature.
    99 !      It is now the temperature. How awful!
    100 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    101             ztv(ig,ll) = zh(ig,ll)/zpspsk(ig,ll)                              ! ztv  set to TRP  env
    102             ztv(ig,ll) = ztv(ig,ll)*(1.+RETV*(zo(ig,ll))-zl(ig,ll))           ! ztv  set to TRPV env
    103             zthl(ig,ll) = pt(ig,ll)/zpspsk(ig,ll)                             ! zthl set to TP   env
    104            
    105          ENDDO
    106       ENDDO
     100         
     101      ELSE
     102         
     103         zt(:,:) = pt(:,:)
     104         ztv(:,:) = pt(:,:) / zpopsk(:,:)
     105         
     106      ENDIF
    107107     
    108108     

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_flux.F90

    r2103 r2127  
    33!
    44SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
    5                           lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
    6                           f,rhobarz,zlev,zw2,fm,entr,                         &
    7                           detr,zqla,lev_out,lunout1,igout)
    8      
    9      
    10 !==============================================================================
    11 ! thermcell_flux: deduction des flux
    12 !==============================================================================
     5                          lalim,lmin,lmax,entr_star,detr_star,                &
     6                          f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
     7     
     8     
     9!===============================================================================
     10!  Purpose: deduction des flux
     11! 
     12!  Modif 2019/04 (AB alexandre.boissinot@lmd.jussieu.fr)
     13! 
     14!===============================================================================
    1315     
    1416      USE print_control_mod, ONLY: prt_level
    … …  
    1820     
    1921     
    20 !==============================================================================
     22!===============================================================================
    2123! Declaration
    22 !==============================================================================
    23      
    24 !      inputs:
    25 !      -------
    26      
    27       INTEGER ngrid
    28       INTEGER nlay
    29       INTEGER igout
    30       INTEGER lev_out
    31       INTEGER lunout1
     24!===============================================================================
     25     
     26!     Inputs:
     27!     -------
     28     
     29      INTEGER ngrid, nlay
    3230      INTEGER lmin(ngrid)
    3331      INTEGER lmax(ngrid)
    3432      INTEGER lalim(ngrid)
    3533     
    36       REAL alim_star(ngrid,nlay)
     34! AB : I remove alimentation, which is in fact included in entr_star. EN COURS
     35!      REAL alim_star(ngrid,nlay)
    3736      REAL entr_star(ngrid,nlay)
    3837      REAL detr_star(ngrid,nlay)
    … …  
    4645      REAL zmax(ngrid)
    4746     
    48 !      outputs:
    49 !      --------     
     47!     Outputs:
     48!     --------     
    5049     
    5150      REAL entr(ngrid,nlay)
    … …  
    5352      REAL fm(ngrid,nlay+1)
    5453     
    55 !      local:
    56 !      ------
    57      
    58       INTEGER ig,l,k
    59       INTEGER lout
    60      
    61       REAL zfm                            ! mass flux such as updraft fraction is equal to its maximal authorized value alphamax
    62       REAL f_old                          ! save initial value of mass flux
    63       REAL eee0                           ! save initial value of entrainment
    64       REAL ddd0                           ! save initial value of detrainment
     54!     Local:
     55!     ------
     56     
     57      INTEGER ig, l, k
     58      INTEGER igout, lout                 ! Error grid point number and level
     59     
     60      REAL zfm                            ! Mass flux such as updraft fraction is equal to its maximal authorized value alphamax
     61      REAL f_old                          ! Save initial value of mass flux
     62      REAL eee0                           ! Save initial value of entrainment
     63      REAL ddd0                           ! Save initial value of detrainment
    6564      REAL eee                            ! eee0 - layer mass * maximal authorized mass fraction / dt
    6665      REAL ddd                            ! ddd0 - eee
    67       REAL zzz                            ! temporary variable set to mass flux
     66      REAL zzz                            ! Temporary variable set to mass flux
    6867      REAL fact
    6968      REAL test
    … …  
    8382      CHARACTER (len=80) :: abort_message
    8483     
    85 !==============================================================================
     84!===============================================================================
    8685! Initialization
    87 !==============================================================================
     86!===============================================================================
    8887     
    8988      ncorecfm1 = 0
    … …  
    10099      fm(:,:)   = 0.
    101100     
    102       IF (prt_level.ge.10) THEN
    103          write(lunout1,*) 'Dans thermcell_flux 0'
    104          write(lunout1,*) 'flux base ', f(igout)
    105          write(lunout1,*) 'lmax ', lmax(igout)
    106          write(lunout1,*) 'lalim ', lalim(igout)
    107          write(lunout1,*) 'ig= ', igout
    108          write(lunout1,*) ' l E*    A*     D*  '
    109          write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),                  &
    110          &                             alim_star(igout,l),detr_star(igout,l), &
    111          &                             l=1,lmax(igout))
    112       ENDIF
    113      
    114 !==============================================================================
     101!===============================================================================
    115102! Calcul de l'entrainement, du detrainement et du flux de masse
    116 !==============================================================================
    117      
    118 !------------------------------------------------------------------------------
     103!===============================================================================
     104     
     105!-------------------------------------------------------------------------------
    119106! Multiplication par la norme issue de la relation de fermeture
    120 !------------------------------------------------------------------------------
     107!-------------------------------------------------------------------------------
    121108     
    122109      DO l=1,nlay
    123          entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
     110! AB : I remove alimentation, which is in fact included in entr_star. EN COURS
     111!         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
     112         entr(:,l) = f(:) * entr_star(:,l)
    124113         detr(:,l) = f(:) * detr_star(:,l)
    125114      ENDDO
    126115     
    127 !------------------------------------------------------------------------------
     116! AB : temporary test
     117      DO l=1,nlay
     118         DO ig=1,ngrid
     119            IF (detr(ig,l)<0.) THEN
     120               print *, 'ERROR: detr is negative in thermcell_flux!'
     121               print *, 'ig,l', ig, l
     122               print *, 'detr_star', detr_star(ig,l)
     123               print *, 'detr', detr(ig,l)
     124            ENDIF
     125         ENDDO
     126      ENDDO
     127     
     128!-------------------------------------------------------------------------------
    128129! Calcul du flux de masse
    129 !------------------------------------------------------------------------------
     130!-------------------------------------------------------------------------------
    130131     
    131132      DO l=1,nlay
    … …  
    144145      ENDDO
    145146     
    146 !==============================================================================
     147!===============================================================================
    147148! Tests de validite
    148 !==============================================================================
     149!===============================================================================
    149150     
    150151      DO l=1,nlay
    151152         
    152 !------------------------------------------------------------------------------
     153!-------------------------------------------------------------------------------
    153154! Verification de la positivite du flux de masse entrant
    154 !------------------------------------------------------------------------------
    155          
    156 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     155!-------------------------------------------------------------------------------
     156         
     157!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    157158! AB : I remove the correction and replace it by an uncompromising test.
    158159!      According to the previous derivations, we MUST have positive mass flux
    … …  
    160161!      Then the only value which can be negative is the mass flux at the top of
    161162!      the plume. However, this value was set to zero a few lines above.
    162 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     163!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    163164         
    164165         labort_physic=.false.
    … …  
    177178            print *, 'lmin,lmax', lmin(igout), lmax(igout)
    178179            print *, '-------------------------------'
    179             print *, 'fm+', fm(igout,lout+1)
    180             print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
    181             print *, 'fm ', fm(igout,lout)
    182             print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
    183             print *, 'fm-', fm(igout,lout-1)
     180            DO k=nlay,1,-1
     181               print *, 'fm ', fm(igout,k+1)
     182               print *, 'entr,detr', entr(igout,k), detr(igout,k)
     183            ENDDO
     184            print *, 'fm-', fm(igout,1)
     185            print *, '-------------------------------'
    184186            abort_message = 'Negative incoming fm in thermcell_flux'
    185187            CALL abort_physic(modname,abort_message,1)
    186188         ENDIF
    187189         
    188 !------------------------------------------------------------------------------
     190!-------------------------------------------------------------------------------
    189191! Test entrainment positif
    190 !------------------------------------------------------------------------------
    191          
    192 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     192!-------------------------------------------------------------------------------
     193         
     194!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    193195! AB : According to the previous derivations, we MUST have positive entrainment
    194196!      and detrainment everywhere!
    195197!      Indeed, they are set to max between zero and a computed value.
    196198!      Then tests are uncompromising.
    197 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    198          
    199          labort_physic=.false.
     199!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    200200         
    201201         DO ig=1,ngrid
    … …  
    212212            print *, 'lmin,lmax', lmin(igout), lmax(igout)
    213213            print *, '-------------------------------'
    214             print *, 'fm+', fm(igout,lout+1)
    215             print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
    216             print *, 'fm ', fm(igout,lout)
    217             print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
    218             print *, 'fm-', fm(igout,lout-1)
     214            DO k=nlay,1,-1
     215               print *, 'fm ', fm(igout,k+1)
     216               print *, 'entr,detr', entr(igout,k), detr(igout,k)
     217            ENDDO
     218            print *, 'fm-', fm(igout,1)
     219            print *, '-------------------------------'
    219220            abort_message = 'Negative entrainment in thermcell_flux'
    220221            CALL abort_physic(modname,abort_message,1)
    221222         ENDIF
    222223         
    223 !------------------------------------------------------------------------------
     224!-------------------------------------------------------------------------------
    224225! Test detrainment positif
    225 !------------------------------------------------------------------------------
     226!-------------------------------------------------------------------------------
    226227         
    227228         DO ig=1,ngrid
    … …  
    238239            print *, 'lmin,lmax', lmin(igout), lmax(igout)
    239240            print *, '-------------------------------'
    240             print *, 'fm+', fm(igout,lout+1)
    241             print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
    242             print *, 'fm ', fm(igout,lout)
    243             print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
    244             print *, 'fm-', fm(igout,lout-1)
     241            DO k=nlay,1,-1
     242               print *, 'fm ', fm(igout,k+1)
     243               print *, 'entr,detr', entr(igout,k), detr(igout,k)
     244            ENDDO
     245            print *, 'fm-', fm(igout,1)
     246            print *, '-------------------------------'
    245247            abort_message = 'Negative detrainment in thermcell_flux'
    246248            CALL abort_physic(modname,abort_message,1)
    247249         ENDIF
    248250         
    249 !------------------------------------------------------------------------------
     251!-------------------------------------------------------------------------------
    250252! Test sur fraca constante ou croissante au-dessus de lalim
    251 !------------------------------------------------------------------------------
     253!-------------------------------------------------------------------------------
    252254         
    253255! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
    … …  
    267269         ENDIF
    268270         
    269 !------------------------------------------------------------------------------
     271!-------------------------------------------------------------------------------
    270272! Test sur flux de masse constant ou decroissant au-dessus de lalim
    271 !------------------------------------------------------------------------------
     273!-------------------------------------------------------------------------------
    272274         
    273275! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
    … …  
    283285         ENDIF
    284286         
    285 !------------------------------------------------------------------------------
     287!-------------------------------------------------------------------------------
    286288! Test detrainment inferieur au flux de masse
    287 !------------------------------------------------------------------------------
    288          
    289 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     289!-------------------------------------------------------------------------------
     290         
     291!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    290292! AB : Even if fm has no negative value, it can be lesser than detr.
    291293!      That's not suitable because when we will mix the plume with the
    … …  
    294296!      otherwise we get a negative outgoing mass flux (cf. above).
    295297!      That is why we correct the detrainment as follows.
    296 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     298!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    297299         
    298300         DO ig=1,ngrid
    … …  
    314316               ncorecfm6  = ncorecfm6 + 1
    315317               
    316 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     318!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    317319! Dans le cas ou on est au dessus de la couche d'alimentation et que le
    318320! detrainement est plus fort que le flux de masse, on stope le thermique.
    … …  
    320322!
    321323! AB : Do we have to stop the plume here?
    322 ! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are
    323 !      set to zero above this level!
    324 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     324! AB : Attention, if lmax is modified, be sure that fm, zw2, entr, detr are set
     325!      to zero above this level!
     326!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    325327!               IF (l.gt.lalim(ig)) THEN
    326328!                  lmax(ig)=l
    … …  
    331333!               ENDIF
    332334            ENDIF
    333 !           
     335           
    334336!            IF (l.gt.lmax(ig)) THEN
    335337!               detr(ig,l) = 0.
    … …  
    360362!         ENDIF
    361363         
    362 !------------------------------------------------------------------------------
     364!-------------------------------------------------------------------------------
    363365! Test fraction masse entrainee inferieure a fomass_max
    364 !------------------------------------------------------------------------------
    365          
    366 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     366!-------------------------------------------------------------------------------
     367         
     368!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    367369! AB : Entrainment is bigger than the maximal authorized value.
    368370!      If we consider that the excess entrainement is in fact plume air which
    … …  
    374376!      detrainment and mass flux profiles such as we do not exceed the maximal
    375377!      authorized entrained mass.
    376 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     378!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    377379         
    378380         labort_physic=.false.
    … …  
    432434            print *, '   ptimestep  :', ptimestep
    433435            print *, 'norm  :', f(igout)
    434             print *, 'alim* :', alim_star(igout,lout)
     436!            print *, 'alim* :', alim_star(igout,lout)
    435437            print *, 'entr* :', entr_star(igout,lout)
    436438            print *, '--------------------'
    … …  
    447449         ENDIF
    448450         
    449 !------------------------------------------------------------------------------
     451!-------------------------------------------------------------------------------
    450452! Test fraction couverte inferieure a alphamax
    451 !------------------------------------------------------------------------------
    452          
    453 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     453!-------------------------------------------------------------------------------
     454         
     455!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    454456! FH : Partie a revisiter.
    455457!      Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1
    … …  
    461463!      semble de toutes facons deraisonable d'avoir des fractions de 1...
    462464!      Ci-dessous, et dans l'etat actuel, le premier des deux if est sans doute inutile. 
    463 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
     465!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
    464466         
    465467         DO ig=1,ngrid
    466468            zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alphamax
     469           
    467470            IF (fm(ig,l+1) .gt. zfm) THEN
    468471               f_old = fm(ig,l+1)
    … …  
    470473               detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
    471474               
    472 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     475!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    473476! AB : The previous change doesn't observe the equation df/dz = e - d.
    474477!      To avoid this issue, what is better to do? I choose to increase the
    475478!      entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and
    476479!      there are never any problems.
    477 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     480!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    478481               IF (l.lt.lmax(ig)) THEN
    479482                  entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
    … …  
    498501      ENDDO
    499502     
    500 !------------------------------------------------------------------------------
     503!-------------------------------------------------------------------------------
    501504! We check if fm, entr and detr vanish correctly in layer lmax
    502 !------------------------------------------------------------------------------
    503      
    504 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     505!-------------------------------------------------------------------------------
     506     
     507!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    505508! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
    506509!      Theoretically, we always get lmax >= lmin >= linf > 0
    507 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     510!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    508511      DO ig=1,ngrid
    509512         IF (lmax(ig).gt.0) THEN
    … …  
    514517      ENDDO
    515518     
    516 !------------------------------------------------------------------------------
     519!-------------------------------------------------------------------------------
    517520! Impression des bidouilles utilisees pour corriger les problemes
    518 !------------------------------------------------------------------------------
     521!-------------------------------------------------------------------------------
    519522     
    520523      IF (prt_level.ge.1) THEN
    … …  
    553556      ENDIF
    554557     
    555 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     558!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    556559! AB : temporary test added to check the validity of eq. df/dz = e - d
    557 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     560!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    558561!      DO l=1,nlay
    559562!         DO ig=1,ngrid

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_height.F90

    r2101 r2127  
    22!
    33!
    4       SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,          &
    5       &                           zw2,zlev,lmax,zmax,zmix,                    &
    6       &                           wmax,f_star,lev_out)
    7      
    8      
    9 !==============================================================================
    10 ! thermcell_height: calcul des caracteristiques du thermique: zmax,wmax,zmix
    11 !==============================================================================
     4SUBROUTINE thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
     5                            zlev,lmax,zmax,zmix,wmax,f_star)
     6     
     7     
     8!===============================================================================
     9!  Purpose: calcul des caracteristiques du thermique: zmax,wmax,zmix
     10!===============================================================================
    1211     
    1312      USE thermcell_mod
    … …  
    1716     
    1817     
    19 !==============================================================================
     18!===============================================================================
    2019! Declaration
    21 !==============================================================================
    22      
    23 !   inputs:
    24 !   -------
    25      
    26       INTEGER ngrid,nlay
    27       INTEGER lalim(ngrid)
     20!===============================================================================
     21     
     22!     Inputs:
     23!     -------
     24     
     25      INTEGER ngrid, nlay
    2826      INTEGER lmin(ngrid)
    29       INTEGER lev_out                           ! niveau pour les print
    3027     
    3128      REAL zlev(ngrid,nlay+1)
    3229      REAL f_star(ngrid,nlay+1)
    3330     
    34 !   outputs:
    35 !   --------
     31!     Outputs:
     32!     --------
    3633     
    3734      INTEGER lmix(ngrid)
    … …  
    3936     
    4037      REAL linter(ngrid)
    41       REAL zmax(ngrid)
    42       REAL zmix(ngrid)
    43       REAL wmax(ngrid)
    44       REAL zw2(ngrid,nlay+1)
    45      
    46 !   local:
    47 !   ------
    48      
    49       INTEGER ig,l
     38      REAL zmax(ngrid)                          ! maximal vertical speed
     39      REAL zmix(ngrid)                          ! altitude of maximal vertical speed
     40      REAL wmax(ngrid)                          ! maximal vertical speed
     41      REAL zw2(ngrid,nlay+1)                    ! square vertical speed (becomes its squere root)
     42     
     43!     Local:
     44!     ------
     45     
     46      INTEGER ig, l
    5047     
    5148      REAL num(ngrid)
    … …  
    5451      REAL zlev2(ngrid,nlay+1)
    5552     
    56 !==============================================================================
     53!===============================================================================
    5754! Initialization
    58 !==============================================================================
    59      
    60       DO ig=1,ngrid
    61          lmax(ig) = lalim(ig)
     55!===============================================================================
     56     
     57      DO ig=1,ngrid
     58         lmax(ig) = lmin(ig)
    6259         lmix(ig) = lmin(ig)
    6360      ENDDO
    … …  
    7673      ENDDO
    7774     
    78 !==============================================================================
     75!===============================================================================
    7976! Calcul de la hauteur max du thermique
    80 !==============================================================================
    81      
    82 !------------------------------------------------------------------------------
     77!===============================================================================
     78     
     79!-------------------------------------------------------------------------------
    8380! Calcul de lmax
    84 !------------------------------------------------------------------------------
     81!-------------------------------------------------------------------------------
    8582     
    8683      DO ig=1,ngrid
    8784         DO l=nlay,lmin(ig)+1,-1
    88             IF (zw2(ig,l).le.0.) THEN
     85            IF (zw2(ig,l).le.0..or.f_star(ig,l).le.0.) THEN
    8986               lmax(ig) = l - 1
    90             ELSEIF (f_star(ig,l).le.0.) THEN
    91                lmax(ig) = l - 1
    9287            ENDIF
    9388         ENDDO
    9489      ENDDO
    9590     
    96 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     91!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    9792! On traite le cas particulier qu'il faudrait eviter ou le thermique
    9893! atteind le haut du modele ...
    99 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    100      
    101       DO ig=1,ngrid
    102          IF (zw2(ig,nlay).gt.1.e-10) THEN
    103             print *, 'WARNING: a thermal plume reaches the highest level!'
    104             lmax(ig)=nlay
    105          ENDIF
    106       ENDDO
    107      
    108 !------------------------------------------------------------------------------
     94!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     95      DO ig=1,ngrid
     96         IF (zw2(ig,nlay).gt.0.) THEN
     97            print *, 'WARNING: a thermal plume reaches the highest layer!'
     98            print *, 'ig,l', ig, nlay
     99            print *, 'zw2', zw2(ig,nlay)
     100            lmax(ig) = nlay
     101         ENDIF
     102      ENDDO
     103     
     104!-------------------------------------------------------------------------------
    109105! Calcul de zmax continu via le linter
    110 !------------------------------------------------------------------------------
    111      
    112 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     106!-------------------------------------------------------------------------------
     107     
     108!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    113109! AB : lmin=lmax means the plume is not active and then zw2=0 at each level.
    114110!      Otherwise we have lmax>lmin.
    115 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     111!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    116112     
    117113      DO ig=1,ngrid
    … …  
    125121      ENDDO
    126122     
    127 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    128 ! AB : zmax is equal to zlevinter-zlev(lmin) because lmin can be different
     123!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     124! AB : zmax is now equal to zlevinter-zlev(lmin) because lmin can be different
    129125!      from 1.
    130 !      We have to substract this level because zmax must be the plume height
    131 !      (to derive the good closure), not the plume highest altitude.
    132 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    133      
    134       DO  ig=1,ngrid
     126!      We have to substract this level because zmax must be the plume height (to
     127!      derive the good closure), not the plume highest altitude.
     128!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     129     
     130      DO ig=1,ngrid
    135131         zlevinter(ig) = (zlev(ig,lmax(ig)+1)-zlev(ig,lmax(ig)))*linter(ig)   &
    136132         &             + zlev(ig,lmax(ig))-lmax(ig)*(zlev(ig,lmax(ig)+1)      &
    … …  
    139135      ENDDO
    140136     
    141 !==============================================================================
     137!===============================================================================
    142138! Calcul de wmax et du niveau d'inversion.
    143 !==============================================================================
    144      
    145 !------------------------------------------------------------------------------
     139!===============================================================================
     140     
     141!-------------------------------------------------------------------------------
    146142! Calcul de lmix et wmax
    147 !------------------------------------------------------------------------------
     143!-------------------------------------------------------------------------------
    148144     
    149145      DO l=1,nlay
    … …  
    152148               IF (zw2(ig,l).lt.0.) THEN
    153149                  print *, 'ERROR: zw2 has negative value(s)!'
     150                  print *, 'ig,l', ig, l
     151                  print *, 'zw2', zw2(ig,l)
    154152                  CALL abort
    155153               ENDIF
    156154               
    157 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     155!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    158156! AB : WARNING zw2 becomes its square root!
    159 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     157!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    160158               zw2(ig,l) = sqrt(zw2(ig,l))
    161159               
    … …  
    170168      ENDDO
    171169     
    172 !------------------------------------------------------------------------------
     170!-------------------------------------------------------------------------------
    173171! Calcul de zmix continu (profil parabolique des vitesses)
    174 !------------------------------------------------------------------------------
    175      
    176 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     172!-------------------------------------------------------------------------------
     173     
     174!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    177175! AB : We substract zlev(lmin) in zmix formula because we have to
    178176!      compare it with zmax which is zlev(linter)-zlev(lmin).
    179 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     177!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    180178     
    181179      DO ig=1,ngrid
    … …  
    197195            ELSE
    198196               zmix(ig) = zlev(ig,lmix(ig)) - zlev(ig,lmin(ig))
    199                print*,'WARNING: problematic zmix value!'
     197               print *, 'WARNING: problematic zmix value!'
    200198            ENDIF
    201199         ELSE
    … …  
    204202      ENDDO
    205203     
    206 !==============================================================================
     204!===============================================================================
    207205! Verification zmax > zmix
    208 !==============================================================================
     206!===============================================================================
    209207     
    210208      DO ig=1,ngrid
    211209         IF ((zmax(ig)-zmix(ig)).lt.0.) THEN
    212            
    213210            zmix(ig) = 0.9 * zmax(ig)
    214            
    215             IF (prt_level.ge.100) THEN
    216                print *, 'WARNING: we have zmix > zmax!'
    217             ENDIF
    218          ENDIF
    219       ENDDO
    220      
    221 !------------------------------------------------------------------------------
     211            print *, 'WARNING: we have zmix > zmax!'
     212            print *, 'zmax,zmix', zmax(ig), zmix(ig)
     213         ENDIF
     214      ENDDO
     215     
     216!-------------------------------------------------------------------------------
    222217! Calcul du nouveau lmix correspondant
    223 !------------------------------------------------------------------------------
     218!-------------------------------------------------------------------------------
    224219     
    225220      DO ig=1,ngrid
    226221         DO l=1,nlay
    227222            IF (zmix(ig).ge.zlev2(ig,l).and.zmix(ig).lt.zlev2(ig,l+1)) THEN
    228                lmix(ig)=l
     223               lmix(ig) = l
    229224            ENDIF
    230225         ENDDO

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_main.F90

    r2104 r2127  
    22!
    33!
    4 SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep,                          &
    5                           pplay,pplev,pphi,firstcall,                         &
    6                           pu,pv,pt,po,                                        &
    7                           pduadj,pdvadj,pdtadj,pdoadj,                        &
    8                           f0,fm0,entr0,detr0,                                 &
    9                           zqta,zqla,ztv,ztva,ztla,zthl,zqsatth,               &
    10                           zw2,fraca,                                          &
    11                           lmin,lmix,lalim,lmax,                               &
    12                           zpspsk,ratqscth,ratqsdiff,                          &
    13                           Ale_bl,Alp_bl,lalim_conv,wght_th,                   &
    14 !!! nrlmd le 10/04/2012
    15                           pbl_tke,pctsrf,omega,airephy,                       &
    16                           zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,  &
    17                           n2,s2,ale_bl_stat,                                  &
    18                           therm_tke_max,env_tke_max,                          &
    19                           alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,         &
    20                           alp_bl_conv,alp_bl_stat)
    21 !!! fin nrlmd le 10/04/2012
    22      
    23      
    24 !==============================================================================
     4SUBROUTINE thermcell_main(ngrid,nlay,nq,ptimestep,firstcall,                  &
     5                          pplay,pplev,pphi,zpopsk,                            &
     6                          pu,pv,pt,pq,                                        &
     7                          pduadj,pdvadj,pdtadj,pdqadj,                        &
     8                          f0,fm0,entr0,detr0,zw2,fraca,                       &
     9                          zqta,zqla,ztv,ztva,zhla,zhl,zqsa,                   &
     10                          lmin,lmix,lalim,lmax)
     11     
     12     
     13!===============================================================================
    2514!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
    2615!   Version du 09.02.07
    … …  
    4433!    Introduction of an implicit computation of vertical advection in
    4534!    the environment of thermal plumes in thermcell_dq
    46 !    impl =     0 : explicit, 1 : implicit, -1 : old version
     35!    impl = 0 : explicit ; impl = 1 : implicit ; impl =-1 : old version
    4736!    controled by iflag_thermals =
    4837!       15, 16 run with impl=-1 : numerical convergence with NPv3
    4938!       17, 18 run with impl=1  : more stable
    5039!    15 and 17 correspond to the activation of the stratocumulus "bidouille"
    51 !
    52 !==============================================================================
     40!
     41! Major changes 2018-19 (AB alexandre.boissinot@lmd.jussieu.fr)
     42!    New detr and entre formulae (no longer alimentation)
     43!    lmin can be greater than 1
     44!    Mix every tracer (EN COURS)
     45!    Old version of thermcell_dq is removed
     46!    Alternative version thermcell_dv2 is removed
     47!
     48!===============================================================================
    5349     
    5450      USE thermcell_mod
     51      USE tracer_h, ONLY: igcm_h2o_vap
    5552      USE print_control_mod, ONLY: lunout, prt_level
    5653     
    … …  
    5855     
    5956     
    60 !==============================================================================
     57!===============================================================================
    6158! Declaration
    62 !==============================================================================
    63      
    64 !   inputs:
    65 !   -------
    66      
    67       INTEGER itap
    68       INTEGER ngrid
    69       INTEGER nlay
     59!===============================================================================
     60     
     61!     Inputs:
     62!     -------
     63     
     64      INTEGER ngrid, nlay, nq
    7065     
    7166      REAL ptimestep
    72       REAL pplay(ngrid,nlay)
    73       REAL pplev(ngrid,nlay+1)
    74       REAL pphi(ngrid,nlay)
    75      
    76       REAL pt(ngrid,nlay)                       !
    77       REAL pu(ngrid,nlay)                       !
    78       REAL pv(ngrid,nlay)                       !
    79       REAL po(ngrid,nlay)                       !
     67      REAL pplay(ngrid,nlay)                    ! Layer pressure
     68      REAL pplev(ngrid,nlay+1)                  ! Level pressure
     69      REAL pphi(ngrid,nlay)                     ! Geopotential
     70     
     71      REAL pu(ngrid,nlay)                       ! Zonal wind
     72      REAL pv(ngrid,nlay)                       ! Meridional wind
     73      REAL pt(ngrid,nlay)                       ! Temperature
     74      REAL pq(ngrid,nlay,nq)                    ! Tracers mass mixing ratio
    8075     
    8176      LOGICAL firstcall
    8277     
    83 !   outputs:
    84 !   --------
    85      
    86       REAL pdtadj(ngrid,nlay)                   ! t convective variations
     78!     Outputs:
     79!     --------
     80     
    8781      REAL pduadj(ngrid,nlay)                   ! u convective variations
    8882      REAL pdvadj(ngrid,nlay)                   ! v convective variations
    89       REAL pdoadj(ngrid,nlay)                   ! water convective variations
    90      
     83      REAL pdtadj(ngrid,nlay)                   ! t convective variations
     84      REAL pdqadj(ngrid,nlay,nq)                ! q convective variations
     85     
     86      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
    9187      REAL fm0(ngrid,nlay+1)                    ! mass flux      (after possible time relaxation)
    9288      REAL entr0(ngrid,nlay)                    ! entrainment    (after possible time relaxation)
    9389      REAL detr0(ngrid,nlay)                    ! detrainment    (after possible time relaxation)
    94       REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
    95      
    96 !   local:
    97 !   ------
    98      
    99       INTEGER, save :: dvimpl=1
    100 !$OMP THREADPRIVATE(dvimpl)
    101      
    102       INTEGER, save :: dqimpl=-1
    103 !$OMP THREADPRIVATE(dqimpl)
    104      
    105       INTEGER, SAVE :: igout=1
    106 !$OMP THREADPRIVATE(igout)
    107      
    108       INTEGER, SAVE :: lunout1=6
    109 !$OMP THREADPRIVATE(lunout1)
    110      
    111       INTEGER, SAVE :: lev_out=10
    112 !$OMP THREADPRIVATE(lev_out)
    113      
    114       INTEGER ig,k,l,ll,ierr
    115       INTEGER lmix_bis(ngrid)
    116       INTEGER lmax(ngrid)                       !
    117       INTEGER lmin(ngrid)                       !
    118       INTEGER lalim(ngrid)                      !
    119       INTEGER lmix(ngrid)                       !
    120      
    121       REAL linter(ngrid)
    122       REAL zmix(ngrid)
    123       REAL zmax(ngrid)
    124       REAL zw2(ngrid,nlay+1)
    125       REAL zw_est(ngrid,nlay+1)
    126       REAL zmax_sec(ngrid)
    127      
    128       REAL zlay(ngrid,nlay)                     ! layers altitude
    129       REAL zlev(ngrid,nlay+1)                   ! levels altitude
    130       REAL rho(ngrid,nlay)                      ! layers density
    131       REAL rhobarz(ngrid,nlay)                  ! levels density
    132       REAL deltaz(ngrid,nlay)                   ! layers heigth
    133       REAL masse(ngrid,nlay)                    ! layers mass
    134       REAL zpspsk(ngrid,nlay)                   ! Exner function
    135      
    136       REAL zu(ngrid,nlay)                       ! environment zonal wind
    137       REAL zv(ngrid,nlay)                       ! environment meridional wind
    138       REAL zo(ngrid,nlay)                       ! environment water tracer
    139       REAL zva(ngrid,nlay)                      ! plume zonal wind
    140       REAL zua(ngrid,nlay)                      ! plume meridional wind
    141       REAL zoa(ngrid,nlay)                      ! plume water tracer
    142      
    143       REAL zt(ngrid,nlay)                       ! T    environment
    144       REAL zh(ngrid,nlay)                       ! T,TP environment
    145       REAL zthl(ngrid,nlay)                     ! TP   environment
    146       REAL ztv(ngrid,nlay)                      ! TPV  environment ?
    147       REAL zl(ngrid,nlay)                       ! ql   environment
    148      
    149       REAL zta(ngrid,nlay)                      !
    150       REAL zha(ngrid,nlay)                      ! TRPV plume
    151       REAL ztla(ngrid,nlay)                     ! TP   plume
    152       REAL ztva(ngrid,nlay)                     ! TRPV plume
    153       REAL ztva_est(ngrid,nlay)                 ! TRPV plume (temporary)
     90     
     91!     Local:
     92!     ------
     93     
     94      INTEGER ig, k, l, iq
     95      INTEGER lmax(ngrid)                       ! Highest layer reached by the plume
     96      INTEGER lmix(ngrid)                       ! Layer in which plume vertical speed is maximal
     97      INTEGER lmin(ngrid)                       ! First unstable layer
     98     
     99      REAL zlay(ngrid,nlay)                     ! Layers altitudes
     100      REAL zlev(ngrid,nlay+1)                   ! Levels altitudes
     101      REAL rho(ngrid,nlay)                      ! Layers densities
     102      REAL rhobarz(ngrid,nlay)                  ! Levels densities
     103      REAL masse(ngrid,nlay)                    ! Layers masses
     104      REAL zpopsk(ngrid,nlay)                   ! Exner function
     105     
     106      REAL zu(ngrid,nlay)                       ! u    environment
     107      REAL zv(ngrid,nlay)                       ! v    environment
     108      REAL zt(ngrid,nlay)                       ! TR   environment
     109      REAL zqt(ngrid,nlay)                      ! qt   environment
     110      REAL zql(ngrid,nlay)                      ! ql   environment
     111      REAL zhl(ngrid,nlay)                      ! TP   environment
     112      REAL ztv(ngrid,nlay)                      ! TRPV environment
     113      REAL zqs(ngrid,nlay)                      ! qsat environment
     114     
     115      REAL zua(ngrid,nlay)                      ! u    plume
     116      REAL zva(ngrid,nlay)                      ! v    plume
     117      REAL zta(ngrid,nlay)                      ! TR   plume
    154118      REAL zqla(ngrid,nlay)                     ! qv   plume
    155119      REAL zqta(ngrid,nlay)                     ! qt   plume
    156      
    157       REAL wmax(ngrid)                          ! maximal vertical speed
    158       REAL wmax_tmp(ngrid)                      ! temporary maximal vertical speed
    159       REAL wmax_sec(ngrid)                      ! maximal vertical speed if dry case
    160      
    161       REAL fraca(ngrid,nlay+1)                  ! updraft fraction
    162       REAL f_star(ngrid,nlay+1)                 ! normalized mass flux
    163       REAL entr_star(ngrid,nlay)                ! normalized entrainment
    164       REAL detr_star(ngrid,nlay)                ! normalized detrainment
    165       REAL alim_star_tot(ngrid)                 ! integrated alimentation
    166       REAL alim_star(ngrid,nlay)                ! normalized alimentation
    167       REAL alim_star_clos(ngrid,nlay)           ! closure alimentation
    168      
    169       REAL fm(ngrid,nlay+1)                     ! mass flux
    170       REAL entr(ngrid,nlay)                     ! entrainment
    171       REAL detr(ngrid,nlay)                     ! detrainment
    172       REAL f(ngrid)                             ! mass flux norm
    173      
    174       REAL zdthladj(ngrid,nlay)                 !
    175       REAL lambda                               ! time relaxation coefficent
    176      
    177       REAL zsortie(ngrid,nlay)
    178       REAL zsortie1d(ngrid)
    179       REAL susqr2pi, Reuler
    180       REAL zf
    181       REAL zf2
    182       REAL thetath2(ngrid,nlay)
    183       REAL wth2(ngrid,nlay)
    184       REAL wth3(ngrid,nlay)
    185       REAL q2(ngrid,nlay)
    186 ! FH : probleme de dimensionnement avec l'allocation dynamique
    187 !     common/comtherm/thetath2,wth2
    188       real wq(ngrid,nlay)
    189       real wthl(ngrid,nlay)
    190       real wthv(ngrid,nlay)
    191       real ratqscth(ngrid,nlay)
    192       real var
    193       real vardiff
    194       real ratqsdiff(ngrid,nlay)
    195 ! niveau de condensation
    196       integer nivcon(ngrid)
    197       real zcon(ngrid)
    198       real CHI
    199       real zcon2(ngrid)
    200       real pcon(ngrid)
    201       real zqsat(ngrid,nlay)
    202       real zqsatth(ngrid,nlay)
    203       real zlevinter(ngrid)
    204       real seuil
    205      
    206 !!! nrlmd le 10/04/2012
    207      
    208 !------Entrees
    209       real pbl_tke(ngrid,nlay+1,nbsrf)
    210       real pctsrf(ngrid,nbsrf)
    211       real omega(ngrid,nlay)
    212       real airephy(ngrid)
    213 !------Sorties
    214       real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
    215       real therm_tke_max0(ngrid)
    216       real env_tke_max0(ngrid)
    217       real n2(ngrid),s2(ngrid)
    218       real ale_bl_stat(ngrid)
    219       real therm_tke_max(ngrid,nlay)
    220       real env_tke_max(ngrid,nlay)
    221       real alp_bl_det(ngrid)
    222       real alp_bl_fluct_m(ngrid)
    223       real alp_bl_fluct_tke(ngrid)
    224       real alp_bl_conv(ngrid)
    225       real alp_bl_stat(ngrid)
    226 !------Local
    227       integer nsrf
    228       real rhobarz0(ngrid)                      ! Densite au LCL
    229       logical ok_lcl(ngrid)                     ! Existence du LCL des thermiques
    230       integer klcl(ngrid)                       ! Niveau du LCL
    231       real interp(ngrid)                        ! Coef d'interpolation pour le LCL
    232 !------Triggering
    233       real,parameter :: Su=4e4                  ! Surface unite: celle d'un updraft elementaire
    234       real,parameter :: hcoef=1.                ! Coefficient directeur pour le calcul de s2
    235       real,parameter :: hmincoef=0.3            ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2 (hmincoef=0.3)
    236       real,parameter :: eps1=0.3                ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
    237       real hmin(ngrid)                          ! Ordonnee a l'origine pour le calcul de s2
    238       real zmax_moy(ngrid)                      ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
    239       real,parameter :: zmax_moy_coef=0.33
    240       real depth(ngrid)                         ! Epaisseur moyenne du cumulus
    241       real w_max(ngrid)                         ! Vitesse max statistique
    242       real s_max(ngrid)
    243 !------Closure
    244       real pbl_tke_max(ngrid,nlay)              ! Profil de TKE moyenne
    245       real pbl_tke_max0(ngrid)                  ! TKE moyenne au LCL
    246       real w_ls(ngrid,nlay)                     ! Vitesse verticale grande echelle (m/s)
    247       real,parameter :: coef_m=1.               ! On considere un rendement pour alp_bl_fluct_m
    248       real,parameter :: coef_tke=1.             ! On considere un rendement pour alp_bl_fluct_tke
    249      
    250 !!! fin nrlmd le 10/04/2012
    251      
    252 ! Nouvelles variables pour la convection
    253       integer lalim_conv(ngrid)
    254       integer n_int(ngrid)
    255       real Ale_bl(ngrid)
    256       real Alp_bl(ngrid)
    257       real alp_int(ngrid)
    258       real dp_int(ngrid),zdp
    259       real ale_int(ngrid)
    260       real fm_tot(ngrid)
    261       real wght_th(ngrid,nlay)
    262      
    263       CHARACTER*2 str2
    264       CHARACTER*10 str10
     120      REAL zhla(ngrid,nlay)                     ! TP   plume
     121      REAL ztva(ngrid,nlay)                     ! TRPV plume
     122      REAL zqsa(ngrid,nlay)                     ! qsat plume
     123     
     124      REAL zqa(ngrid,nlay,nq)                   ! q    plume (ql=0, qv=qt)
     125     
     126      REAL linter(ngrid)                        ! Level (continuous) of maximal vertical speed
     127      REAL zmix(ngrid)                          ! Altitude of maximal vertical speed
     128      REAL zmax(ngrid)                          ! Maximal altitude reached by the plume
     129      REAL wmax(ngrid)                          ! Maximal vertical speed
     130      REAL zw2(ngrid,nlay+1)                    ! Plume vertical speed
     131     
     132      REAL fraca(ngrid,nlay+1)                  ! Updraft fraction
     133     
     134      REAL f_star(ngrid,nlay+1)                 ! Normalized mass flux
     135      REAL entr_star(ngrid,nlay)                ! Normalized entrainment
     136      REAL detr_star(ngrid,nlay)                ! Normalized detrainment
     137     
     138      REAL f(ngrid)                             ! Mass flux norm
     139      REAL fm(ngrid,nlay+1)                     ! Mass flux
     140      REAL entr(ngrid,nlay)                     ! Entrainment
     141      REAL detr(ngrid,nlay)                     ! Detrainment
     142     
     143      REAL lambda                               ! Time relaxation coefficent
     144     
     145      REAL zdthladj(ngrid,nlay)                 ! Potential temperature variations
     146      REAL dummy(ngrid,nlay)                    ! Dummy argument for thermcell_dq()
    265147     
    266148      CHARACTER (len=20) :: modname='thermcell_main'
    267149      CHARACTER (len=80) :: abort_message
    268150     
    269       EXTERNAL SCOPY
    270      
    271 !==============================================================================
     151! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
     152      INTEGER lalim(ngrid)                      ! Highest alimentation level
     153      REAL alim_star(ngrid,nlay)                ! Normalized alimentation
     154      REAL alim_star_tot(ngrid)                 ! Integrated alimentation
     155      REAL alim_star_clos(ngrid,nlay)           ! Closure alimentation
     156     
     157!===============================================================================
    272158! Initialization
    273 !==============================================================================
    274      
    275       seuil = 0.25
     159!===============================================================================
    276160     
    277161      IF (firstcall) THEN
    278          IF (iflag_thermals==15.or.iflag_thermals==16) THEN
    279             dvimpl = 0
    280             dqimpl = -1
    281          ELSE
    282             dvimpl = 1
    283             dqimpl = 1
    284          ENDIF
    285          
    286162         fm0(:,:) = 0.
    287163         entr0(:,:) = 0.
    … …  
    289165      ENDIF
    290166     
     167      f_star(:,:) = 0.
     168      entr_star(:,:) = 0.
     169      detr_star(:,:) = 0.
     170     
     171      f(:) = 0.
     172     
    291173      fm(:,:) = 0.
    292174      entr(:,:) = 0.
    293175      detr(:,:) = 0.
    294       f(:) = 0.
    295      
     176     
     177      lmax(:) = 1
     178      lmix(:) = 1
     179      lmin(:) = 1
     180     
     181      pduadj(:,:) = 0.0
     182      pdvadj(:,:) = 0.0
     183      pdtadj(:,:) = 0.0
     184      pdqadj(:,:,:) = 0.0
     185     
     186! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
     187      alim_star(:,:) = 0.
     188      alim_star_tot(:) = 0.
     189      alim_star_clos(:,:) = 0.
     190      lalim(:) = 1
     191     
     192! AB: Careful, hard-coded value from Earth tuned version of the thermal plume model!
    296193      DO ig=1,ngrid
    297194         f0(ig) = max(f0(ig), 1.e-2)
    … …  
    304201      ENDIF
    305202     
    306       wmax_tmp(:) = 0.
    307      
    308 !------------------------------------------------------------------------------
     203!===============================================================================
     204! Environment settings
     205!===============================================================================
     206     
     207!-------------------------------------------------------------------------------
    309208! Calcul de T,q,ql a partir de Tl et qt dans l environnement
    310 !------------------------------------------------------------------------------
    311      
    312       CALL thermcell_env(ngrid,nlay,po,pt,pu,pv,pplay,                        &
    313       &                  pplev,zo,zh,zl,ztv,zthl,zu,zv,zpspsk,zqsat,lev_out)
    314      
    315 !------------------------------------------------------------------------------
    316 !
    317 !                       --------------------
    318 !
    319 !
    320 !                       + + + + + + + + + + +
    321 !
    322 !
    323 !  wa, fraca, wd, fracd --------------------   zlev(2), rhobarz
    324 !  wh,wt,wo ...
    325 !
    326 !                       + + + + + + + + + + +  zh,zu,zv,zo,rho
    327 !
    328 !
    329 !                       --------------------   zlev(1)
    330 !                       \\\\\\\\\\\\\\\\\\\\
    331 !
    332 !------------------------------------------------------------------------------
    333 ! Calcul des altitudes des couches
    334 !------------------------------------------------------------------------------
     209!-------------------------------------------------------------------------------
     210     
     211      CALL thermcell_env(ngrid,nlay,nq,pq,pt,pu,pv,pplay,pplev,               &
     212      &                  zqt,zql,zt,ztv,zhl,zu,zv,zpopsk,zqs)
     213     
     214!-------------------------------------------------------------------------------
     215! Levels and layers altitudes
     216!-------------------------------------------------------------------------------
    335217     
    336218      DO l=2,nlay
    … …  
    345227      ENDDO
    346228     
    347 !------------------------------------------------------------------------------
    348 ! Calcul de l'epaisseur des couches
    349 !------------------------------------------------------------------------------
    350      
    351       DO l=1,nlay
    352          deltaz(:,l) = zlev(:,l+1)-zlev(:,l)
    353       ENDDO
    354      
    355 !------------------------------------------------------------------------------
    356 ! Calcul des densites
    357 !------------------------------------------------------------------------------
    358      
    359       rho(:,:) = pplay(:,:) / (zpspsk(:,:) * RD * ztv(:,:))
     229!-------------------------------------------------------------------------------
     230! Levels and layers densities
     231!-------------------------------------------------------------------------------
     232     
     233      rho(:,:) = pplay(:,:) / (zpopsk(:,:) * RD * ztv(:,:))
    360234     
    361235      IF (prt_level.ge.10) THEN
    … …  
    369243      ENDDO
    370244     
    371 !------------------------------------------------------------------------------
    372 ! Calcul de la masse
    373 !------------------------------------------------------------------------------
     245!-------------------------------------------------------------------------------
     246! Layers masses
     247!-------------------------------------------------------------------------------
    374248     
    375249      DO l=1,nlay
    … …  
    377251      ENDDO
    378252     
    379       IF (prt_level.ge.1) print *, 'thermcell_main apres initialisation'
    380      
    381 !------------------------------------------------------------------------------
    382 !             
    383 !             /|\
    384 !    --------  |  F_k+1 -------   
    385 !                              ----> D_k
    386 !             /|\              <---- E_k , A_k
    387 !    --------  |  F_k ---------
    388 !                              ----> D_k-1
    389 !                              <---- E_k-1 , A_k-1
    390 !
    391 !
    392 !
    393 !
    394 !
    395 !    ---------------------------
    396 !
    397 !    ----- F_lmax+1=0 ----------         \
    398 !            lmax     (zmax)              |
    399 !    ---------------------------          |
    400 !                                         |
    401 !    ---------------------------          |
    402 !                                         |
    403 !    ---------------------------          |
    404 !                                         |
    405 !    ---------------------------          |
    406 !                                         |
    407 !    ---------------------------          |
    408 !                                         |  E
    409 !    ---------------------------          |  D
    410 !                                         |
    411 !    ---------------------------          |
    412 !                                         |
    413 !    ---------------------------  \       |
    414 !            lalim                 |      |
    415 !    ---------------------------   |      |
    416 !                                  |      |
    417 !    ---------------------------   |      |
    418 !                                  | A    |
    419 !    ---------------------------   |      |
    420 !                                  |      |
    421 !    ---------------------------   |      |
    422 !    lmin                          |      |
    423 !    ----- F_lmin=0 ------------  /      /
    424 !
    425 !    ---------------------------
    426 !   ////////////////////////////
    427 !
    428 !------------------------------------------------------------------------------
    429      
    430 !==============================================================================
    431 ! Calculs initiaux ne faisant pas intervenir les changements de phase
    432 !==============================================================================
    433      
    434 !------------------------------------------------------------------------------
    435 !  1. alim_star est le profil vertical de l'alimentation a la base du
    436 !     panache thermique, calcule a partir de la flotabilite de l'air sec
    437 !  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
    438 !  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
    439 !     panache sec conservatif (e=d=0) alimente selon alim_star
    440 !     Il s'agit d'un calcul de type CAPE
    441 !     zmax_sec est utilise pour determiner la geometrie du thermique.
    442 !------------------------------------------------------------------------------
    443      
    444       entr_star(:,:) = 0.
    445       detr_star(:,:) = 0.
    446       alim_star(:,:) = 0.
    447      
    448       alim_star_tot(:) = 0.
    449      
    450       lmin(:) = 1
    451      
    452 !==============================================================================
    453 ! Calcul du melange et des variables dans le thermique
    454 !==============================================================================
    455      
    456       CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,                &
    457       &    po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,      &
    458       &    lalim,f0,detr_star,entr_star,f_star,ztva,                          &
    459       &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,      &
    460       &    lmin,lev_out,lunout1,igout)
    461      
    462 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
    463 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
    464      
    465 !==============================================================================
     253!===============================================================================
     254! Explicative schemes
     255!===============================================================================
     256
     257!-------------------------------------------------------------------------------
     258! Thermal plume variables
     259!-------------------------------------------------------------------------------
     260     
     261!           top of the model     
     262!     ===========================
     263!                               
     264!     ---------------------------
     265!                                        _
     266!     ----- F_lmax+1=0 ------zmax         \
     267!     lmax                                 |
     268!     ------F_lmax>0-------------          |
     269!                                          |
     270!     ---------------------------          |
     271!                                          |
     272!     ---------------------------          |
     273!                                          |
     274!     ------------------wmax,zmix          |
     275!     lmix                                 |
     276!     ---------------------------          |
     277!                                          |
     278!     ---------------------------          |
     279!                                          | E, D
     280!     ---------------------------          |
     281!                                          |
     282!     --------------------------- rhobarz, f_star, fm, fm0, zw2, fraca
     283!         zt, zu, zv, zo, rho              |
     284!     ---------------------------          |
     285!                                          |
     286!     ---------------------------          |
     287!                                          |
     288!     ---------------------------          |
     289!                                          |
     290!     ------F_lmin+1>0-----------          |
     291!     lmin                                 |
     292!     ----- F_lmin=0 ------------        _/
     293!                               
     294!     ---------------------------
     295!                               
     296!     ===========================
     297!         bottom of the model   
     298     
     299!-------------------------------------------------------------------------------
     300! Zoom on layers k and k-1
     301!-------------------------------------------------------------------------------
     302     
     303!     |     /|\                  |                          |
     304!     |----  |  F_k+1 -----------|--------------------------|   level k+1
     305!     |      |  w_k+1         |                             |
     306!     |                     --|--> D_k                      |
     307!     |                       |                             |   layer k
     308!     |                    <--|--  E_k                      |
     309!     |     /|\               |                             |
     310!     |----  |  F_k ----------|-----------------------------|   level k
     311!     |      |  w_k        |                                |
     312!     |                  --|--> D_k-1                       |
     313!     |                    |                                |   layer k-1
     314!     |                 <--|--  E_k-1                       |
     315!     |     /|\            |                                |
     316!     |----  |  F_k-1 -----|--------------------------------|   level k-1
     317!            |  w_k-1                                       
     318!     0                  fraca                              1
     319!      \__________________/ \______________________________/
     320!          plume (fraca)          environment (1-fraca)   
     321     
     322!===============================================================================
     323! Thermal plumes computation
     324!===============================================================================
     325     
     326!-------------------------------------------------------------------------------
     327! Thermal plumes speeds, fluxes, tracers and temperatures
     328!-------------------------------------------------------------------------------
     329     
     330      CALL thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
     331      &                    zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
     332      &                    detr_star,entr_star,f_star,                        &
     333      &                    ztva,zhla,zqla,zqta,zta,                           &
     334      &                    zw2,zqsa,lmix,lmin)
     335     
     336!-------------------------------------------------------------------------------
    466337! Thermal plumes characteristics: zmax, zmix, wmax
    467 !==============================================================================
    468      
    469       CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,            &
    470       &                     zlev,lmax,zmax,zmix,wmax,f_star,lev_out)
    471      
    472 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    473 ! AB : WARNING: zw2 became its square root in thermcell_height!
    474 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    475      
    476 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
    477 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
    478 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
    479 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
    480      
    481 !==============================================================================
     338!-------------------------------------------------------------------------------
     339     
     340! AB: Careful, zw2 became its square root in thermcell_height!
     341      CALL thermcell_height(ngrid,nlay,lmin,linter,lmix,zw2,                  &
     342      &                     zlev,lmax,zmax,zmix,wmax,f_star)
     343     
     344!===============================================================================
    482345! Closure and mass fluxes computation
    483 !==============================================================================
    484      
    485       CALL thermcell_dry(ngrid,nlay,zlev,pphi,ztv,alim_star,                  &
    486       &                  lalim,lmin,zmax_sec,wmax_sec,lev_out)
    487      
    488 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
    489 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
    490      
    491 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    492 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
    493 ! Apparemment sans importance
    494 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    495       alim_star_clos(:,:) = alim_star(:,:)
    496       alim_star_clos(:,:) = entr_star(:,:) + alim_star(:,:)
    497      
    498 !------------------------------------------------------------------------------
    499 ! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2)
    500 !------------------------------------------------------------------------------
    501      
    502       IF (iflag_thermals_closure.eq.1) THEN
    503          CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                &
    504          &                      lalim,alim_star_clos,f_star,                  &
    505          &                      zmax_sec,wmax_sec,f,lev_out)
    506       ELSEIF (iflag_thermals_closure.eq.2) THEN
    507          CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                &
    508          &                      lalim,alim_star,f_star,                       &
    509          &                      zmax,wmax,f,lev_out)
    510       ELSE
    511          print *, 'ERROR: no closure had been selected!'
    512          call abort
    513       ENDIF
     346!===============================================================================
     347     
     348!-------------------------------------------------------------------------------
     349! Closure
     350!-------------------------------------------------------------------------------
     351     
     352      CALL thermcell_closure(ngrid,nlay,ptimestep,rho,zlev,                   &
     353      &                      lmax,entr_star,zmax,wmax,f)
    514354     
    515355      IF (tau_thermals>1.) THEN
    … …  
    520360      ENDIF
    521361     
    522 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     362!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    523363! Test valable seulement en 1D mais pas genant
    524 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     364!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    525365      IF (.not. (f0(1).ge.0.) ) THEN
    526366         abort_message = '.not. (f0(1).ge.0.)'
    … …  
    529369      ENDIF
    530370     
    531 !------------------------------------------------------------------------------
     371!-------------------------------------------------------------------------------
    532372! Mass fluxes
    533 !------------------------------------------------------------------------------
     373!-------------------------------------------------------------------------------
    534374     
    535375      CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
    536       &                   lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
    537       &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla,               &
    538       &                   lev_out,lunout1,igout)
    539      
    540 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
    541 !      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
    542      
    543 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     376! AB: I remove alimentation, which is in fact included in entr_star. EN COURS
     377!     That is not already done for thermcell_flux.
     378      &                   lalim,lmin,lmax,entr_star,detr_star,                &
     379      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla)
     380     
     381!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    544382! On ne prend pas directement les profils issus des calculs precedents mais on
    545383! s'autorise genereusement une relaxation vers ceci avec une constante de temps
    546 ! tau_thermals (typiquement 1800s).
    547 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     384! tau_thermals (typiquement 1800s sur Terre).
     385!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    548386     
    549387      IF (tau_thermals>1.) THEN
    … …  
    558396      ENDIF
    559397     
    560 !------------------------------------------------------------------------------
     398!-------------------------------------------------------------------------------
    561399! Updraft fraction
    562 !------------------------------------------------------------------------------
     400!-------------------------------------------------------------------------------
    563401     
    564402      DO ig=1,ngrid
    … …  
    577415      ENDDO
    578416     
    579 !==============================================================================
     417!===============================================================================
    580418! Transport vertical
    581 !==============================================================================
    582      
    583 !------------------------------------------------------------------------------
    584 ! Calcul du transport vertical (de qt et tp)
    585 !------------------------------------------------------------------------------
    586      
    587       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
    588       &                 zthl,zdthladj,zta,lmin,lev_out)
    589      
    590       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse,    &
    591       &                 po,pdoadj,zoa,lmin,lev_out)
     419!===============================================================================
     420     
     421!-------------------------------------------------------------------------------
     422! Calcul du transport vertical de la temperature potentielle
     423!-------------------------------------------------------------------------------
     424     
     425      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
     426      &                 zhl,zdthladj,dummy,lmin)
    592427     
    593428      DO l=1,nlay
    594429         DO ig=1,ngrid
    595            pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l) 
     430            pdtadj(ig,l) = zdthladj(ig,l) * zpopsk(ig,l) 
    596431         ENDDO
    597432      ENDDO
    598433     
    599 !------------------------------------------------------------------------------
     434!-------------------------------------------------------------------------------
     435! Calcul du transport vertical des traceurs
     436!-------------------------------------------------------------------------------
     437     
     438      DO iq=1,nq
     439         CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,        &
     440         &                 pq(:,:,iq),pdqadj(:,:,iq),zqa(:,:,iq),lmin)
     441      ENDDO
     442     
     443!-------------------------------------------------------------------------------
    600444! Calcul du transport vertical du moment horizontal
    601 !------------------------------------------------------------------------------
    602      
    603       IF (dvimpl.eq.0) THEN
    604 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    605 ! Calcul du transport de V tenant compte d'echange par gradient
    606 ! de pression horizontal avec l'environnement
    607 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    608          CALL thermcell_dv2(ngrid,nlay,ptimestep,fm0,entr0,masse,fraca,       &
    609          &                  zmax,zu,zv,pduadj,pdvadj,zua,zva,lev_out)
    610       ELSE
    611 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    612 ! Calcul purement conservatif pour le transport de V=(u,v)
    613 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    614          CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
    615          &                 zu,pduadj,zua,lmin,lev_out)
    616          
    617          CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,detr0,masse, &
    618          &                 zv,pdvadj,zva,lmin,lev_out)
    619       ENDIF
    620      
    621 !==============================================================================
    622 ! Calculs de diagnostiques pour les sorties
    623 !==============================================================================
    624      
    625       IF (sorties) THEN
    626          
    627 !------------------------------------------------------------------------------
    628 ! Calcul du niveau de condensation
    629 !------------------------------------------------------------------------------
    630          
    631          DO ig=1,ngrid
    632             nivcon(ig) = 0
    633             zcon(ig) = 0.
    634          ENDDO
    635 !nouveau calcul
    636          do ig=1,ngrid
    637             CHI = zh(ig,1)/(1669.0-122.0*zo(ig,1)/zqsat(ig,1)-zh(ig,1))
    638             pcon(ig) = pplay(ig,1)*(zo(ig,1)/zqsat(ig,1))**CHI
    639          ENDDO
    640          
    641 !IM       do k=1,nlay
    642          do k=1,nlay-1
    643             do ig=1,ngrid
    644                IF ((pcon(ig).le.pplay(ig,k)).and.(pcon(ig).gt.pplay(ig,k+1))) then
    645                   zcon2(ig) = zlay(ig,k) - (pcon(ig) - pplay(ig,k))           &
    646                             / (RG * rho(ig,k)) / 100.
    647                ENDIF
    648             ENDDO
    649          ENDDO
    650          
    651          ierr = 0
    652          
    653          do ig=1,ngrid
    654            IF (pcon(ig).le.pplay(ig,nlay)) then
    655               zcon2(ig) = zlay(ig,nlay) - (pcon(ig) - pplay(ig,nlay))         &
    656                         / (RG * rho(ig,nlay)) / 100.
    657               ierr = 1
    658             ENDIF
    659          ENDDO
    660          
    661          IF (ierr==1) then
    662               write(*,*) 'ERROR: thermal plumes rise the model top'
    663               CALL abort
    664          ENDIF
    665          
    666          IF (prt_level.ge.1) print*,'14b OK convect8'
    667          
    668          do k=nlay,1,-1
    669             do ig=1,ngrid
    670                IF (zqla(ig,k).gt.1e-10) then
    671                   nivcon(ig) = k
    672                   zcon(ig) = zlev(ig,k)
    673                ENDIF
    674             ENDDO
    675          ENDDO
    676          
    677          IF (prt_level.ge.1) print*,'14c OK convect8'
    678          
    679 !------------------------------------------------------------------------------
    680 ! Calcul des moments
    681 !------------------------------------------------------------------------------
    682          
    683          do l=1,nlay
    684             do ig=1,ngrid
    685                q2(ig,l) = 0.
    686                wth2(ig,l) = 0.
    687                wth3(ig,l) = 0.
    688                ratqscth(ig,l) = 0.
    689                ratqsdiff(ig,l) = 0.
    690             ENDDO
    691          ENDDO
    692          
    693          IF (prt_level.ge.1) print*,'14d OK convect8'
    694          
    695          do l=1,nlay
    696            do ig=1,ngrid
    697                zf = fraca(ig,l)
    698                zf2 = zf/(1.-zf)
    699                thetath2(ig,l) = zf2*(ztla(ig,l)-zthl(ig,l))**2
    700                
    701                IF (zw2(ig,l).gt.1.e-10) then
    702                   wth2(ig,l) = zf2*(zw2(ig,l))**2
    703                else
    704                   wth2(ig,l) = 0.
    705                ENDIF
    706                
    707                wth3(ig,l) = zf2*(1-2.*fraca(ig,l))/(1-fraca(ig,l))            &
    708                &                *zw2(ig,l)*zw2(ig,l)*zw2(ig,l)
    709                q2(ig,l) = zf2*(zqta(ig,l)*1000.-po(ig,l)*1000.)**2
    710 !test: on calcul q2/po=ratqsc
    711                ratqscth(ig,l) = sqrt(max(q2(ig,l),1.e-6)/(po(ig,l)*1000.))
    712             ENDDO
    713          ENDDO
    714          
    715 !------------------------------------------------------------------------------
    716 ! Calcul des flux: q, thetal et thetav
    717 !------------------------------------------------------------------------------
    718          
    719          do l=1,nlay
    720             do ig=1,ngrid
    721                wq(ig,l)=fraca(ig,l)*zw2(ig,l)*(zqta(ig,l)*1000.-po(ig,l)*1000.)
    722                wthl(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztla(ig,l)-zthl(ig,l))
    723                wthv(ig,l)=fraca(ig,l)*zw2(ig,l)*(ztva(ig,l)-ztv(ig,l))
    724             ENDDO
    725          ENDDO
    726          
    727          CALL thermcell_alp(ngrid,nlay,ptimestep,                             &
    728          &                  pplay,pplev,                                      &
    729          &                  fm0,entr0,lmax,                                   &
    730          &                  Ale_bl,Alp_bl,lalim_conv,wght_th,                 &
    731          &                  zw2,fraca,pcon,                                   &
    732          &                  rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,    &
    733          &                  pbl_tke,pctsrf,omega,airephy,                     &
    734          &                  zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,&
    735          &                  n2,s2,ale_bl_stat,                                &
    736          &                  therm_tke_max,env_tke_max,                        &
    737          &                  alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,       &
    738          &                  alp_bl_conv,alp_bl_stat)
    739          
    740 !------------------------------------------------------------------------------
    741 ! Calcul du ratqscdiff
    742 !------------------------------------------------------------------------------
    743      
    744          var = 0.
    745          vardiff = 0.
    746          ratqsdiff(:,:) = 0.
    747          
    748          DO l=1,nlay
    749             DO ig=1,ngrid
    750                IF (l<=lalim(ig)) THEN
    751                   var = var + alim_star(ig,l) * zqta(ig,l) * 1000.
    752                ENDIF
    753             ENDDO
    754          ENDDO
    755          
    756          IF (prt_level.ge.1) print*,'14f OK convect8'
    757      
    758          DO l=1,nlay
    759             DO ig=1,ngrid
    760                IF (l<=lalim(ig)) THEN
    761                   zf  = fraca(ig,l)
    762                   zf2 = zf / (1.-zf)
    763                   vardiff = vardiff + alim_star(ig,l)                         &
    764                           * (zqta(ig,l) * 1000. - var)**2
    765                ENDIF
    766             ENDDO
    767          ENDDO
    768      
    769          IF (prt_level.ge.1) print*,'14g OK convect8'
    770          
    771          DO l=1,nlay
    772             DO ig=1,ngrid
    773                ratqsdiff(ig,l) = sqrt(vardiff) / (po(ig,l) * 1000.)   
    774             ENDDO
    775          ENDDO
    776          
    777       ENDIF
     445!-------------------------------------------------------------------------------
     446     
     447      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
     448      &                 zu,pduadj,zua,lmin)
     449     
     450      CALL thermcell_dq(ngrid,nlay,ptimestep,fm0,entr0,detr0,masse,           &
     451      &                 zv,pdvadj,zva,lmin)
    778452     
    779453     
    780454RETURN
    781455END
    782 
    783 
    784 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    785 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    786 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    787 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    788 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    789 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    790 
    791 
    792 SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,              &
    793                        ztva,zqla,f_star,zw2,comment)
    794      
    795      
    796       USE print_control_mod, ONLY: prt_level
    797      
    798       IMPLICIT NONE
    799      
    800      
    801 !==============================================================================
    802 ! Declaration
    803 !==============================================================================
    804      
    805 !      inputs:
    806 !      -------
    807      
    808       INTEGER ngrid
    809       INTEGER nlay
    810       INTEGER long(ngrid)
    811      
    812       REAL pplev(ngrid,nlay+1)
    813       REAL pplay(ngrid,nlay)
    814       REAL ztv(ngrid,nlay)
    815       REAL po(ngrid,nlay)
    816       REAL ztva(ngrid,nlay)
    817       REAL zqla(ngrid,nlay)
    818       REAL f_star(ngrid,nlay)
    819       REAL zw2(ngrid,nlay)
    820       REAL seuil
    821      
    822       CHARACTER*21 comment
    823      
    824 !      local:
    825 !      ------
    826      
    827       INTEGER i, k
    828      
    829 !==============================================================================
    830 ! Test
    831 !==============================================================================
    832      
    833       IF (prt_level.ge.1) THEN
    834          write(*,*) 'WARNING: in test, ', comment
    835       ENDIF
    836            
    837       return
    838      
    839 !  test sur la hauteur des thermiques ...
    840       do i=1,ngrid
    841 !IMtemp IF (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
    842         IF (prt_level.ge.10) then
    843             print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
    844             print *, '  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
    845             do k=1,nlay
    846                write(6,'(i3,7f10.3)') k,pplay(i,k),ztv(i,k),1000*po(i,k),ztva(i,k),1000*zqla(i,k),f_star(i,k),zw2(i,k)
    847             ENDDO
    848         ENDIF
    849       ENDDO
    850      
    851      
    852 RETURN
    853 END
    854 
    855 
    856 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    857 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    858 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    859 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    860 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    861 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    862 
    863 
    864 SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,            &
    865                                    rg,pplev,therm_tke_max)
    866      
    867      
    868 !==============================================================================
    869 !   Calcul du transport verticale dans la couche limite en presence
    870 !   de "thermiques" explicitement representes
    871 !   calcul du dq/dt une fois qu'on connait les ascendances
    872 !
    873 !   Transport de la TKE par le thermique moyen pour la fermeture en ALP
    874 !   On transporte pbl_tke pour donner therm_tke
    875 !==============================================================================
    876      
    877       USE print_control_mod, ONLY: prt_level
    878      
    879       IMPLICIT NONE
    880      
    881      
    882 !==============================================================================
    883 ! Declarations
    884 !==============================================================================
    885      
    886       integer ngrid,nlay,nsrf
    887      
    888       real ptimestep
    889       real masse0(ngrid,nlay),fm0(ngrid,nlay+1),pplev(ngrid,nlay+1)
    890       real entr0(ngrid,nlay),rg
    891       real therm_tke_max(ngrid,nlay)
    892       real detr0(ngrid,nlay)
    893      
    894       real masse(ngrid,nlay),fm(ngrid,nlay+1)
    895       real entr(ngrid,nlay)
    896       real q(ngrid,nlay)
    897      
    898       real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
    899      
    900       real zzm
    901      
    902       integer ig,k
    903       integer isrf
    904      
    905 !------------------------------------------------------------------------------
    906 ! Calcul du detrainement
    907 !------------------------------------------------------------------------------
    908      
    909       do k=1,nlay
    910          detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
    911          masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
    912       ENDDO
    913      
    914 !------------------------------------------------------------------------------
    915 ! Decalage vertical des entrainements et detrainements.
    916 !------------------------------------------------------------------------------
    917      
    918       masse(:,1)=0.5*masse0(:,1)
    919       entr(:,1)=0.5*entr0(:,1)
    920       detr(:,1)=0.5*detr0(:,1)
    921       fm(:,1)=0.
    922      
    923       do k=1,nlay-1
    924          masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
    925          entr(:,k+1)=0.5*(entr0(:,k)+entr0(:,k+1))
    926          detr(:,k+1)=0.5*(detr0(:,k)+detr0(:,k+1))
    927          fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
    928       ENDDO
    929      
    930       fm(:,nlay+1)=0.
    931      
    932 !!! nrlmd le 16/09/2010
    933 !   calcul de la valeur dans les ascendances
    934 !      do ig=1,ngrid
    935 !         qa(ig,1) = q(ig,1)
    936 !      ENDDO
    937      
    938       q(:,:)=therm_tke_max(:,:)
    939      
    940       do ig=1,ngrid
    941          qa(ig,1)=q(ig,1)
    942       ENDDO
    943      
    944       IF (1==1) then
    945          do k=2,nlay
    946             do ig=1,ngrid
    947                IF ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
    948                   qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))         &
    949                   &        / (fm(ig,k+1)+detr(ig,k))
    950                else
    951                   qa(ig,k)=q(ig,k)
    952                ENDIF
    953                
    954 !               IF (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
    955 !               IF (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
    956             ENDDO
    957          ENDDO
    958          
    959 !------------------------------------------------------------------------------
    960 ! Calcul du flux subsident
    961 !------------------------------------------------------------------------------
    962          
    963          do k=2,nlay
    964             do ig=1,ngrid
    965                wqd(ig,k) = fm(ig,k) * q(ig,k)
    966                IF (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
    967             ENDDO
    968          ENDDO
    969          
    970          do ig=1,ngrid
    971             wqd(ig,1) = 0.
    972             wqd(ig,nlay+1) = 0.
    973          ENDDO
    974          
    975 !------------------------------------------------------------------------------
    976 ! Calcul des tendances
    977 !------------------------------------------------------------------------------
    978          
    979          do k=1,nlay
    980             do ig=1,ngrid
    981                q(ig,k) = q(ig,k) + ptimestep / masse(ig,k)                    &
    982                &       * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k)        &
    983                &       - wqd(ig,k) + wqd(ig,k+1))
    984             ENDDO
    985          ENDDO
    986       ENDIF
    987      
    988       therm_tke_max(:,:) = q(:,:)
    989      
    990      
    991 RETURN
    992 END
    993 

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_mod.F90

    r2092 r2127  
    66! Flags for computations
    77                                                                  !  default
    8       INTEGER,PARAMETER :: iflag_thermals_optflux     = 1         !  0        !
    9       INTEGER,PARAMETER :: iflag_thermals_closure     = 2         !  2        !
    10       INTEGER,PARAMETER :: iflag_thermals             = 18        !  18       !
    11      
    12 ! Flags for (terrestrial) diagnoses
    13      
    14       LOGICAL,PARAMETER :: sorties                    = .false.   !  false
    15       INTEGER,PARAMETER :: iflag_trig_bl              = 1         !  1
    16       INTEGER,PARAMETER :: iflag_clos_bl              = 1         !  1
    17       INTEGER,PARAMETER :: iflag_coupl                = 5         !  5
     8      INTEGER,PARAMETER :: iflag_thermals_optflux     = 1         !  0        -
     9      INTEGER,PARAMETER :: dqimpl                     = 1         !  1        flag for thermcell_dq version (1 : implicit scheme || 0 : explicit scheme)
    1810     
    1911! Physical parameters
    2012     
    21       REAL,PARAMETER :: fact_thermals_ed_dz           = 0.007     !  0.007    !
    2213      REAL,PARAMETER :: r_aspect_thermals             = 1.0       !           Aspect ratio of the thermals (width / height)
    2314      REAL,PARAMETER :: tau_thermals                  = 0.        !  0.       Relaxation time
    24       REAL,PARAMETER :: betalpha                      = 0.9       !  0.9      !
    25       REAL,PARAMETER :: afact                         = 2./3.     !  2./3.    !
    26       REAL,PARAMETER :: fact_epsilon                  = 0.000     !  0.002    !
    27       REAL,PARAMETER :: detr_q_power                  = 0.5       !  0.5      !
    28       REAL,PARAMETER :: detr_q_coef                   = 0.012     !  0.012    !
    29       REAL,PARAMETER :: mix0                          = 0.        !  0.       !
    30       REAL,PARAMETER :: detr_min                      = 1.d-5     !  1.e-5    Minimal detrainment value
    31       REAL,PARAMETER :: entr_min                      = 1.d-5     !  1.e-5    Maximal detrainment value
    32       REAL,PARAMETER :: alphamax                      = 0.7       !           Maximal permitted updraft fraction
    33       REAL,PARAMETER :: fomass_max                    = 0.5       !           Maximal permitted outgoing layer mass fraction
    34       REAL,PARAMETER :: pres_limit                    = 1.e5      !  1.e5     !
     15      REAL,PARAMETER :: betalpha                      = 1.0       !  0.9      - between 0 (e=d) and 1 (rho*fraca=cst)
     16      REAL,PARAMETER :: afact                         = 1.        !  2./3.    - buoyancy contribution, between 0 and 1
     17      REAL,PARAMETER :: fact_epsilon                  = 1.e-4     !  2.e-3    - friction
     18      REAL,PARAMETER :: nu                            = 0.000     !           Geometrical contributions to entrainment and detrainment
     19      REAL,PARAMETER :: alphamax                      = 0.7       !  0.7      Maximal permitted updraft fraction
     20      REAL,PARAMETER :: fomass_max                    = 0.5       !  0.5      Maximal permitted outgoing layer mass fraction
     21      REAL,PARAMETER :: pres_limit                    = 2.e5      !  1.e5     -
    3522     
    3623!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    4330!      ngrid.
    4431!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    45       INTEGER,PARAMETER :: linf                       = 2         !     1
     32      INTEGER,PARAMETER :: linf                       = 1         !     1
    4633     
    4734!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    5037!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    5138      REAL,PARAMETER :: d_temp                        = 0.        !     0.
    52      
    53 ! Parameters for diagnoses
    54      
    55       REAL,PARAMETER :: alp_bl_k                      = 0.5       !     0.5
    5639     
    5740! Physical constants
    … …  
    6447     
    6548!$OMP THREADPRIVATE(RTT, RG, RKAPPA, RPI, RD)
    66      
    67 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    68 ! AB : Parameters needed only for a loop in thermcell_alp (diagnoses).
    69 !      Maybe to be removed.
    70 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    71       INTEGER,PARAMETER :: nbsrf                      = 1
    7249     
    7350     

  • trunk/LMDZ.GENERIC/libf/phystd/thermcell_plume.F90

    r2113 r2127  
    22!
    33!
    4 SUBROUTINE thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,                     &
    5                            zthl,po,zl,rhobarz,zlev,pplev,pphi,zpspsk,         &
    6                            alim_star,alim_star_tot,lalim,f0,detr_star,        &
    7                            entr_star,f_star,ztva,ztla,zqla,zqta,zha,          &
    8                            zw2,w_est,ztva_est,zqsatth,lmix,lmix_bis,          &
    9                            lmin,lev_out,lunout1,igout)
    10      
    11      
    12 !==============================================================================
    13 ! thermcell_plume: calcule les valeurs de qt, thetal et w dans l ascendance
    14 ! AB : ql means "liquid water mass mixing ratio"
    15 !      qt means "total water mass mixing ratio"
    16 !      TP means "potential temperature"
    17 !      TRPV means "virtual potential temperature with latent heat release" 
    18 !      TPV means "virtual potential temperature"
    19 !      TR means "temperature with latent heat release"
    20 !==============================================================================
     4SUBROUTINE thermcell_plume(ngrid,nlay,nq,ptimestep,ztv,                       &
     5                           zhl,zqt,zql,rhobarz,zlev,pplev,pphi,zpopsk,        &
     6                           detr_star,entr_star,f_star,                        &
     7                           ztva,zhla,zqla,zqta,zta,                           &
     8                           zw2,zqsa,lmix,lmin)
     9     
     10     
     11!===============================================================================
     12!  Purpose: calcule les valeurs de qt, thetal et w dans l ascendance
     13! 
     14!  Nota Bene
     15!     ql   means "non-gaseous water mass mixing ratio" (liquid and solid)
     16!     qv   means "vapor mass mixing ratio"
     17!     qt   means "total water mass mixing ratio"
     18!     TP   means "potential temperature"
     19!     TRPV means "virtual potential temperature with latent heat release" 
     20!     TPV  means "virtual potential temperature"
     21!     TR   means "temperature with latent heat release"
     22!===============================================================================
    2123     
    2224      USE print_control_mod, ONLY: prt_level
    2325      USE watercommon_h, ONLY: RLvCp, RETV, Psat_water
     26      USE tracer_h, ONLY: igcm_h2o_vap
    2427      USE thermcell_mod
    2528     
    … …  
    2730     
    2831     
    29 !==============================================================================
     32!===============================================================================
    3033! Declaration
    31 !==============================================================================
     34!===============================================================================
    3235     
    3336!     Inputs:
    3437!     -------
    3538     
    36       INTEGER itap
    37       INTEGER ngrid
    38       INTEGER nlay
    39       INTEGER lunout1
    40       INTEGER igout
    41       INTEGER lev_out                           ! niveau pour les print
    42      
    43       REAL ptimestep                            ! time step
     39      INTEGER ngrid, nlay, nq
     40     
     41      REAL ptimestep
     42      REAL rhobarz(ngrid,nlay)                  ! Levels density
     43      REAL zlev(ngrid,nlay+1)                   ! Levels altitude
     44      REAL pplev(ngrid,nlay+1)                  ! Levels pressure
     45      REAL pphi(ngrid,nlay)                     ! Geopotential
     46      REAL zpopsk(ngrid,nlay)                   ! Exner function
     47     
    4448      REAL ztv(ngrid,nlay)                      ! TRPV environment
    45       REAL zthl(ngrid,nlay)                     ! TP   environment
    46       REAL po(ngrid,nlay)                       ! qt   environment
    47       REAL zl(ngrid,nlay)                       ! ql   environment
    48       REAL rhobarz(ngrid,nlay)                  ! levels density
    49       REAL zlev(ngrid,nlay+1)                   ! levels altitude
    50       REAL pplev(ngrid,nlay+1)                  ! levels pressure
    51       REAL pphi(ngrid,nlay)                     ! geopotential
    52       REAL zpspsk(ngrid,nlay)                   ! Exner function
     49      REAL zhl(ngrid,nlay)                      ! TP   environment
     50      REAL zqt(ngrid,nlay)                      ! qt   environment
     51      REAL zql(ngrid,nlay)                      ! ql   environment
    5352     
    5453!     Outputs:
    … …  
    5655     
    5756      INTEGER lmin(ngrid)                       ! plume base level (first unstable level)
    58       INTEGER lalim(ngrid)                      ! higher alimentation level
    5957      INTEGER lmix(ngrid)                       ! maximum vertical speed level
    60       INTEGER lmix_bis(ngrid)                   ! maximum vertical speed level (modified)
    61      
    62       REAL alim_star(ngrid,nlay)                ! normalized alimentation
    63       REAL alim_star_tot(ngrid)                 ! integrated alimentation
    64      
    65       REAL f0(ngrid)                            ! previous time step mass flux norm
    6658     
    6759      REAL detr_star(ngrid,nlay)                ! normalized detrainment
    … …  
    7062     
    7163      REAL ztva(ngrid,nlay)                     ! TRPV plume (after mixing)
    72       REAL ztva_est(ngrid,nlay)                 ! TRPV plume (before mixing)
    73       REAL ztla(ngrid,nlay)                     ! TP   plume
     64      REAL zhla(ngrid,nlay)                     ! TP   plume ?
    7465      REAL zqla(ngrid,nlay)                     ! ql   plume (after mixing)
    75       REAL zqta(ngrid,nlay)                     ! qt   plume
    76       REAL zha(ngrid,nlay)                      ! TRPV plume
    77      
    78       REAL w_est(ngrid,nlay+1)                  ! updraft square vertical speed (before mixing)
    79       REAL zw2(ngrid,nlay+1)                    ! updraft square vertical speed (after mixing)
    80      
    81       REAL zqsatth(ngrid,nlay)                  ! saturation vapor pressure (after mixing)
     66      REAL zqta(ngrid,nlay)                     ! qt   plume ?
     67      REAL zqsa(ngrid,nlay)                     ! qsat plume (after mixing)
     68      REAL zw2(ngrid,nlay+1)                    ! w    plume (after mixing)
    8269     
    8370!     Local:
    … …  
    8572     
    8673      INTEGER ig, l, k
    87       INTEGER lt
    88       INTEGER lm
    89      
     74     
     75      REAL ztva_est(ngrid,nlay)                 ! TRPV plume (before mixing)
    9076      REAL zqla_est(ngrid,nlay)                 ! ql   plume (before mixing)
    9177      REAL zta_est(ngrid,nlay)                  ! TR   plume (before mixing)
    92       REAL zbuoy(ngrid,nlay)                    ! plume buoyancy
    93       REAL zbuoyjam(ngrid,nlay)                 ! plume buoyancy (modified)
    94      
    95       REAL ztemp(ngrid)                         ! temperature for saturation vapor pressure computation in plume
    96       REAL zqsat(ngrid)                         ! saturation vapor pressure (before mixing)
    97       REAL zdz                                  ! layers height
    98       REAL ztv2(ngrid,nlay)                     ! ztv + d_temp * Dirac(l=lmin)
    99      
    100       REAL zalpha                               !
    101       REAL zdqt(ngrid,nlay)                     !
     78      REAL zqsa_est(ngrid)                      ! qsat plume (before mixing)
     79      REAL zw2_est(ngrid,nlay+1)                ! w    plume (before mixing)
     80     
     81      REAL zta(ngrid,nlay)                      ! TR   plume (after mixing)
     82     
     83      REAL zbuoy(ngrid,nlay)                    ! Plume buoyancy
     84      REAL ztemp(ngrid)                         ! Temperature for saturation vapor pressure computation in plume
     85      REAL zdz                                  ! Layers heights
     86      REAL ztv2(ngrid,nlay)                     ! ztv + d_temp * Dirac(l=linf)
     87     
    10288      REAL zbetalpha                            !
    10389      REAL zdw2                                 !
    10490      REAL zdw2bis                              !
    10591      REAL zw2fact                              !
    106       REAL zw2factbis                           !
    107       REAL zw2m                                 !
    108       REAL zdzbis                               !
    109       REAL coefzlmel                            !
    110       REAL zdz2                                 !
    111       REAL zdz3                                 !
    112       REAL lmel                                 !
    113       REAL zlmel                                !
    114       REAL zlmelup                              !
    115       REAL zlmeldwn                             !
    116       REAL zlt                                  !
    117       REAL zltdwn                               !
    118       REAL zltup                                ! useless here
    119      
    120       REAL psat                                 ! dummy argument for Psat_water()
    121      
    122       LOGICAL active(ngrid)                     ! if the plume is active at ig,l (speed and incoming mass flux > 0 or l=lmin)
    123       LOGICAL activetmp(ngrid)                  ! if the plus is active at ig,l (active=true and outgoing mass flux > 0)
    124       LOGICAL, SAVE :: first = .true.           ! if it is the first time step
    125      
    126 !$OMP THREADPRIVATE(first)
    127      
    128 !==============================================================================
     92      REAL zw2m                                 ! Average vertical velocity between two successive levels
     93      REAL gamma                                ! Plume acceleration term (to compute vertical velocity)
     94      REAL test                                 ! Test to know how to compute entrainment and detrainment
     95       
     96      REAL psat                                 ! Dummy argument for Psat_water()
     97     
     98      LOGICAL active(ngrid)                     ! If the plume is active at ig (speed and incoming mass flux > 0 or l=lmin)
     99      LOGICAL activetmp(ngrid)                  ! If the plume is active at ig (active=true and outgoing mass flux > 0)
     100     
     101!===============================================================================
    129102! Initialization
    130 !==============================================================================
     103!===============================================================================
    131104     
    132105      zbetalpha = betalpha / (1. + betalpha)
    133106     
    134       ztva(:,:)        = ztv(:,:)               ! ztva     is set to the virtual potential temperature without latent heat release
    135       ztva_est(:,:)    = ztva(:,:)              ! ztva_est is set to the virtual potential temperature without latent heat release
    136       ztla(:,:)        = zthl(:,:)              ! ztla     is set to the potential temperature
    137       zqta(:,:)        = po(:,:)                ! zqta     is set to qt
    138       zqla(:,:)        = 0.                     ! zqla     is set to ql
    139       zqla_est(:,:)    = 0.                     ! zqla_est is set to ql
    140       zha(:,:)         = ztva(:,:)              ! zha      is set to the plume virtual potential temperature without latent heat release
    141      
    142       zqsat(:)         = 0.
    143       zqsatth(:,:)     = 0.
    144      
    145       w_est(:,:)       = 0.
     107      ztva(:,:)        = ztv(:,:)                                                ! ztva     is set to TPV environment
     108      zhla(:,:)        = zhl(:,:)                                                ! zhla     is set to TP  environment
     109      zqta(:,:)        = zqt(:,:)                                                ! zqta     is set to qt  environment
     110      zqla(:,:)        = zql(:,:)                                                ! zqla     is set to ql  environment
     111     
     112      zqsa_est(:)      = 0.
     113      zqsa(:,:)        = 0.
     114     
     115      zw2_est(:,:)     = 0.
    146116      zw2(:,:)         = 0.
    147117     
    148118      zbuoy(:,:)       = 0.
    149       zbuoyjam(:,:)    = 0.
    150119     
    151120      f_star(:,:)      = 0.
    152121      detr_star(:,:)   = 0.
    153122      entr_star(:,:)   = 0.
    154       alim_star(:,:)   = 0.
    155       alim_star_tot(:) = 0.
    156123     
    157124      lmix(:)          = 1
    158       lmix_bis(:)      = 2
    159       lalim(:)         = 1
    160       lmin(:)          = linf
     125      lmin(:)          = 1
    161126     
    162127      ztv2(:,:)        = ztv(:,:)
    163128      ztv2(:,linf)     = ztv(:,linf) + d_temp
    164129     
    165 !==============================================================================
    166 ! 0. Calcul de l'alimentation
    167 !==============================================================================
    168      
    169 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    170 ! AB : Convective plumes can go off from every layer above the linf-th and
    171 !      where pressure is lesser than pres_limit (cf. thermcell_plume).
    172 !      The second constraint is added to avoid the parametrization occurs too
    173 !      high when the low atmosphere is stable.
    174 !      However, once there is a triggered plume, it can rise as high as its
    175 !      velocity allows it (it can overshoot).
    176 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     130!===============================================================================
     131! First layer computation
     132!===============================================================================
     133     
    177134      DO ig=1,ngrid
    178135         active(ig) = .false.
    179136         l = linf
    180          DO WHILE ((.not.active(ig)) .and. pplev(ig,l+1).gt.pres_limit .and. l.lt.nlay)
    181             IF (ztv2(ig,l).gt.ztv2(ig,l+1)) THEN
     137         DO WHILE (.not.active(ig).and.(pplev(ig,l+1).GT.pres_limit).and.(l.LT.nlay))
     138            zdz = (zlev(ig,l+1) - zlev(ig,l))
     139            zw2(ig,l+1) = 2. * afact * RG * zdz                               &  ! Do we have to divide by 1+betalpha (consider entrainment) ?
     140            &           * (ztv2(ig,l) - ztv2(ig,l+1))                         &
     141            &           / (ztv2(ig,l+1) * (1. + betalpha))
     142            zw2m = zw2(ig,l+1) / 2.
     143            entr_star(ig,l) = zdz * zbetalpha * (afact * RG * (ztv2(ig,l)     &
     144            &               - ztv2(ig,l+1)) / ztv2(ig,l+1) / zw2m - fact_epsilon)
     145            IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(entr_star(ig,l).GT.0.)) THEN
    182146               active(ig) = .true.
    183147               lmin(ig) = l
     148               f_star(ig,l+1) = entr_star(ig,l)
     149               zw2_est(ig,l+1) = zw2(ig,l+1)
     150            ELSE
     151               zw2(ig,l+1) = 0.
     152               entr_star(ig,l) = 0.
    184153            ENDIF
    185154            l = l + 1
    … …  
    187156      ENDDO
    188157     
    189 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    190 ! AB : On pourrait n'appeler thermcell_alim que si la plume est active
    191 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    192       CALL thermcell_alim(ngrid,nlay,ztv2,zlev,alim_star,lalim,lmin)
    193      
    194 !==============================================================================
    195 ! 1. Calcul dans la premiere couche
    196 !==============================================================================
    197      
    198       DO ig=1,ngrid
    199          IF (active(ig)) THEN
    200            
    201 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    202 ! AB : plume takes the environment features for every layer below lmin.
    203 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    204             DO l=1,lmin(ig)
    205                ztla(ig,l) = zthl(ig,l)
    206                zqta(ig,l) = po(ig,l)
    207                zqla(ig,l) = zl(ig,l)
    208             ENDDO
    209            
    210             l = lmin(ig)
    211             f_star(ig,l+1) = alim_star(ig,l)
    212            
    213             zw2(ig,l+1) = 2. * RG * (zlev(ig,l+1) - zlev(ig,l))               &
    214             &                * (ztv2(ig,l) - ztv2(ig,l+1)) / ztv2(ig,l+1)
    215            
    216             w_est(ig,l+1) = zw2(ig,l+1)
    217          ENDIF
    218       ENDDO
    219      
    220 !==============================================================================
    221 ! 2. Boucle de calcul de la vitesse verticale dans le thermique
    222 !==============================================================================
     158!===============================================================================
     159! Thermal plumes computations
     160!===============================================================================
    223161     
    224162      DO l=2,nlay-1
    225163         
    226 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    227 ! AB : we decide here if the plume is still active or not. When the plume's
    228 !      first level is reached, we set active to "true". Otherwise, it is given
    229 !      by zw2, f_star, alim_star and entr_star.
    230 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     164!-------------------------------------------------------------------------------
     165! Check if thermal plume is (still) active
     166!-------------------------------------------------------------------------------
     167         
     168!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     169! AB: we decide here if the plume is still active or not. When the plume's
     170!     first level is reached, we set active to "true". Otherwise, it is given
     171!     by zw2, f_star and entr_star.
     172!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    231173         DO ig=1,ngrid
    232174            IF (l==lmin(ig)+1) THEN
    … …  
    235177           
    236178            active(ig) = active(ig)                                           &
    237             &            .and. zw2(ig,l)>1.e-10                               &
    238             &            .and. f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)>1.e-10
    239          ENDDO
    240          
    241          ztemp(:) = zpspsk(:,l) * ztla(:,l-1)
    242          
    243          DO ig=1,ngrid
    244             CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsat(ig))
    245          ENDDO
    246          
    247 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    248 ! AB : we compute thermodynamical values and speed in the plume in the layer l
    249 !      without mixing with environment.
    250 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     179            &          .and. (zw2(ig,l).GT.1.e-10)                            &
     180            &          .and. (f_star(ig,l) + entr_star(ig,l)).GT.1.e-10
     181         ENDDO
     182         
     183         ztemp(:) = zpopsk(:,l) * zhla(:,l-1)
     184         
     185         DO ig=1,ngrid
     186            CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa_est(ig))
     187         ENDDO
     188         
     189!-------------------------------------------------------------------------------
     190! Vertical speed (before mixing between plume and environment)
     191!-------------------------------------------------------------------------------
    251192         
    252193         DO ig=1,ngrid
    253194            IF (active(ig)) THEN
    254                zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsat(ig))                   ! zqla_est set to ql   plume
    255                ztva_est(ig,l) = ztla(ig,l-1)*zpspsk(ig,l)+RLvCp*zqla_est(ig,l)   ! ztva_est set to TR   plume
    256                zta_est(ig,l)  = ztva_est(ig,l)                                   ! zta_est  set to TR   plume
    257                ztva_est(ig,l) = ztva_est(ig,l)/zpspsk(ig,l)                      ! ztva_est set to TRP  plume
    258                ztva_est(ig,l) = ztva_est(ig,l)*(1.+RETV*(zqta(ig,l-1)         &  ! ztva_est set to TRPV plume
    259                &              - zqla_est(ig,l))-zqla_est(ig,l))
    260                
    261                zbuoy(ig,l) = RG * (ztva_est(ig,l)-ztv(ig,l)) / ztv(ig,l)
     195               zqla_est(ig,l) = max(0.,zqta(ig,l-1)-zqsa_est(ig))                ! zqla_est set to ql   plume
     196               zta_est(ig,l)  = zhla(ig,l-1) * zpopsk(ig,l)                   &  ! zta_est  set to TR   plume
     197               &              + RLvCp * zqla_est(ig,l)
     198               ztva_est(ig,l) = zta_est(ig,l) / zpopsk(ig,l)                  &  ! ztva_est set to TRPV plume
     199               &              * (1. + RETV * (zqta(ig,l-1)-zqla_est(ig,l)) - zqla_est(ig,l))
     200               
     201               zbuoy(ig,l) = RG * (ztva_est(ig,l) - ztv(ig,l)) / ztv(ig,l)
    262202               zdz = zlev(ig,l+1) - zlev(ig,l)
    263203               
    264 !==============================================================================
    265 ! 3. Calcul de la flotabilite modifiee par melange avec l'air au dessus
    266 !==============================================================================
    267                
    268                lmel = fact_thermals_ed_dz * zlev(ig,l)
    269                zlmel = zlev(ig,l) + lmel
    270                lt = l + 1
    271                zlt = zlev(ig,lt)
    272                zdz2 = zlev(ig,lt) - zlev(ig,l)
    273                
    274                DO while (lmel.gt.zdz2)
    275                   lt   = lt + 1
    276                   zlt  = zlev(ig,lt)
    277                   zdz2 = zlev(ig,lt) - zlev(ig,l)
    278                ENDDO
    279                
    280 !               IF (lt-l.gt.1) THEN
    281 !                  print *, 'WARNING: lt is greater than l+1!'
    282 !                  print *, 'l,lt', l, lt
    283 !               ENDIF
    284                
    285                zdz3 = zlev(ig,lt+1) - zlt
    286                zltdwn = zlev(ig,lt) - zdz3 / 2
    287                zlmelup = zlmel + (zdz / 2)
    288                coefzlmel = Min(1.,(zlmelup - zltdwn) / zdz)
    289                
    290                zbuoyjam(ig,l) = 1.* RG * (coefzlmel *                         &
    291                &                (ztva_est(ig,l) - ztv(ig,lt)) / ztv(ig,lt)    &
    292                &              + (1. - coefzlmel) *                            &
    293                &                (ztva_est(ig,l) - ztv(ig,lt-1)) / ztv(ig,lt-1)) &
    294                &              + 0. * zbuoy(ig,l)
    295                
    296 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    297 ! AB : initial formulae
    298 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     204!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     205! AB: initial formulae
     206!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    299207!               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
    300208!               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
    301209!               zdw2bis = afact * zbuoy(ig,l-1) / fact_epsilon
    302 !               w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2)+zdw2)
    303 !               w_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(w_est(ig,l)-zdw2bis)+zdw2)
    304 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    305 ! AB : own derivation for w_est (Rio 2010 formula with e=d=0)
    306 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     210!               zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2)+zdw2)
     211!               zw2_est(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2_est(ig,l)-zdw2bis)+zdw2)
     212!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     213! AB: own derivation for zw2_est (Rio et al. 2010)
     214!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    307215!               zw2fact = 2. * fact_epsilon * zdz
    308216!               zdw2 = 2. * afact * zbuoy(ig,l) * zdz
    309217               zw2fact = 2. * fact_epsilon * zdz / (1. + betalpha)
    310218               zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
    311                w_est(ig,l+1) = Max(0., exp(-zw2fact) * w_est(ig,l) + zdw2)
    312                
    313 !               IF (w_est(ig,l+1).le.0.) THEN
    314 !                  print *, 'WARNING: w_est is negative!'
    315 !                  print *, 'l,w_est', l+1, w_est(ig,l+1)
    316 !               ENDIF
    317             ENDIF
    318          ENDDO
    319          
    320 !==============================================================================
    321 ! 4. Calcul de l'entrainement et du detrainement
    322 !==============================================================================
     219               zw2_est(ig,l+1) = Max(0., exp(-zw2fact) * zw2_est(ig,l) + zdw2)
     220            ENDIF
     221         ENDDO
     222         
     223!-------------------------------------------------------------------------------
     224! Mass flux, entrainment and detrainment
     225!-------------------------------------------------------------------------------
    323226         
    324227         DO ig=1,ngrid
    … …  
    326229               
    327230               zdz = zlev(ig,l+1) - zlev(ig,l)
    328                
    329 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    330 ! AB : The next test is added to avoid divisions by zero when w_est vanishes.
    331 !      Indeed, entr and detr computed here are of no importance because w_est
    332 !      <= 0 means it will be the last layer reached by the plume and then they
    333 !      will be reset in thermcell_flux.
    334 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    335                IF (w_est(ig,l+1).eq.0.) THEN
    336                   zw2m = 1.
    337                   zalpha = 0.
     231               zw2m = (zw2_est(ig,l+1) + zw2_est(ig,l)) / 2.                     ! AB: est-ce la bonne vitesse a utiliser ?
     232               gamma = afact * zbuoy(ig,l) - fact_epsilon * zw2m
     233               
     234               IF (zw2_est(ig,l).GT.0.) THEN
     235                  test = gamma / zw2m - nu
    338236               ELSE
    339                   zw2m = w_est(ig,l+1)
    340                   zalpha = f0(ig) * f_star(ig,l) / sqrt(w_est(ig,l+1)) / rhobarz(ig,l)
     237                  print *, 'ERROR: zw2_est is negative while plume is active!'
     238                  print *, 'ig,l', ig, l
     239                  print *, 'zw2_est', zw2_est(ig,l)
     240                  call abort
    341241               ENDIF
    342242               
    343 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    344 ! AB : The next test is added to avoid a division by zero if there is no water
    345 !      in the environment.
    346 !      In the case where there is no water in the env. but water in the plume
    347 !      (ascending from depth) we set the effect on detrainment equal to zero
    348 !      but at the next time step, po will be positive thanks to the mixing and
    349 !      then the physical effect of the water gradient will be taken on board.
    350 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    351                IF (po(ig,l).lt.1.e-6) THEN
    352 !                  print *, 'WARNING: po=0 in layer',l,'!'
    353 !                  print *, 'po,zqta', po(ig,l), zqta(ig,l-1)
    354                   zdqt(ig,l) = 0.0
     243               IF (test.gt.0.) THEN
     244                  detr_star(ig,l) = zdz * nu
     245                  entr_star(ig,l) = zdz * (zbetalpha * gamma / zw2m + nu)
    355246               ELSE
    356                   zdqt(ig,l) = max(zqta(ig,l-1)-po(ig,l),0.) / po(ig,l)
     247                  detr_star(ig,l) = zdz * (nu - betalpha * gamma / zw2m)
     248                  entr_star(ig,l) = zdz * nu
    357249               ENDIF
    358                
    359 !------------------------------------------------------------------------------
    360 ! Detrainment
    361 !------------------------------------------------------------------------------
    362                
    363                detr_star(ig,l) = f_star(ig,l) * zdz * (                       &
    364                &                 mix0 * 0.1 / (zalpha + 0.001)                &
    365                &               + MAX(detr_min,                                &
    366                &                 -afact * zbetalpha * zbuoyjam(ig,l) / zw2m   &
    367                &               + detr_q_coef*(zdqt(ig,l)/zw2m)**detr_q_power) )
    368250               
    369251!               IF (detr_star(ig,l).lt.0.) THEN
    … …  
    372254!               ENDIF
    373255               
    374 !------------------------------------------------------------------------------
    375 ! Entrainment
    376 !------------------------------------------------------------------------------
    377                
    378                entr_star(ig,l) = f_star(ig,l) * zdz * (                       &
    379                &                 mix0 * 0.1 / (zalpha+0.001)                  &
    380                &               + MAX(entr_min,                                &
    381                &                 zbetalpha * afact * zbuoyjam(ig,l) / zw2m    &
    382                &               - zbetalpha * fact_epsilon) )
    383                
    384256!               IF (entr_star(ig,l).lt.0.) THEN
    385257!                  print *, 'WARNING: entrainment is negative!'
    … …  
    387259!               ENDIF
    388260               
    389 !------------------------------------------------------------------------------
    390 ! Alimentation and entrainment
    391 !------------------------------------------------------------------------------
    392                
    393                IF (l.lt.lalim(ig)) THEN
    394                   alim_star(ig,l) = max(alim_star(ig,l),entr_star(ig,l))
    395                   entr_star(ig,l) = 0.
    396                ENDIF
    397                
    398 !------------------------------------------------------------------------------
    399 ! Mass flux
    400 !------------------------------------------------------------------------------
    401                
    402                f_star(ig,l+1) = f_star(ig,l) + alim_star(ig,l)                &
    403                &              + entr_star(ig,l) - detr_star(ig,l)
     261               f_star(ig,l+1) = f_star(ig,l) + entr_star(ig,l) - detr_star(ig,l)
    404262               
    405263!               IF (f_star(ig,l+1).le.0.) THEN
    … …  
    411269         ENDDO
    412270         
    413 !==============================================================================
    414 ! 5. Calcul de la vitesse verticale en melangeant Tl et qt du thermique
    415 !==============================================================================
    416          
    417          activetmp(:) = active(:) .and. f_star(:,l+1)>1.e-10
    418          
    419 !------------------------------------------------------------------------------
    420 ! Calcul du melange avec l'environnement
    421 !------------------------------------------------------------------------------
     271!-------------------------------------------------------------------------------
     272! Mixing between thermal plume and environment
     273!-------------------------------------------------------------------------------
     274         
     275         activetmp(:) = active(:) .and. (f_star(:,l+1).GT.1.e-10)
    422276         
    423277         DO ig=1,ngrid
    424278            IF (activetmp(ig)) THEN
    425                ztla(ig,l) = (f_star(ig,l) * ztla(ig,l-1)                      &  ! ztla is set to TP in plume (mixed)
    426                &          + (alim_star(ig,l) + entr_star(ig,l)) * zthl(ig,l)) &
     279               zhla(ig,l) = (f_star(ig,l) * zhla(ig,l-1)                      &  ! zhla is set to TP in plume (mixed)
     280               &          + entr_star(ig,l) * zhl(ig,l))                      &
    427281               &          / (f_star(ig,l+1) + detr_star(ig,l))
    428282               zqta(ig,l) = (f_star(ig,l) * zqta(ig,l-1) +                    &  ! zqta is set to qt in plume (mixed)
    429                &          + (alim_star(ig,l) + entr_star(ig,l)) * po(ig,l))   &
     283               &          + entr_star(ig,l) * zqt(ig,l))                      &
    430284               &          / (f_star(ig,l+1) + detr_star(ig,l))
    431285            ENDIF
    432286         ENDDO
    433287         
    434          ztemp(:) = zpspsk(:,l) * ztla(:,l)
     288         ztemp(:) = zpopsk(:,l) * zhla(:,l)
    435289         
    436290         DO ig=1,ngrid
    437291            IF (activetmp(ig)) THEN
    438                CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsatth(ig,l))
    439             ENDIF
    440          ENDDO
    441          
    442 !------------------------------------------------------------------------------
    443 ! Calcul de la vitesse verticale zw2 apres le melange
    444 !------------------------------------------------------------------------------
     292               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsa(ig,l))
     293            ENDIF
     294         ENDDO
     295         
     296!-------------------------------------------------------------------------------
     297! Vertical speed (after mixing between plume and environment)
     298!-------------------------------------------------------------------------------
    445299         
    446300         DO ig=1,ngrid
    447301            IF (activetmp(ig)) THEN
    448                zqla(ig,l) = max(0.,zqta(ig,l)-zqsatth(ig,l))                     ! zqla is set to ql   plume (mixed)
    449                ztva(ig,l) = ztla(ig,l) * zpspsk(ig,l)+RLvCp*zqla(ig,l)           ! ztva is set to TR   plume (mixed)
    450                ztva(ig,l) = ztva(ig,l) / zpspsk(ig,l)                            ! ztva is set to TRP  plume (mixed)
    451                zha(ig,l)  = ztva(ig,l)                                           ! zha  is set to TRP  plume (mixed)
    452                ztva(ig,l) = ztva(ig,l) * (1. + RETV*(zqta(ig,l)-zqla(ig,l))   &  ! ztva is set to TRPV plume (mixed)
    453                &                             - zqla(ig,l))
     302               zqla(ig,l) = max(0.,zqta(ig,l)-zqsa(ig,l))                        ! zqla is set to ql   plume (mixed)
     303               
     304               zta(ig,l) = zhla(ig,l) * zpopsk(ig,l) + RLvCp * zqla(ig,l)        ! ztva is set to TR   plume (mixed)
     305               ztva(ig,l) = zta(ig,l) / zpopsk(ig,l)                          &  ! ztva is set to TRPV plume (mixed)           
     306               &          * (1. + RETV*(zqta(ig,l)-zqla(ig,l)) - zqla(ig,l))
    454307               
    455308               zbuoy(ig,l) = RG * (ztva(ig,l) - ztv(ig,l)) / ztv(ig,l)
    456309               zdz = zlev(ig,l+1) - zlev(ig,l)
    457310               
    458 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    459 ! AB : initial formula
    460 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     311!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     312! AB: initial formula
     313!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    461314!               zw2fact = fact_epsilon * 2. * zdz / (1. + betalpha)
    462315!               zdw2 = afact * zbuoy(ig,l) / fact_epsilon
    463316!               zw2(ig,l+1) = Max(0.0001,exp(-zw2fact)*(zw2(ig,l)-zdw2)+zdw2)
    464 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    465 ! AB : own derivation for zw2 (Rio 2010 formula)
    466 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     317!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     318! AB: own derivation for zw2 (Rio et al. 2010)
     319!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    467320!               zw2fact = 2. * (fact_epsilon * zdz  + entr_star(ig,l) / f_star(ig,l))
    468321!               zdw2 = 2. * afact * zbuoy(ig,l) * zdz
    … …  
    470323               zdw2 = 2. * afact * zbuoy(ig,l) * zdz / (1. + betalpha)
    471324               zw2(ig,l+1) = Max(0., exp(-zw2fact) * zw2(ig,l) + zdw2)
    472                
    473 !               IF (zw2(ig,l+1).le.0.) THEN
    474 !                  print *, 'WARNING: zw2 is negative!'
    475 !                  print *, 'l,zw2', l+1, zw2(ig,l+1)
    476 !               ENDIF
    477             ENDIF
    478          ENDDO
    479          
    480       ENDDO
    481      
    482 !==============================================================================
    483 ! 6. New computation of alim_star_tot
    484 !==============================================================================
    485      
    486       DO ig=1,ngrid
    487          alim_star_tot(ig) = 0.
    488       ENDDO
    489      
    490       DO ig=1,ngrid
    491          DO l=1,lalim(ig)-1
    492             alim_star_tot(ig) = alim_star_tot(ig) + alim_star(ig,l)
    493          ENDDO
     325            ENDIF
     326         ENDDO
     327         
    494328      ENDDO
    495329     
Note: See TracChangeset for help on using the changeset viewer.