Ignore:
Timestamp:
Feb 15, 2019, 2:43:57 PM (6 years ago)
Author:
aboissinot
Message:

Fix a bug in thermcell_alim.F90 where vertical and horizontal loops must be inverted.
In thermal plume model, arrays size declaration is standardised (no longer done thanks to dimphy module but by way of arguments).
Clean up some thermal plume model routines (remove uselesss variables...)

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

Legend:

Unmodified
Added
Removed

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

    r2069 r2101  
    267267      real entr0(ngrid, nlayer)           ! Entrainment
    268268      real detr0(ngrid, nlayer)           ! Detrainment
    269       real zmax0(ngrid)                   ! Plume height
    270269      real fraca(ngrid, nlayer+1)         ! Fraction of the surface that plumes occupies
    271270      real zw2(ngrid, nlayer+1)           ! Squared vertical speed or its square root
    … …  
    11951194                             f0, fm0, entr0, detr0,                              &
    11961195                             zqta, zqla, ztv, ztva, ztla, zthl, zqsatth,         &
    1197                              zmax0, zw2, fraca,                                  &
     1196                             zw2, fraca,                                         &
    11981197                             lmin, lmix, lalim, lmax,                            &
    11991198                             zpopsk, ratqscth, ratqsdiff,                        &

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

    r2093 r2101  
    22!
    33!
    4       SUBROUTINE thermcell_alim(ngrid,klev,ztv2,zlev,alim_star,lalim,lmin)
     4SUBROUTINE thermcell_alim(ngrid,nlay,ztv2,zlev,alim_star,lalim,lmin)
    55     
    6      
    7       USE thermcell_mod, ONLY: linf, d_temp
    8      
    9       IMPLICIT NONE
    106     
    117!==============================================================================
    … …  
    1410! laterale a la base des panaches thermiques
    1511!==============================================================================
     12     
     13      IMPLICIT NONE
    1614     
    1715     
    … …  
    2422     
    2523      INTEGER, INTENT(IN) :: ngrid
    26       INTEGER, INTENT(IN) :: klev
    27       INTEGER, INTENT(IN) :: lmin(ngrid)        ! plume initial level
     24      INTEGER, INTENT(IN) :: nlay
     25      INTEGER, INTENT(IN) :: lmin(ngrid)           ! plume initial level
    2826     
    29       REAL, INTENT(IN) :: ztv2(ngrid,klev)      ! Virtual potential temperature
    30       REAL, INTENT(IN) :: zlev(ngrid,klev+1)    ! levels altitude
     27      REAL, INTENT(IN) :: ztv2(ngrid,nlay)         ! Virtual potential temperature
     28      REAL, INTENT(IN) :: zlev(ngrid,nlay+1)       ! levels altitude
    3129     
    3230!      outputs:
    3331!      --------
    3432     
    35       INTEGER, INTENT(OUT) :: lalim(ngrid)      ! alimentation maximal level
     33      INTEGER, INTENT(OUT) :: lalim(ngrid)         ! alimentation maximal level
    3634     
    37       REAL, INTENT(OUT) :: alim_star(ngrid,klev)
     35      REAL, INTENT(OUT) :: alim_star(ngrid,nlay)   ! Normalized alimentation
    3836     
    3937!      local:
    … …  
    4240      INTEGER :: ig, l
    4341     
    44       REAL :: alim_star_tot(ngrid)              ! integrated alimentation
     42      REAL :: alim_star_tot(ngrid)                 ! integrated alimentation
    4543     
    4644!==============================================================================
    … …  
    5654!==============================================================================
    5755     
    58       DO l=lmin(ig),klev-1
    59          DO ig=1,ngrid
    60             IF ((ztv2(ig,l)>ztv2(ig,l+1)).and.(ztv2(ig,lmin(ig))>=ztv2(ig,l))) THEN
     56      DO ig=1,ngrid
     57         DO l=lmin(ig),nlay-1
     58            IF ((ztv2(ig,l).gt.ztv2(ig,l+1)).and.(ztv2(ig,lmin(ig)).ge.ztv2(ig,l))) THEN
    6159               alim_star(ig,l) = MAX( (ztv2(ig,l) - ztv2(ig,l+1)), 0.)  &
    6260               &               * sqrt(zlev(ig,l+1))
    … …  
    7371!------------------------------------------------------------------------------
    7472     
    75       DO l=1,klev
     73      DO l=1,nlay
    7674         DO ig=1,ngrid
    7775            IF (alim_star_tot(ig) > 1.e-10 ) THEN

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

    r2060 r2101  
    22!
    33!
    4       SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep  &
    5      &                  ,pplay,pplev  &
    6      &                  ,fm0,entr0,lmax  &
    7      &                  ,ale_bl,alp_bl,lalim_conv,wght_th &
    8      &                  ,zw2,fraca &
    9 !!! ncessaire en plus
    10      &                  ,pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax &
    11 !!! nrlmd le 10/04/2012
    12      &                  ,pbl_tke,pctsrf,omega,airephy &
    13      &                  ,zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
    14      &                  ,n2,s2,ale_bl_stat &
    15      &                  ,therm_tke_max,env_tke_max &
    16      &                  ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
    17      &                  ,alp_bl_conv,alp_bl_stat )
    18 !!! fin nrlmd le 10/04/2012
    19 
    20       USE dimphy
    21 ! AB : why use dimphy to get klon, klev while subroutine arguments ngrid, nlay are the same?
    22       USE thermcell_mod
    23      
    24       IMPLICIT NONE
    25 
    26 !=======================================================================
     4SUBROUTINE thermcell_alp(ngrid,nlay,ptimestep,pplay,pplev,                    &
     5                         fm0,entr0,lmax,ale_bl,alp_bl,lalim_conv,             &
     6                         wght_th,zw2,fraca,                                   &
     7                         pcon,rhobarz,wth3,wmax_sec,lalim,fm,alim_star,zmax,  &
     8                         pbl_tke,pctsrf,omega,airephy,                        &
     9                         zlcl,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0,   &
     10                         n2,s2,ale_bl_stat,                                   &
     11                         therm_tke_max,env_tke_max,                           &
     12                         alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke,          &
     13                         alp_bl_conv,alp_bl_stat)
     14     
     15     
     16!==============================================================================
    2717!   Auteurs: Frederic Hourdin, Catherine Rio, Anne Mathieu
    2818!   Version du 09.02.07
    … …  
    5040!       15, 16 run with impl=-1 : numerical convergence with NPv3
    5141!       17, 18 run with impl=1  : more stable
    52 !    15 and 17 correspond to the activation of the stratocumulus "bidouille"
    53 !
    54 !=======================================================================
    55 
    56 !-----------------------------------------------------------------------
    57 !   declarations:
    58 !   -------------
    59 
    60 !   arguments:
    61 !   ----------
    62 
    63 !IM 140508
    64 
    65       INTEGER ngrid,nlay
    66       real ptimestep
    67       REAL pplay(ngrid,nlay),pplev(ngrid,nlay+1)
    68 
    69 !   local:
    70 !   ------
    71 
    72 
    73       REAL susqr2pi, reuler
    74 
    75       INTEGER ig,k,l
    76       INTEGER lmax(klon),lalim(klon)
    77       real zmax(klon),zw2(klon,klev+1)
    78 
    79 !on garde le zmax du pas de temps precedent
    80 
    81 
    82       real fraca(klon,klev+1)
    83       real wth3(klon,klev)
     42!       15 and 17 correspond to the activation of the stratocumulus "bidouille"
     43!
     44!==============================================================================
     45     
     46      USE thermcell_mod
     47     
     48      IMPLICIT NONE
     49     
     50     
     51!==============================================================================
     52! Declarations
     53!==============================================================================
     54     
     55!      Inputs:
     56!      -------
     57     
     58      INTEGER,INTENT(in) :: ngrid
     59      INTEGER,INTENT(in) :: nlay
     60      INTEGER lmax(ngrid)
     61      INTEGER lalim(ngrid)
     62     
     63      REAL,INTENT(in) :: ptimestep
     64      REAL,INTENT(in) :: pplay(ngrid,nlay)
     65      REAL,INTENT(in) :: pplev(ngrid,nlay+1)
     66! On garde le zmax du pas de temps precedent
     67      REAL zmax(ngrid)
     68      REAL zw2(ngrid,nlay+1)
     69      REAL fraca(ngrid,nlay+1)
     70      real fm0(ngrid,nlay+1)
     71      real entr0(ngrid,nlay)
     72      real pbl_tke(ngrid,nlay+1,nbsrf)
     73      real pctsrf(ngrid,nbsrf)
     74      real omega(ngrid,nlay)
     75      real airephy(ngrid)
     76      real alim_star(ngrid,nlay)
     77     
     78!      Outputs:
     79!     ---------
     80     
     81      REAL susqr2pi
     82      REAL reuler
     83      real wth3(ngrid,nlay)
    8484! FH probleme de dimensionnement avec l'allocation dynamique
    8585!     common/comtherm/thetath2,wth2
    86       real rhobarz(klon,klev)
    87 
    88       real wmax_sec(klon)
    89       real fm0(klon,klev+1),entr0(klon,klev)
    90       real fm(klon,klev+1)
    91 
    92 !niveau de condensation
    93       real pcon(klon)
    94 
    95       real alim_star(klon,klev)
    96 
    97 !!! nrlmd le 10/04/2012
    98 
    99 !------Entrées
    100       real pbl_tke(klon,klev+1,nbsrf)
    101       real pctsrf(klon,nbsrf)
    102       real omega(klon,klev)
    103       real airephy(klon)
    104 !------Sorties
    105       real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
    106       real therm_tke_max0(klon),env_tke_max0(klon)
    107       real n2(klon),s2(klon)
    108       real ale_bl_stat(klon)
    109       real therm_tke_max(klon,klev),env_tke_max(klon,klev)
    110       real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
    111 !------Local
     86      real rhobarz(ngrid,nlay)
     87      real wmax_sec(ngrid)
     88      real fm(ngrid,nlay+1)
     89     
     90! niveau de condensation
     91      real pcon(ngrid)
     92     
     93      real zlcl(ngrid)
     94      real fraca0(ngrid)
     95      real w0(ngrid)
     96      real w_conv(ngrid)
     97      real therm_tke_max0(ngrid)
     98      real env_tke_max0(ngrid)
     99      real n2(ngrid)
     100      real s2(ngrid)
     101      real ale_bl_stat(ngrid)
     102      real therm_tke_max(ngrid,nlay)
     103      real env_tke_max(ngrid,nlay)
     104      real alp_bl_det(ngrid)
     105      real alp_bl_fluct_m(ngrid)
     106      real alp_bl_fluct_tke(ngrid)
     107      real alp_bl_conv(ngrid)
     108      real alp_bl_stat(ngrid)
     109     
     110!      local:
     111!      ------
     112     
     113      INTEGER ig,k,l
    112114      integer nsrf
    113       real rhobarz0(klon)                    ! Densité au LCL
    114       logical ok_lcl(klon)                   ! Existence du LCL des thermiques
    115       integer klcl(klon)                     ! Niveau du LCL
    116       real interp(klon)                      ! Coef d'interpolation pour le LCL
     115      integer klcl(ngrid)                     ! Niveau du LCL
     116     
     117      real rhobarz0(ngrid)                    ! Densité au LCL
     118      logical ok_lcl(ngrid)                   ! Existence du LCL des thermiques
     119      real interp(ngrid)                      ! Coef d'interpolation pour le LCL
    117120!--Triggering
    118       real Su                                ! Surface unité: celle d'un updraft élémentaire
    119       parameter(Su=4e4)
    120       real hcoef                             ! Coefficient directeur pour le calcul de s2
    121       parameter(hcoef=1)
    122       real hmincoef                          ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
    123       parameter(hmincoef=0.3)
    124       real eps1                              ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
    125       parameter(eps1=0.3)
    126       real hmin(ngrid)                       ! Ordonnée à l'origine pour le calcul de s2
    127       real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
    128       real zmax_moy_coef
    129       parameter(zmax_moy_coef=0.33)
    130       real depth(klon)                       ! Epaisseur moyenne du cumulus
    131       real w_max(klon)                       ! Vitesse max statistique
    132       real s_max(klon)
     121      real,parameter :: Su = 4.e4             ! Surface unité: celle d'un updraft élémentaire
     122      real,parameter :: hcoef = 1.            ! Coefficient directeur pour le calcul de s2
     123      real,parameter :: hmincoef = 0.3        ! Coefficient directeur pour l'ordonnée à l'origine pour le calcul de s2
     124      real,parameter :: eps1 = 0.3            ! Fraction de surface occupée par la population 1 : eps1=n1*s1/(fraca0*Sd)
     125      real,parameter :: zmax_moy_coef=0.33
     126      real hmin(ngrid)                        ! Ordonnée à l'origine pour le calcul de s2
     127      real zmax_moy(ngrid)                    ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
     128      real depth(ngrid)                       ! Epaisseur moyenne du cumulus
     129      real w_max(ngrid)                       ! Vitesse max statistique
     130      real s_max(ngrid)
    133131!--Closure
    134       real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne
    135       real pbl_tke_max0(klon)                ! TKE moyenne au LCL
    136       real w_ls(klon,klev)                   ! Vitesse verticale grande échelle (m/s)
    137       real coef_m                            ! On considère un rendement pour alp_bl_fluct_m
    138       parameter(coef_m=1.)
    139       real coef_tke                          ! On considère un rendement pour alp_bl_fluct_tke
    140       parameter(coef_tke=1.)
    141 
    142 !!! fin nrlmd le 10/04/2012
    143 
    144 !
    145       !nouvelles variables pour la convection
    146       real ale_bl(klon)
    147       real alp_bl(klon)
    148       real alp_int(klon),dp_int(klon),zdp
    149       real fm_tot(klon)
    150       real wght_th(klon,klev)
    151       integer lalim_conv(klon)
     132      real pbl_tke_max(ngrid,nlay)            ! Profil de TKE moyenne
     133      real pbl_tke_max0(ngrid)                ! TKE moyenne au LCL
     134      real w_ls(ngrid,nlay)                   ! Vitesse verticale grande échelle (m/s)
     135      real,parameter :: coef_m=1.             ! On considère un rendement pour alp_bl_fluct_m
     136      real,parameter :: coef_tke = 1.         ! On considère un rendement pour alp_bl_fluct_tke
     137     
     138! nouvelles variables pour la convection
     139      real ale_bl(ngrid)
     140      real alp_bl(ngrid)
     141      real alp_int(ngrid),dp_int(ngrid),zdp
     142      real fm_tot(ngrid)
     143      real wght_th(ngrid,nlay)
     144      integer lalim_conv(ngrid)
    152145!v1d     logical therm
    153146!v1d     save therm
    154147     
    155      
    156 !------------------------------------------------------------
    157 !  Initialize output arrays related to stochastic triggering
    158 !------------------------------------------------------------
    159      
    160       DO ig = 1,klon
     148!==============================================================================
     149! Initialization
     150!==============================================================================
     151     
     152      DO ig = 1,ngrid
    161153         zlcl(ig) = 0.
    162154         fraca0(ig) = 0.
    … …  
    175167      ENDDO
    176168     
    177       DO l = 1,klev
    178          DO ig = 1,klon
     169      DO l = 1,nlay
     170         DO ig = 1,ngrid
    179171            therm_tke_max(ig,l) = 0.
    180172            env_tke_max(ig,l) = 0.
    … …  
    185177!------------Test sur le LCL des thermiques
    186178      do ig=1,ngrid
    187       ok_lcl(ig)=.false.
    188       if ( (pcon(ig) .gt. pplay(ig,klev-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
    189    enddo
    190 
     179         ok_lcl(ig)=.false.
     180         if ( (pcon(ig) .gt. pplay(ig,nlay-1)) .and. (pcon(ig) .lt. pplay(ig,1)) ) ok_lcl(ig)=.true.
     181      enddo
     182     
    191183!------------Localisation des niveaux entourant le LCL et du coef d'interpolation
    192     do l=1,nlay-1
    193       do ig=1,ngrid
    194         if (ok_lcl(ig)) then
     184      do l=1,nlay-1
     185         do ig=1,ngrid
     186            if (ok_lcl(ig)) then
    195187!ATTENTION,zw2 calcule en pplev
    196 !          if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
    197 !          klcl(ig)=l
    198 !          interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
    199 !          endif
    200           if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
    201           klcl(ig)=l
    202           interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
    203           endif
    204         endif
    205       enddo
    206     enddo
    207 
     188!               if ((pplay(ig,l) .ge. pcon(ig)) .and. (pplay(ig,l+1) .le. pcon(ig))) then
     189!                  klcl(ig)=l
     190!                  interp(ig)=(pcon(ig)-pplay(ig,klcl(ig)))/(pplay(ig,klcl(ig)+1)-pplay(ig,klcl(ig)))
     191!               endif
     192               if ((pplev(ig,l) .ge. pcon(ig)) .and. (pplev(ig,l+1) .le. pcon(ig))) then
     193                  klcl(ig)=l
     194                  interp(ig)=(pcon(ig)-pplev(ig,klcl(ig)))/(pplev(ig,klcl(ig)+1)-pplev(ig,klcl(ig)))
     195               endif
     196            endif
     197         enddo
     198      enddo
     199     
    208200!------------Hauteur des thermiques
     201!jyg le 27/04/2012
     202!      do ig =1,ngrid
     203!         rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
     204!         &               -rhobarz(ig,klcl(ig)))*interp(ig)
     205!         zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
     206!         if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
     207!!     enddo
     208     
     209      do ig =1,ngrid
     210!CR:REHABILITATION ZMAX CONTINU
     211         if (ok_lcl(ig)) then
     212            rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
     213            &               -rhobarz(ig,klcl(ig)))*interp(ig)
     214            zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
     215            zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
     216         else
     217            rhobarz0(ig)=0.
     218            zlcl(ig)=zmax(ig)
     219         endif
     220      enddo
     221!!jyg fin
     222     
     223!------------Calcul des propriétés du thermique au LCL
     224      IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
     225         
     226!-----Initialisation de la TKE moyenne
     227         do l=1,nlay
     228            do ig=1,ngrid
     229               pbl_tke_max(ig,l)=0.
     230            enddo
     231         enddo
     232         
     233!-----Calcul de la TKE moyenne
     234         do nsrf=1,nbsrf
     235            do l=1,nlay
     236               do ig=1,ngrid
     237                  pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
     238               enddo
     239            enddo
     240         enddo
     241         
     242!-----Initialisations des TKE dans et hors des thermiques
     243         do l=1,nlay
     244            do ig=1,ngrid
     245               therm_tke_max(ig,l)=pbl_tke_max(ig,l)
     246               env_tke_max(ig,l)=pbl_tke_max(ig,l)
     247            enddo
     248         enddo
     249         
     250!-----Calcul de la TKE transportée par les thermiques : therm_tke_max
     251         call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
     252         &           rg,pplev,therm_tke_max)
     253         
     254!         print *,' thermcell_tke_transport -> '   !!jyg
     255         
     256!-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
     257         do l=1,nlay
     258            do ig=1,ngrid
     259               pbl_tke_max(ig,l) = fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l)    !  Recalcul de TKE moyenne aprés transport de TKE_TH
     260               env_tke_max(ig,l) = (pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l))  !  Recalcul de TKE dans  l'environnement aprés transport de TKE_TH
     261               w_ls(ig,l) = -1.*omega(ig,l)/(RG*rhobarz(ig,l))                                           !  Vitesse verticale de grande échelle
     262            enddo
     263         enddo
     264         
     265!         print *,' apres w_ls = '   !!jyg
     266         
     267         do ig=1,ngrid
     268            if (ok_lcl(ig)) then
     269               fraca0(ig) = fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1)          &
     270               &          - fraca(ig,klcl(ig)))*interp(ig)
     271               w0(ig) = zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1)                  &
     272               &      - zw2(ig,klcl(ig)))*interp(ig)
     273               w_conv(ig) = w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1)            &
     274               &          - w_ls(ig,klcl(ig)))*interp(ig)
     275               therm_tke_max0(ig) = therm_tke_max(ig,klcl(ig))                &
     276               &                  + (therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
     277               env_tke_max0(ig) = env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
     278               &                - env_tke_max(ig,klcl(ig)))*interp(ig)
     279               pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
     280               &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
     281               
     282               if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig) = 20.
     283               
     284               if (env_tke_max0(ig).ge.20.) env_tke_max0(ig) = 20.
     285               
     286               if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig) = 20.
     287            else
     288               fraca0(ig) = 0.
     289               w0(ig) = 0.
     290! jyg le 27/04/2012
     291!               zlcl(ig) = 0.
     292            endif
     293         enddo
     294     
     295      ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
     296     
     297!      print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
     298     
     299!------------Triggering------------------
     300      IF (iflag_trig_bl.ge.1) THEN
     301         
     302!-----Initialisations
     303         depth(:)=0.
     304         n2(:)=0.
     305         s2(:)=100. ! some low value, arbitrary
     306         s_max(:)=0.
     307         
     308!-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
     309         do ig=1,ngrid
     310            zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
     311            depth(ig)=zmax_moy(ig)-zlcl(ig)
     312            hmin(ig)=hmincoef*zlcl(ig)
     313           
     314            if (depth(ig).ge.10.) then
     315               s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
     316               n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
     317               
    209318!!jyg le 27/04/2012
    210 !!    do ig =1,ngrid
    211 !!    rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
    212 !!    &               -rhobarz(ig,klcl(ig)))*interp(ig)
    213 !!    zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
    214 !!      if ( (.not.ok_lcl(ig)) .or. (zlcl(ig).gt.zmax(ig)) ) zlcl(ig)=zmax(ig) ! Si zclc > zmax alors on pose zlcl = zmax
    215 !!    enddo
    216     do ig =1,ngrid
    217 !CR:REHABILITATION ZMAX CONTINU
    218       if (ok_lcl(ig)) then
    219          rhobarz0(ig)=rhobarz(ig,klcl(ig))+(rhobarz(ig,klcl(ig)+1) &
    220          &               -rhobarz(ig,klcl(ig)))*interp(ig)
    221          zlcl(ig)=(pplev(ig,1)-pcon(ig))/(rhobarz0(ig)*RG)
    222          zlcl(ig)=min(zlcl(ig),zmax(ig))   ! Si zlcl > zmax alors on pose zlcl = zmax
    223       else
    224          rhobarz0(ig)=0.
    225          zlcl(ig)=zmax(ig)
    226       endif
    227    enddo
    228 !!jyg fin
    229 
    230 !------------Calcul des propriétés du thermique au LCL
    231   IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) THEN
    232 
    233 !-----Initialisation de la TKE moyenne
    234    do l=1,nlay
    235       do ig=1,ngrid
    236          pbl_tke_max(ig,l)=0.
    237       enddo
    238    enddo
    239 
    240 !-----Calcul de la TKE moyenne
    241    do nsrf=1,nbsrf
    242       do l=1,nlay
    243          do ig=1,ngrid
    244             pbl_tke_max(ig,l)=pctsrf(ig,nsrf)*pbl_tke(ig,l,nsrf)+pbl_tke_max(ig,l)
    245          enddo
    246       enddo
    247    enddo
    248 
    249 !-----Initialisations des TKE dans et hors des thermiques
    250    do l=1,nlay
    251       do ig=1,ngrid
    252          therm_tke_max(ig,l)=pbl_tke_max(ig,l)
    253          env_tke_max(ig,l)=pbl_tke_max(ig,l)
    254       enddo
    255    enddo
    256 
    257 !-----Calcul de la TKE transportée par les thermiques : therm_tke_max
    258    call thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,  &
    259    &           rg,pplev,therm_tke_max)
    260    
    261 !   print *,' thermcell_tke_transport -> '   !!jyg
    262 
    263 !-----Calcul des profils verticaux de TKE hors thermiques : env_tke_max, et de la vitesse verticale grande échelle : W_ls
    264    do l=1,nlay
    265       do ig=1,ngrid
    266          pbl_tke_max(ig,l)=fraca(ig,l)*therm_tke_max(ig,l)+(1.-fraca(ig,l))*env_tke_max(ig,l)         !  Recalcul de TKE moyenne aprés transport de TKE_TH
    267          env_tke_max(ig,l)=(pbl_tke_max(ig,l)-fraca(ig,l)*therm_tke_max(ig,l))/(1.-fraca(ig,l))       !  Recalcul de TKE dans  l'environnement aprés transport de TKE_TH
    268          w_ls(ig,l)=-1.*omega(ig,l)/(RG*rhobarz(ig,l))                                                !  Vitesse verticale de grande échelle
    269       enddo
    270    enddo
    271    
    272 !   print *,' apres w_ls = '   !!jyg
    273    
    274    do ig=1,ngrid
    275       if (ok_lcl(ig)) then
    276          fraca0(ig)=fraca(ig,klcl(ig))+(fraca(ig,klcl(ig)+1) &
    277          &             -fraca(ig,klcl(ig)))*interp(ig)
    278          w0(ig)=zw2(ig,klcl(ig))+(zw2(ig,klcl(ig)+1) &
    279          &         -zw2(ig,klcl(ig)))*interp(ig)
    280          w_conv(ig)=w_ls(ig,klcl(ig))+(w_ls(ig,klcl(ig)+1) &
    281          &             -w_ls(ig,klcl(ig)))*interp(ig)
    282          therm_tke_max0(ig)=therm_tke_max(ig,klcl(ig)) &
    283          &                     +(therm_tke_max(ig,klcl(ig)+1)-therm_tke_max(ig,klcl(ig)))*interp(ig)
    284          env_tke_max0(ig)=env_tke_max(ig,klcl(ig))+(env_tke_max(ig,klcl(ig)+1) &
    285          &                   -env_tke_max(ig,klcl(ig)))*interp(ig)
    286          pbl_tke_max0(ig)=pbl_tke_max(ig,klcl(ig))+(pbl_tke_max(ig,klcl(ig)+1) &
    287          &                   -pbl_tke_max(ig,klcl(ig)))*interp(ig)
    288          
    289          if (therm_tke_max0(ig).ge.20.) therm_tke_max0(ig)=20.
    290          
    291          if (env_tke_max0(ig).ge.20.) env_tke_max0(ig)=20.
    292          
    293          if (pbl_tke_max0(ig).ge.20.) pbl_tke_max0(ig)=20.
    294       else
    295          fraca0(ig)=0.
    296          w0(ig)=0.
    297 ! jyg le 27/04/2012
    298 !         zlcl(ig)=0.
    299       endif
    300    enddo
    301    
    302    ENDIF ! IF ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) )
    303 
    304 !   print *,'ENDIF  ( (iflag_trig_bl.ge.1) .or. (iflag_clos_bl.ge.1) ) '    !!jyg
    305 
    306 !------------Triggering------------------
    307   IF (iflag_trig_bl.ge.1) THEN
    308 
    309 !-----Initialisations
    310    depth(:)=0.
    311    n2(:)=0.
    312    s2(:)=100. ! some low value, arbitrary
    313    s_max(:)=0.
    314 
    315 !-----Epaisseur du nuage (depth) et détermination de la queue du spectre de panaches (n2,s2) et du panache le plus gros (s_max)
    316    do ig=1,ngrid
    317       zmax_moy(ig)=zlcl(ig)+zmax_moy_coef*(zmax(ig)-zlcl(ig))
    318       depth(ig)=zmax_moy(ig)-zlcl(ig)
    319       hmin(ig)=hmincoef*zlcl(ig)
    320      
    321       if (depth(ig).ge.10.) then
    322          s2(ig)=(hcoef*depth(ig)+hmin(ig))**2
    323          n2(ig)=(1.-eps1)*fraca0(ig)*airephy(ig)/s2(ig)
    324          
    325 !!jyg le 27/04/2012
    326 !!        s_max(ig)=s2(ig)*log(n2(ig))
    327 !!        if (n2(ig) .lt. 1) s_max(ig)=0.
    328           s_max(ig)=s2(ig)*log(max(n2(ig),1.))
     319!!              s_max(ig) = s2(ig)*log(n2(ig))
     320!!              if (n2(ig) .lt. 1) s_max(ig) = 0.
     321               s_max(ig) = s2(ig)*log(max(n2(ig),1.))
    329322!!fin jyg
    330       else
    331          n2(ig)=0.
    332          s_max(ig)=0.
    333       endif
    334    enddo
    335    
    336 !   print *,'avant Calcul de Wmax '    !!jyg
    337 
     323            else
     324               n2(ig) = 0.
     325               s_max(ig) = 0.
     326            endif
     327         enddo
     328        
     329!         print *,'avant Calcul de Wmax '    !!jyg
     330         
    338331!-----Calcul de Wmax et ALE_BL_STAT associée
    339332!!jyg le 30/04/2012
    … …  
    347340!!     endif
    348341!!   enddo
    349    
    350    susqr2pi=su*sqrt(2.*Rpi)
    351    reuler=exp(1.)
    352    
    353    do ig=1,ngrid
    354       if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
    355          w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
    356          ale_bl_stat(ig)=0.5*w_max(ig)**2
    357       else
    358          w_max(ig)=0.
    359          ale_bl_stat(ig)=0.
    360       endif
    361    enddo
    362 
    363    ENDIF ! iflag_trig_bl
    364    
    365 !   print *,'ENDIF  iflag_trig_bl'    !!jyg
    366 
     342        
     343         susqr2pi=su*sqrt(2.*Rpi)
     344         reuler=exp(1.)
     345        
     346         do ig=1,ngrid
     347            if ( (depth(ig).ge.10.) .and. (s_max(ig).gt.susqr2pi*reuler) ) then
     348               w_max(ig)=w0(ig)*(1.+sqrt(2.*log(s_max(ig)/susqr2pi)-log(2.*log(s_max(ig)/susqr2pi))))
     349               ale_bl_stat(ig)=0.5*w_max(ig)**2
     350            else
     351               w_max(ig)=0.
     352               ale_bl_stat(ig)=0.
     353            endif
     354         enddo
     355         
     356      ENDIF ! iflag_trig_bl
     357      
     358!      print *,'ENDIF  iflag_trig_bl'    !!jyg
     359     
    367360!------------Closure------------------
    368 
    369    IF (iflag_clos_bl.ge.2) THEN
    370      
     361     
     362      IF (iflag_clos_bl.ge.2) THEN
     363         
    371364!-----Calcul de ALP_BL_STAT
    372       do ig=1,ngrid
    373          alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
    374          alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
    375          &                   (w0(ig)**2)
    376          alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
    377          &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
    378          
    379          if (iflag_clos_bl.ge.2) then
    380             alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
     365         do ig=1,ngrid
     366            alp_bl_det(ig)=0.5*coef_m*rhobarz0(ig)*(w0(ig)**3)*fraca0(ig)*(1.-2.*fraca0(ig))/((1.-fraca0(ig))**2)
     367            alp_bl_fluct_m(ig)=1.5*rhobarz0(ig)*fraca0(ig)*(w_conv(ig)+coef_m*w0(ig))* &
    381368            &                   (w0(ig)**2)
    382          else
    383             alp_bl_conv(ig)=0.
    384          endif
    385          
    386          alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
    387       enddo
    388      
     369            alp_bl_fluct_tke(ig)=3.*coef_m*rhobarz0(ig)*w0(ig)*fraca0(ig)*(therm_tke_max0(ig)-env_tke_max0(ig)) &
     370            &                    +3.*rhobarz0(ig)*w_conv(ig)*pbl_tke_max0(ig)
     371           
     372            if (iflag_clos_bl.ge.2) then
     373               alp_bl_conv(ig)=1.5*coef_m*rhobarz0(ig)*fraca0(ig)*(fraca0(ig)/(1.-fraca0(ig)))*w_conv(ig)* &
     374               &                   (w0(ig)**2)
     375            else
     376               alp_bl_conv(ig)=0.
     377            endif
     378           
     379            alp_bl_stat(ig)=alp_bl_det(ig)+alp_bl_fluct_m(ig)+alp_bl_fluct_tke(ig)+alp_bl_conv(ig)
     380         enddo
     381         
    389382!-----Sécurité ALP infinie
    390       do ig=1,ngrid
    391          if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
    392       enddo
    393      
    394   ENDIF ! (iflag_clos_bl.ge.2)
    395 
     383         do ig=1,ngrid
     384            if (fraca0(ig).gt.0.98) alp_bl_stat(ig)=2.
     385         enddo
     386         
     387      ENDIF ! (iflag_clos_bl.ge.2)
     388     
    396389!!! fin nrlmd le 10/04/2012
    397 
    398 !   print*,'avant calcul ale et alp'
    399 
    400 !calcul de ALE et ALP pour la convection
    401    alp_bl(:)=0.
    402    ale_bl(:)=0.
    403      
    404 !   print*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
    405 
    406    do l=1,nlay
    407       do ig=1,ngrid
    408          alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
    409          ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
    410 !         print*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
    411       enddo
    412    enddo
    413 
    414 !   ale sec (max de wmax/2 sous la zone d'inhibition) dans
    415 !   le cas iflag_trig_bl=3
    416    
    417    IF (iflag_trig_bl==3) then
    418       ale_bl(:)=0.5*wmax_sec(:)**2
    419    ENDIF
    420 
     390     
     391!      print*,'avant calcul ale et alp'
     392     
     393! calcul de ALE et ALP pour la convection
     394      alp_bl(:)=0.
     395      ale_bl(:)=0.
     396     
     397!      print*,'ALE,ALP ,l,zw2(ig,l),ale_bl(ig),alp_bl(ig)'
     398     
     399      do l=1,nlay
     400         do ig=1,ngrid
     401            alp_bl(ig)=max(alp_bl(ig),0.5*rhobarz(ig,l)*wth3(ig,l) )
     402            ale_bl(ig)=max(ale_bl(ig),0.5*zw2(ig,l)**2)
     403!            print*,'ALE,ALP',l,zw2(ig,l),ale_bl(ig),alp_bl(ig)
     404         enddo
     405      enddo
     406     
     407! ale sec (max de wmax/2 sous la zone d'inhibition) dans le cas iflag_trig_bl=3
     408     
     409      IF (iflag_trig_bl==3) then
     410         ale_bl(:)=0.5*wmax_sec(:)**2
     411      ENDIF
     412     
    421413!test:calcul de la ponderation des couches pour KE
    422414!initialisations
    423 
     415     
    424416      fm_tot(:)=0.
    425417      wght_th(:,:)=1.
    426418      lalim_conv(:)=lalim(:)
    427 
    428       do k=1,klev
     419     
     420      do k=1,nlay
    429421         do ig=1,ngrid
    430422            if (k<=lalim_conv(ig)) fm_tot(ig)=fm_tot(ig)+fm(ig,k)
    431423         enddo
    432424      enddo
    433 
     425     
    434426! assez bizarre car, si on est dans la couche d'alim et que alim_star et
    435427! plus petit que 1.e-10, on prend wght_th=1.
    436       do k=1,klev
     428      do k=1,nlay
    437429         do ig=1,ngrid
    438430            if (k<=lalim_conv(ig).and.alim_star(ig,k)>1.e-10) then
    … …  
    441433         enddo
    442434      enddo
    443 
     435     
    444436!      print*,'apres wght_th'
    445437!test pour prolonger la convection
    446438      do ig=1,ngrid
    447 !v1d  if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
    448       if ((alim_star(ig,1).lt.1.e-10)) then
    449       lalim_conv(ig)=1
    450       wght_th(ig,1)=1.
    451 !      print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
    452       endif
    453       enddo
    454 
     439!v1d      if ((alim_star(ig,1).lt.1.e-10).and.(therm)) then
     440         if ((alim_star(ig,1).lt.1.e-10)) then
     441            lalim_conv(ig)=1
     442            wght_th(ig,1)=1.
     443!            print*,'lalim_conv ok',lalim_conv(ig),wght_th(ig,1)
     444         endif
     445      enddo
     446     
    455447!------------------------------------------------------------------------
    456448! Modif CR/FH 20110310 : alp integree sur la verticale.
    … …  
    459451! couches
    460452!------------------------------------------------------------------------
    461 
     453     
    462454      alp_int(:)=0.
    463455      dp_int(:)=0.
     456     
    464457      do l=2,nlay
    465458        do ig=1,ngrid
    … …  
    471464        enddo
    472465      enddo
    473 
     466     
    474467      if (iflag_coupl>=3 .and. iflag_coupl<=5) then
    475468         do ig=1,ngrid
    … …  
    480473         enddo
    481474      endif
    482 
    483 
     475     
    484476! Facteur multiplicatif sur alp_bl
    485477      alp_bl(:)=alp_bl_k*alp_bl(:)
    486 
    487 !------------------------------------------------------------------------
    488 
    489 
    490 
    491       return
    492       end
     478     
     479     
     480return
     481end

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

    r2092 r2101  
    22!
    33!
    4       SUBROUTINE thermcell_flux(ngrid,klev,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)
     4SUBROUTINE 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)
    88     
    99     
    … …  
    2626     
    2727      INTEGER ngrid
    28       INTEGER klev
     28      INTEGER nlay
    2929      INTEGER igout
    3030      INTEGER lev_out
    … …  
    3434      INTEGER lalim(ngrid)
    3535     
    36       REAL alim_star(ngrid,klev)
    37       REAL entr_star(ngrid,klev)
    38       REAL detr_star(ngrid,klev)
    39       REAL zw2(ngrid,klev+1)
    40       REAL zlev(ngrid,klev+1)
    41       REAL masse(ngrid,klev)
     36      REAL alim_star(ngrid,nlay)
     37      REAL entr_star(ngrid,nlay)
     38      REAL detr_star(ngrid,nlay)
     39      REAL zw2(ngrid,nlay+1)
     40      REAL zlev(ngrid,nlay+1)
     41      REAL masse(ngrid,nlay)
    4242      REAL ptimestep
    43       REAL rhobarz(ngrid,klev)
     43      REAL rhobarz(ngrid,nlay)
    4444      REAL f(ngrid)
    45       REAL zqla(ngrid,klev)
     45      REAL zqla(ngrid,nlay)
    4646      REAL zmax(ngrid)
    4747     
    … …  
    4949!      --------     
    5050     
    51       REAL entr(ngrid,klev)
    52       REAL detr(ngrid,klev)
    53       REAL fm(ngrid,klev+1)
     51      REAL entr(ngrid,nlay)
     52      REAL detr(ngrid,nlay)
     53      REAL fm(ngrid,nlay+1)
    5454     
    5555!      local:
    … …  
    124124!------------------------------------------------------------------------------
    125125     
    126       DO l=1,klev
     126      DO l=1,nlay
    127127         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
    128128         detr(:,l) = f(:) * detr_star(:,l)
    … …  
    133133!------------------------------------------------------------------------------
    134134     
    135       DO l=1,klev
     135      DO l=1,nlay
    136136         DO ig=1,ngrid
    137137            IF (l.lt.lmax(ig) .and. l.ge.lmin(ig)) THEN
    … …  
    152152!==============================================================================
    153153     
    154       DO l=1,klev
     154      DO l=1,nlay
    155155         
    156156!------------------------------------------------------------------------------
    … …  
    186186! Test sur fraca constante ou croissante au-dessus de lalim
    187187!------------------------------------------------------------------------------
     188         
    188189! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
    189190         IF (iflag_thermals_optflux==0) THEN
    … …  
    205206! Test sur flux de masse constant ou decroissant au-dessus de lalim
    206207!------------------------------------------------------------------------------
     208         
    207209! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
    208210         IF (iflag_thermals_optflux==0) THEN
    … …  
    544546     
    545547!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    546 ! AB : temporary test added to check the equation df/dz = e - d validity
    547 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    548 !      DO l=1,klev
     548! AB : temporary test added to check the validity of eq. df/dz = e - d
     549!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     550!      DO l=1,nlay
    549551!         DO ig=1,ngrid
    550552!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))

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

    r2060 r2101  
    33!
    44      SUBROUTINE thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,          &
    5       &                           zw2,zlev,lmax,zmax,zmax0,zmix,              &
     5      &                           zw2,zlev,lmax,zmax,zmix,                    &
    66      &                           wmax,f_star,lev_out)
    77     
    … …  
    4040      REAL linter(ngrid)
    4141      REAL zmax(ngrid)
    42       REAL zmax0(ngrid)
    4342      REAL zmix(ngrid)
    4443      REAL wmax(ngrid)
    … …  
    138137         &             - zlev(ig,lmax(ig)))
    139138         zmax(ig) = max(zmax(ig), zlevinter(ig) - zlev(ig,lmin(ig)))
    140          zmax0(ig) = zmax(ig)
    141139      ENDDO
    142140     

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

    r2072 r2101  
    22!
    33!
    4       SUBROUTINE thermcell_main(itap,ngrid,nlay,ptimestep                     &
    5       &                  ,pplay,pplev,pphi,debut                              &
    6       &                  ,pu,pv,pt,po                                         &
    7       &                  ,pduadj,pdvadj,pdtadj,pdoadj                         &
    8       &                  ,f0,fm0,entr0,detr0                                  &
    9       &                  ,zqta,zqla,ztv,ztva,ztla,zthl,zqsatth                &
    10       &                  ,zmax0,zw2,fraca                                     &
    11       &                  ,lmin,lmix,lalim,lmax                                &
    12       &                  ,zpspsk,ratqscth,ratqsdiff                           &
    13       &                  ,Ale_bl,Alp_bl,lalim_conv,wght_th                    &
     4SUBROUTINE 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                    &
    1414!!! 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)
     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)
    2121!!! fin nrlmd le 10/04/2012
     22     
    2223     
    2324!==============================================================================
    … …  
    5152!==============================================================================
    5253     
    53       USE dimphy
    54 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    55 ! AB : why use dimphy to get klon, klev while subroutine arguments ngrid and
    56 !      nlay are the same?
    57 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    5854      USE thermcell_mod
    59       USE print_control_mod, ONLY: lunout,prt_level
    60       USE planete_mod, ONLY: preff
     55      USE print_control_mod, ONLY: lunout, prt_level
    6156     
    6257      IMPLICIT NONE
    … …  
    8479      REAL po(ngrid,nlay)                       !
    8580     
    86       LOGICAL debut
     81      LOGICAL firstcall
    8782     
    8883!   outputs:
    … …  
    9489      REAL pdoadj(ngrid,nlay)                   ! water convective variations
    9590     
    96       REAL fm0(klon,klev+1)                     ! mass flux      (after possible time relaxation)
    97       REAL entr0(klon,klev)                     ! entrainment    (after possible time relaxation)
    98       REAL detr0(klon,klev)                     ! detrainment    (after possible time relaxation)
    99       REAL f0(klon)                             ! mass flux norm (after possible time relaxation)
     91      REAL fm0(ngrid,nlay+1)                    ! mass flux      (after possible time relaxation)
     92      REAL entr0(ngrid,nlay)                    ! entrainment    (after possible time relaxation)
     93      REAL detr0(ngrid,nlay)                    ! detrainment    (after possible time relaxation)
     94      REAL f0(ngrid)                            ! mass flux norm (after possible time relaxation)
    10095     
    10196!   local:
    … …  
    118113     
    119114      INTEGER ig,k,l,ll,ierr
    120       INTEGER lmix_bis(klon)
    121       INTEGER lmax(klon)                        !
    122       INTEGER lmin(klon)                        !
    123       INTEGER lalim(klon)                       !
    124       INTEGER lmix(klon)                        !
    125      
    126       REAL linter(klon)
    127       REAL zmix(klon)
    128       REAL zmax(klon)
    129       REAL zw2(klon,klev+1)
    130       REAL zw_est(klon,klev+1)
    131       REAL zmax_sec(klon)
    132       REAL zmax0(klon)
    133      
    134       REAL zlay(klon,klev)                      ! layers altitude
    135       REAL zlev(klon,klev+1)                    ! levels altitude
    136       REAL rho(klon,klev)                       ! layers density
    137       REAL rhobarz(klon,klev)                   ! levels density
    138       REAL deltaz(klon,klev)                    ! layers heigth
    139       REAL masse(klon,klev)                     ! layers mass
    140       REAL zpspsk(klon,klev)                    ! Exner function
    141      
    142       REAL zu(klon,klev)                        ! environment zonal wind
    143       REAL zv(klon,klev)                        ! environment meridional wind
    144       REAL zo(klon,klev)                        ! environment water tracer
    145       REAL zva(klon,klev)                       ! plume zonal wind
    146       REAL zua(klon,klev)                       ! plume meridional wind
    147       REAL zoa(klon,klev)                       ! plume water tracer
    148      
    149       REAL zt(klon,klev)                        ! T    environment
    150       REAL zh(klon,klev)                        ! T,TP environment
    151       REAL zthl(klon,klev)                      ! TP   environment
    152       REAL ztv(klon,klev)                       ! TPV  environment ?
    153       REAL zl(klon,klev)                        ! ql   environment
    154      
    155       REAL zta(klon,klev)                       !
    156       REAL zha(klon,klev)                       ! TRPV plume
    157       REAL ztla(klon,klev)                      ! TP   plume
    158       REAL ztva(klon,klev)                      ! TRPV plume
    159       REAL ztva_est(klon,klev)                  ! TRPV plume (temporary)
    160       REAL zqla(klon,klev)                      ! qv   plume
    161       REAL zqta(klon,klev)                      ! qt   plume
    162      
    163       REAL wmax(klon)                           ! maximal vertical speed
    164       REAL wmax_tmp(klon)                       ! temporary maximal vertical speed
    165       REAL wmax_sec(klon)                       ! maximal vertical speed if dry case
    166      
    167       REAL fraca(klon,klev+1)                   ! updraft fraction
    168       REAL f_star(klon,klev+1)                  ! normalized mass flux
    169       REAL entr_star(klon,klev)                 ! normalized entrainment
    170       REAL detr_star(klon,klev)                 ! normalized detrainment
    171       REAL alim_star_tot(klon)                  ! integrated alimentation
    172       REAL alim_star(klon,klev)                 ! normalized alimentation
    173       REAL alim_star_clos(klon,klev)            ! closure alimentation
    174      
    175       REAL fm(klon,klev+1)                      ! mass flux
    176       REAL entr(klon,klev)                      ! entrainment
    177       REAL detr(klon,klev)                      ! detrainment
    178       REAL f(klon)                              ! mass flux norm
    179      
    180       REAL zdthladj(klon,klev)                  !
     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)
     154      REAL zqla(ngrid,nlay)                     ! qv   plume
     155      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)                 !
    181175      REAL lambda                               ! time relaxation coefficent
    182176     
    183       REAL zsortie(klon,klev)
    184       REAL zsortie1d(klon)
     177      REAL zsortie(ngrid,nlay)
     178      REAL zsortie1d(ngrid)
    185179      REAL susqr2pi, Reuler
    186180      REAL zf
    187181      REAL zf2
    188       REAL thetath2(klon,klev)
    189       REAL wth2(klon,klev)
    190       REAL wth3(klon,klev)
    191       REAL q2(klon,klev)
     182      REAL thetath2(ngrid,nlay)
     183      REAL wth2(ngrid,nlay)
     184      REAL wth3(ngrid,nlay)
     185      REAL q2(ngrid,nlay)
    192186! FH : probleme de dimensionnement avec l'allocation dynamique
    193187!     common/comtherm/thetath2,wth2
    194       real wq(klon,klev)
    195       real wthl(klon,klev)
    196       real wthv(klon,klev)
    197       real ratqscth(klon,klev)
     188      real wq(ngrid,nlay)
     189      real wthl(ngrid,nlay)
     190      real wthv(ngrid,nlay)
     191      real ratqscth(ngrid,nlay)
    198192      real var
    199193      real vardiff
    200       real ratqsdiff(klon,klev)
     194      real ratqsdiff(ngrid,nlay)
    201195! niveau de condensation
    202       integer nivcon(klon)
    203       real zcon(klon)
     196      integer nivcon(ngrid)
     197      real zcon(ngrid)
    204198      real CHI
    205       real zcon2(klon)
    206       real pcon(klon)
    207       real zqsat(klon,klev)
    208       real zqsatth(klon,klev)
    209 ! FH/IM :   save f0
    210       real zlevinter(klon)
     199      real zcon2(ngrid)
     200      real pcon(ngrid)
     201      real zqsat(ngrid,nlay)
     202      real zqsatth(ngrid,nlay)
     203      real zlevinter(ngrid)
    211204      real seuil
    212205     
    … …  
    214207     
    215208!------Entrees
    216       real pbl_tke(klon,klev+1,nbsrf)
    217       real pctsrf(klon,nbsrf)
    218       real omega(klon,klev)
    219       real airephy(klon)
     209      real pbl_tke(ngrid,nlay+1,nbsrf)
     210      real pctsrf(ngrid,nbsrf)
     211      real omega(ngrid,nlay)
     212      real airephy(ngrid)
    220213!------Sorties
    221       real zlcl(klon),fraca0(klon),w0(klon),w_conv(klon)
    222       real therm_tke_max0(klon),env_tke_max0(klon)
    223       real n2(klon),s2(klon)
    224       real ale_bl_stat(klon)
    225       real therm_tke_max(klon,klev),env_tke_max(klon,klev)
    226       real alp_bl_det(klon),alp_bl_fluct_m(klon),alp_bl_fluct_tke(klon),alp_bl_conv(klon),alp_bl_stat(klon)
     214      real zlcl(ngrid),fraca0(ngrid),w0(ngrid),w_conv(ngrid)
     215      real therm_tke_max0(ngrid),env_tke_max0(ngrid)
     216      real n2(ngrid),s2(ngrid)
     217      real ale_bl_stat(ngrid)
     218      real therm_tke_max(ngrid,nlay),env_tke_max(ngrid,nlay)
     219      real alp_bl_det(ngrid),alp_bl_fluct_m(ngrid),alp_bl_fluct_tke(ngrid),alp_bl_conv(ngrid),alp_bl_stat(ngrid)
    227220!------Local
    228221      integer nsrf
    229       real rhobarz0(klon)                    ! Densite au LCL
    230       logical ok_lcl(klon)                   ! Existence du LCL des thermiques
    231       integer klcl(klon)                     ! Niveau du LCL
    232       real interp(klon)                      ! Coef d'interpolation pour le LCL
     222      real rhobarz0(ngrid)                      ! Densite au LCL
     223      logical ok_lcl(ngrid)                     ! Existence du LCL des thermiques
     224      integer klcl(ngrid)                       ! Niveau du LCL
     225      real interp(ngrid)                        ! Coef d'interpolation pour le LCL
    233226!------Triggering
    234       real Su                                ! Surface unite: celle d'un updraft elementaire
     227      real Su                                   ! Surface unite: celle d'un updraft elementaire
    235228      parameter(Su=4e4)
    236       real hcoef                             ! Coefficient directeur pour le calcul de s2
     229      real hcoef                                ! Coefficient directeur pour le calcul de s2
    237230      parameter(hcoef=1)
    238       real hmincoef                          ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2
     231      real hmincoef                             ! Coefficient directeur pour l'ordonnee a l'origine pour le calcul de s2
    239232      parameter(hmincoef=0.3)
    240       real eps1                              ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
     233      real eps1                                 ! Fraction de surface occupee par la population 1 : eps1=n1*s1/(fraca0*Sd)
    241234      parameter(eps1=0.3)
    242       real hmin(ngrid)                       ! Ordonnee a l'origine pour le calcul de s2
    243       real zmax_moy(ngrid)                   ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
     235      real hmin(ngrid)                          ! Ordonnee a l'origine pour le calcul de s2
     236      real zmax_moy(ngrid)                      ! Hauteur moyenne des thermiques : zmax_moy = zlcl + 0.33 (zmax-zlcl)
    244237      real zmax_moy_coef
    245238      parameter(zmax_moy_coef=0.33)
    246       real depth(klon)                       ! Epaisseur moyenne du cumulus
    247       real w_max(klon)                       ! Vitesse max statistique
    248       real s_max(klon)
     239      real depth(ngrid)                         ! Epaisseur moyenne du cumulus
     240      real w_max(ngrid)                         ! Vitesse max statistique
     241      real s_max(ngrid)
    249242!------Closure
    250       real pbl_tke_max(klon,klev)            ! Profil de TKE moyenne
    251       real pbl_tke_max0(klon)                ! TKE moyenne au LCL
    252       real w_ls(klon,klev)                   ! Vitesse verticale grande echelle (m/s)
    253       real coef_m                            ! On considere un rendement pour alp_bl_fluct_m
     243      real pbl_tke_max(ngrid,nlay)              ! Profil de TKE moyenne
     244      real pbl_tke_max0(ngrid)                  ! TKE moyenne au LCL
     245      real w_ls(ngrid,nlay)                     ! Vitesse verticale grande echelle (m/s)
     246      real coef_m                               ! On considere un rendement pour alp_bl_fluct_m
    254247      parameter(coef_m=1.)
    255       real coef_tke                          ! On considere un rendement pour alp_bl_fluct_tke
     248      real coef_tke                             ! On considere un rendement pour alp_bl_fluct_tke
    256249      parameter(coef_tke=1.)
    257250     
    … …  
    259252     
    260253! Nouvelles variables pour la convection
    261       real Ale_bl(klon)
    262       real Alp_bl(klon)
    263       real alp_int(klon),dp_int(klon),zdp
    264       real ale_int(klon)
    265       integer n_int(klon)
    266       real fm_tot(klon)
    267       real wght_th(klon,klev)
    268       integer lalim_conv(klon)
     254      real Ale_bl(ngrid)
     255      real Alp_bl(ngrid)
     256      real alp_int(ngrid),dp_int(ngrid),zdp
     257      real ale_int(ngrid)
     258      integer n_int(ngrid)
     259      real fm_tot(ngrid)
     260      real wght_th(ngrid,nlay)
     261      integer lalim_conv(ngrid)
    269262     
    270263      CHARACTER*2 str2
    … …  
    282275      seuil = 0.25
    283276     
    284       IF (debut) THEN
     277      IF (firstcall) THEN
    285278         IF (iflag_thermals==15.or.iflag_thermals==16) THEN
    286279            dvimpl = 0
    … …  
    299292      entr(:,:) = 0.
    300293      detr(:,:) = 0.
    301      
    302       IF (ngrid.ne.klon) THEN
    303          print *, 'ERROR: ngrid and klon are different!'
    304          print *, 'ngrid =', ngrid
    305          print *, 'klon  =', klon
    306       ENDIF
    307      
    308       DO ig=1,klon
     294      f(:) = 0.
     295     
     296      DO ig=1,ngrid
    309297         f0(ig) = max(f0(ig), 1.e-2)
    310          zmax0(ig) = max(zmax0(ig),40.)
    311298      ENDDO
    312299     
    … …  
    352339     
    353340      zlev(:,1) = 0.
    354       zlev(:,nlay+1) = (2. * pphi(:,klev) - pphi(:,klev-1)) / RG
     341      zlev(:,nlay+1) = (2. * pphi(:,nlay) - pphi(:,nlay-1)) / RG
    355342     
    356343      DO l=1,nlay
    … …  
    467454!==============================================================================
    468455     
    469       CALL thermcell_plume(itap,ngrid,nlay,ptimestep,ztv,zthl,             &
    470       &    po,zl,rhobarz,zlev,pplev,pphi,zpspsk,alim_star,alim_star_tot,   &
    471       &    lalim,f0,detr_star,entr_star,f_star,ztva,                       &
    472       &    ztla,zqla,zqta,zha,zw2,zw_est,ztva_est,zqsatth,lmix,lmix_bis,   &
     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,      &
    473460      &    lmin,lev_out,lunout1,igout)
    474461     
    475       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lalim ')
    476       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_plum lmix  ')
    477      
    478 !==============================================================================
    479 ! Calcul des caracteristiques du thermique:zmax,zmix,wmax
    480 !==============================================================================
    481      
    482       CALL thermcell_height(ngrid,nlay,lalim,lmin,linter,lmix,zw2,   &
    483       &                              zlev,lmax,zmax,zmax0,zmix,wmax,f_star,   &
    484       &                              lev_out)
     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!==============================================================================
     466! 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)
    485471     
    486472!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    488474!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    489475     
    490       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lalim ')
    491       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmin  ')
    492       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmix ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmix  ')
    493       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_heig lmax  ')
    494      
    495 !==============================================================================
    496 ! Closure and mass flux determining
     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!==============================================================================
     482! Closure and mass fluxes computation
    497483!==============================================================================
    498484     
    … …  
    500486      &                  lalim,lmin,zmax_sec,wmax_sec,lev_out)
    501487     
    502       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmin ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lmin  ')
    503       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_dry  lalim ')
     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 ')
    504490     
    505491!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    509495      alim_star_clos(:,:) = alim_star(:,:)
    510496      alim_star_clos(:,:) = entr_star(:,:) + alim_star(:,:)
    511          
     497     
    512498!------------------------------------------------------------------------------
    513499! Closure (dry if iflag_thermals_closure=1, moist if iflag_thermals_closure=2)
    … …  
    531517         f0 = (1.-lambda) * f + lambda * f0
    532518      ELSE
    533          f0 = f
     519         f0(:) = f(:)
    534520      ENDIF
    535521     
    … …  
    539525      IF (.not. (f0(1).ge.0.) ) THEN
    540526         abort_message = '.not. (f0(1).ge.0.)'
    541          write(*,*) 'f0 =', f0(1)
     527         print *, 'f0 =', f0(1)
    542528         CALL abort_physic(modname,abort_message,1)
     529      ELSE
     530         print *, 'f0 =', f0(1)
    543531      ENDIF
    544532     
    545 !==============================================================================
    546 ! Deduction des flux
    547 !==============================================================================
     533!------------------------------------------------------------------------------
     534! Mass fluxes
     535!------------------------------------------------------------------------------
    548536     
    549537      CALL thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
    550       &       lalim,lmin,lmax,alim_star,entr_star,detr_star,                  &
    551       &       f,rhobarz,zlev,zw2,fm,entr,                                     &
    552       &       detr,zqla,lev_out,lunout1,igout)
    553      
    554       CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
    555       CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
     538      &                   lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
     539      &                   f,rhobarz,zlev,zw2,fm,entr,detr,zqla,               &
     540      &                   lev_out,lunout1,igout)
     541     
     542!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lalim,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lalim ')
     543!      CALL test_ltherm(ngrid,nlay,pplev,pplay,lmax ,seuil,ztv,po,ztva,zqla,f_star,zw2,'thermcell_flux lmax  ')
    556544     
    557545!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    567555         detr0 = (1.-lambda) * detr + lambda * detr0
    568556      ELSE
    569          fm0   = fm
    570          entr0 = entr
    571          detr0 = detr
     557         fm0(:,:)   = fm(:,:)
     558         entr0(:,:) = entr(:,:)
     559         detr0(:,:) = detr(:,:)
    572560      ENDIF
    573561     
    574 !==============================================================================
    575 ! Calcul du transport vertical (de qt et tp)
    576 !==============================================================================
    577      
    578       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
    579       &                 zthl,zdthladj,zta,lmin,lev_out)
    580      
    581       CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
    582       &                 po,pdoadj,zoa,lmin,lev_out)
    583      
    584       DO l=1,nlay
     562!------------------------------------------------------------------------------
     563! Updraft fraction
     564!------------------------------------------------------------------------------
     565     
     566      DO ig=1,ngrid
     567         fraca(ig,1) = 0.
     568         fraca(ig,nlay+1) = 0.
     569      ENDDO
     570     
     571      DO l=2,nlay
    585572         DO ig=1,ngrid
    586            pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l) 
    587          ENDDO
    588       ENDDO
    589      
    590 !==============================================================================
    591 ! Calcul de la fraction de l'ascendance
    592 !==============================================================================
    593      
    594       DO ig=1,klon
    595          fraca(ig,1)=0.
    596          fraca(ig,nlay+1)=0.
    597       ENDDO
    598      
    599       DO l=2,nlay
    600          DO ig=1,klon
    601             IF (zw2(ig,l).gt.1.e-10) THEN
     573            IF (zw2(ig,l).gt.0.) THEN
    602574               fraca(ig,l) = fm(ig,l) / (rhobarz(ig,l) * zw2(ig,l))
    603575            ELSE
    … …  
    608580     
    609581!==============================================================================
     582! Transport vertical
     583!==============================================================================
     584     
     585!------------------------------------------------------------------------------
     586! Calcul du transport vertical (de qt et tp)
     587!------------------------------------------------------------------------------
     588     
     589      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
     590      &                 zthl,zdthladj,zta,lmin,lev_out)
     591     
     592      CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,          &
     593      &                 po,pdoadj,zoa,lmin,lev_out)
     594     
     595      DO l=1,nlay
     596         DO ig=1,ngrid
     597           pdtadj(ig,l) = zdthladj(ig,l) * zpspsk(ig,l) 
     598         ENDDO
     599      ENDDO
     600     
     601!------------------------------------------------------------------------------
    610602! Calcul du transport vertical du moment horizontal
    611 !==============================================================================
     603!------------------------------------------------------------------------------
    612604     
    613605      IF (dvimpl.eq.0) THEN
    … …  
    620612      ELSE
    621613!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    622 ! Calcul purement conservatif pour le transport de V
     614! Calcul purement conservatif pour le transport de V=(u,v)
    623615!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    624616         CALL thermcell_dq(ngrid,nlay,dqimpl,ptimestep,fm0,entr0,masse,       &
    … …  
    756748         ratqsdiff(:,:) = 0.
    757749         
    758          DO l=1,klev
     750         DO l=1,nlay
    759751            DO ig=1,ngrid
    760752               IF (l<=lalim(ig)) THEN
    … …  
    766758         if (prt_level.ge.1) print*,'14f OK convect8'
    767759     
    768          DO l=1,klev
     760         DO l=1,nlay
    769761            DO ig=1,ngrid
    770762               IF (l<=lalim(ig)) THEN
    … …  
    799791
    800792
    801       SUBROUTINE test_ltherm(klon,klev,pplev,pplay,long,seuil,ztv,po,         &
    802                              ztva,zqla,f_star,zw2,comment)
     793SUBROUTINE test_ltherm(ngrid,nlay,pplev,pplay,long,seuil,ztv,po,              &
     794                       ztva,zqla,f_star,zw2,comment)
    803795     
    804796     
    … …  
    815807!      -------
    816808     
    817       INTEGER klon
    818       INTEGER klev
    819       INTEGER long(klon)
    820      
    821       REAL pplev(klon,klev+1)
    822       REAL pplay(klon,klev)
    823       REAL ztv(klon,klev)
    824       REAL po(klon,klev)
    825       REAL ztva(klon,klev)
    826       REAL zqla(klon,klev)
    827       REAL f_star(klon,klev)
    828       REAL zw2(klon,klev)
     809      INTEGER ngrid
     810      INTEGER nlay
     811      INTEGER long(ngrid)
     812     
     813      REAL pplev(ngrid,nlay+1)
     814      REAL pplay(ngrid,nlay)
     815      REAL ztv(ngrid,nlay)
     816      REAL po(ngrid,nlay)
     817      REAL ztva(ngrid,nlay)
     818      REAL zqla(ngrid,nlay)
     819      REAL f_star(ngrid,nlay)
     820      REAL zw2(ngrid,nlay)
    829821      REAL seuil
    830822     
    … …  
    847839     
    848840!  test sur la hauteur des thermiques ...
    849       do i=1,klon
     841      do i=1,ngrid
    850842!IMtemp if (pplay(i,long(i)).lt.seuil*pplev(i,1)) then
    851843        if (prt_level.ge.10) then
    852844            print *, 'WARNING ',comment,' au point ',i,' K= ',long(i)
    853845            print *, '  K  P(MB)  THV(K)     Qenv(g/kg)THVA        QLA(g/kg)   F*        W2'
    854             do k=1,klev
     846            do k=1,nlay
    855847               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)
    856848            enddo
    … …  
    871863
    872864
    873       SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,      &
    874       &                                  rg,pplev,therm_tke_max)
    875      
    876      
    877       USE print_control_mod, ONLY: prt_level
    878      
    879       IMPLICIT NONE
     865SUBROUTINE thermcell_tke_transport(ngrid,nlay,ptimestep,fm0,entr0,            &
     866                                   rg,pplev,therm_tke_max)
    880867     
    881868     
    … …  
    887874!   Transport de la TKE par le thermique moyen pour la fermeture en ALP
    888875!   On transporte pbl_tke pour donner therm_tke
     876!==============================================================================
     877     
     878      USE print_control_mod, ONLY: prt_level
     879     
     880      IMPLICIT NONE
     881     
     882     
     883!==============================================================================
     884! Declarations
    889885!==============================================================================
    890886     
    … …  
    900896      real entr(ngrid,nlay)
    901897      real q(ngrid,nlay)
    902       integer lev_out                           ! niveau pour les print
    903898     
    904899      real qa(ngrid,nlay),detr(ngrid,nlay),wqd(ngrid,nlay+1)
    … …  
    909904      integer isrf
    910905     
    911       lev_out=0
    912      
    913       if (prt_level.ge.1) print*,'Q2 THERMCEL_DQ 0'
    914      
    915 !   calcul du detrainement
     906!------------------------------------------------------------------------------
     907! Calcul du detrainement
     908!------------------------------------------------------------------------------
     909     
    916910      do k=1,nlay
    917          detr0(:,k)=fm0(:,k)-fm0(:,k+1)+entr0(:,k)
    918          masse0(:,k)=(pplev(:,k)-pplev(:,k+1))/RG
     911         detr0(:,k) = fm0(:,k) - fm0(:,k+1) + entr0(:,k)
     912         masse0(:,k) = (pplev(:,k) - pplev(:,k+1)) / RG
    919913      enddo
    920 
    921 
     914     
     915!------------------------------------------------------------------------------
    922916! Decalage vertical des entrainements et detrainements.
     917!------------------------------------------------------------------------------
     918     
    923919      masse(:,1)=0.5*masse0(:,1)
    924920      entr(:,1)=0.5*entr0(:,1)
    925921      detr(:,1)=0.5*detr0(:,1)
    926922      fm(:,1)=0.
     923     
    927924      do k=1,nlay-1
    928925         masse(:,k+1)=0.5*(masse0(:,k)+masse0(:,k+1))
    … …  
    931928         fm(:,k+1)=fm(:,k)+entr(:,k)-detr(:,k)
    932929      enddo
     930     
    933931      fm(:,nlay+1)=0.
    934 
     932     
    935933!!! nrlmd le 16/09/2010
    936934!   calcul de la valeur dans les ascendances
    937 !       do ig=1,ngrid
    938 !          qa(ig,1)=q(ig,1)
    939 !       enddo
    940 !!!
    941 
    942 !do isrf=1,nsrf
    943 
    944 !   q(:,:)=therm_tke(:,:,isrf)
    945    q(:,:)=therm_tke_max(:,:)
    946 !!! nrlmd le 16/09/2010
     935!      do ig=1,ngrid
     936!         qa(ig,1) = q(ig,1)
     937!      enddo
     938     
     939      q(:,:)=therm_tke_max(:,:)
     940     
    947941      do ig=1,ngrid
    948942         qa(ig,1)=q(ig,1)
    949943      enddo
    950 !!!
    951 
    952     if (1==1) then
    953       do k=2,nlay
     944     
     945      if (1==1) then
     946         do k=2,nlay
     947            do ig=1,ngrid
     948               if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.1.e-5*masse(ig,k)) then
     949                  qa(ig,k) = (fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))         &
     950                  &        / (fm(ig,k+1)+detr(ig,k))
     951               else
     952                  qa(ig,k)=q(ig,k)
     953               endif
     954               
     955!               if (qa(ig,k).lt.0.) print *, 'WARNING: qa is negative!'
     956!               if (q(ig,k).lt.0.) print *, 'WARNING: q is negative!'
     957            enddo
     958         enddo
     959         
     960!------------------------------------------------------------------------------
     961! Calcul du flux subsident
     962!------------------------------------------------------------------------------
     963         
     964         do k=2,nlay
     965            do ig=1,ngrid
     966               wqd(ig,k) = fm(ig,k) * q(ig,k)
     967               if (wqd(ig,k).lt.0.) print*,'WARNING: wqd is negative!'
     968            enddo
     969         enddo
     970         
    954971         do ig=1,ngrid
    955             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
    956      &         1.e-5*masse(ig,k)) then
    957          qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
    958      &         /(fm(ig,k+1)+detr(ig,k))
    959             else
    960                qa(ig,k)=q(ig,k)
    961             endif
    962             if (qa(ig,k).lt.0.) then
    963 !               print*,'qa<0!!!'
    964             endif
    965             if (q(ig,k).lt.0.) then
    966 !               print*,'q<0!!!'
    967             endif
    968          enddo
    969       enddo
    970 
    971 ! Calcul du flux subsident
    972 
    973       do k=2,nlay
    974          do ig=1,ngrid
    975             wqd(ig,k)=fm(ig,k)*q(ig,k)
    976             if (wqd(ig,k).lt.0.) then
    977 !               print*,'wqd<0!!!'
    978             endif
    979          enddo
    980       enddo
    981       do ig=1,ngrid
    982          wqd(ig,1)=0.
    983          wqd(ig,nlay+1)=0.
    984       enddo
    985 
     972            wqd(ig,1) = 0.
     973            wqd(ig,nlay+1) = 0.
     974         enddo
     975         
     976!------------------------------------------------------------------------------
    986977! Calcul des tendances
    987       do k=1,nlay
    988          do ig=1,ngrid
    989             q(ig,k)=q(ig,k)+(detr(ig,k)*qa(ig,k)-entr(ig,k)*q(ig,k)  &
    990      &               -wqd(ig,k)+wqd(ig,k+1))  &
    991      &               *ptimestep/masse(ig,k)
    992          enddo
    993       enddo
    994 
    995  endif
    996 
    997    therm_tke_max(:,:)=q(:,:)
     978!------------------------------------------------------------------------------
     979         
     980         do k=1,nlay
     981            do ig=1,ngrid
     982               q(ig,k) = q(ig,k) + ptimestep / masse(ig,k)                    &
     983               &       * (detr(ig,k) * qa(ig,k) - entr(ig,k) * q(ig,k)        &
     984               &       - wqd(ig,k) + wqd(ig,k+1))
     985            enddo
     986         enddo
     987      endif
     988     
     989      therm_tke_max(:,:) = q(:,:)
    998990     
    999991     

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

    r2093 r2101  
    22!
    33!
    4       SUBROUTINE thermcell_plume(itap,ngrid,klev,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)
     4SUBROUTINE 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)
    1010     
    1111     
    … …  
    3636      INTEGER itap
    3737      INTEGER ngrid
    38       INTEGER klev
     38      INTEGER nlay
    3939      INTEGER lunout1
    4040      INTEGER igout
    … …  
    4242     
    4343      REAL ptimestep                            ! time step
    44       REAL ztv(ngrid,klev)                      ! TRPV environment
    45       REAL zthl(ngrid,klev)                     ! TP   environment
    46       REAL po(ngrid,klev)                       ! qt   environment
    47       REAL zl(ngrid,klev)                       ! ql   environment
    48       REAL rhobarz(ngrid,klev)                  ! levels density
    49       REAL zlev(ngrid,klev+1)                   ! levels altitude
    50       REAL pplev(ngrid,klev+1)                  ! levels pressure
    51       REAL pphi(ngrid,klev)                     ! geopotential
    52       REAL zpspsk(ngrid,klev)                   ! Exner function
     44      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
    5353     
    5454!      outputs:
    … …  
    6060      INTEGER lmix_bis(ngrid)                   ! maximum vertical speed level (modified)
    6161     
    62       REAL alim_star(ngrid,klev)                ! normalized alimentation
     62      REAL alim_star(ngrid,nlay)                ! normalized alimentation
    6363      REAL alim_star_tot(ngrid)                 ! integrated alimentation
    6464     
    6565      REAL f0(ngrid)                            ! previous time step mass flux norm
    6666     
    67       REAL detr_star(ngrid,klev)                ! normalized detrainment
    68       REAL entr_star(ngrid,klev)                ! normalized entrainment
    69       REAL f_star(ngrid,klev+1)                 ! normalized mass flux
    70      
    71       REAL ztva(ngrid,klev)                     ! TRPV plume (after mixing)
    72       REAL ztva_est(ngrid,klev)                 ! TRPV plume (before mixing)
    73       REAL ztla(ngrid,klev)                     ! TP   plume
    74       REAL zqla(ngrid,klev)                     ! ql   plume (after mixing)
    75       REAL zqta(ngrid,klev)                     ! qt   plume
    76       REAL zha(ngrid,klev)                      ! TRPV plume
    77      
    78       REAL w_est(ngrid,klev+1)                  ! updraft square vertical speed (before mixing)
    79       REAL zw2(ngrid,klev+1)                    ! updraft square vertical speed (after mixing)
    80      
    81       REAL zqsatth(ngrid,klev)                  ! saturation vapor pressure (after mixing)
     67      REAL detr_star(ngrid,nlay)                ! normalized detrainment
     68      REAL entr_star(ngrid,nlay)                ! normalized entrainment
     69      REAL f_star(ngrid,nlay+1)                 ! normalized mass flux
     70     
     71      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
     74      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)
    8282     
    8383!      local:
    … …  
    8888      INTEGER lm
    8989     
    90       REAL zqla_est(ngrid,klev)                 ! ql   plume (before mixing)
    91       REAL zta_est(ngrid,klev)                  ! TR   plume (before mixing)
    92       REAL zbuoy(ngrid,klev)                    ! plume buoyancy
    93       REAL zbuoyjam(ngrid,klev)                 ! plume buoyancy (modified)
     90      REAL zqla_est(ngrid,nlay)                 ! ql   plume (before mixing)
     91      REAL zta_est(ngrid,nlay)                  ! TR   plume (before mixing)
     92      REAL zbuoy(ngrid,nlay)                    ! plume buoyancy
     93      REAL zbuoyjam(ngrid,nlay)                 ! plume buoyancy (modified)
    9494     
    9595      REAL ztemp(ngrid)                         ! temperature for saturation vapor pressure computation in plume
    9696      REAL zqsat(ngrid)                         ! saturation vapor pressure (before mixing)
    9797      REAL zdz                                  ! layers height
    98       REAL ztv2(ngrid,klev)                     ! ztv + d_temp * Dirac(l=lmin)
     98      REAL ztv2(ngrid,nlay)                     ! ztv + d_temp * Dirac(l=lmin)
    9999     
    100100      REAL zalpha                               !
    101       REAL zdqt(ngrid,klev)                     !
     101      REAL zdqt(ngrid,nlay)                     !
    102102      REAL zbetalpha                            !
    103103      REAL zdw2                                 !
    … …  
    118118      REAL zltup                                ! useless here
    119119     
    120       REAL dummy
     120      REAL psat                                 ! dummy argument for Psat_water()
    121121     
    122122      LOGICAL active(ngrid)                     ! if the plume is active at ig,l (speed and incoming mass flux > 0 or l=lmin)
    … …  
    161161     
    162162      ztv2(:,:)        = ztv(:,:)
    163       ztv2(:,linf)     = ztv2(:,linf) + d_temp
     163      ztv2(:,linf)     = ztv(:,linf) + d_temp
    164164     
    165165!==============================================================================
    … …  
    178178         active(ig) = .false.
    179179         l = linf
    180          DO WHILE ((.not.active(ig)) .and. pplev(ig,l+1).gt.pres_limit .and. l.lt.klev)
     180         DO WHILE ((.not.active(ig)) .and. pplev(ig,l+1).gt.pres_limit .and. l.lt.nlay)
    181181            IF (ztv2(ig,l).gt.ztv2(ig,l+1)) THEN
    182182               active(ig) = .true.
    … …  
    190190! AB : On pourrait n'appeler thermcell_alim que si la plume est active
    191191!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    192       CALL thermcell_alim(ngrid,klev,ztv2,zlev,alim_star,lalim,lmin)
     192      CALL thermcell_alim(ngrid,nlay,ztv2,zlev,alim_star,lalim,lmin)
    193193     
    194194!==============================================================================
    … …  
    222222!==============================================================================
    223223     
    224       DO l=2,klev-1
     224      DO l=2,nlay-1
    225225         
    226226!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    242242         
    243243         DO ig=1,ngrid
    244             CALL Psat_water(ztemp(ig), pplev(ig,l), dummy, zqsat(ig))
     244            CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsat(ig))
    245245         ENDDO
    246246         
    … …  
    434434         DO ig=1,ngrid
    435435            IF (activetmp(ig)) THEN
    436                CALL Psat_water(ztemp(ig), pplev(ig,l), dummy, zqsatth(ig,l))
     436               CALL Psat_water(ztemp(ig), pplev(ig,l), psat, zqsatth(ig,l))
    437437            ENDIF
    438438         ENDDO
Note: See TracChangeset for help on using the changeset viewer.