Changeset 1032 for trunk


Ignore:
Timestamp:
Sep 6, 2013, 8:51:56 AM (12 years ago)
Author:
aslmd
Message:

LMDZ.MARS cleaned and commented version of the thermal plume model with automatic arrays.

Location:
trunk/LMDZ.MARS
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r1013 r1032  
    18881888- Upgrade on the thermospheric photochemical reaction rates. These are now
    18891889  read in from a file "chemthermos_reactionrates.def".
     1890
     1891== 06/09/2013 == AC + EM + AS
     1892- Cleaned and commented version of thermal plume model with automatic arrays
     1893- Checked: exact same results than before modifications

  • trunk/LMDZ.MARS/libf/phymars/calltherm_interface.F90

    r1028 r1032  
    1 ! ---------------------------------------------------------------------
    2 ! Arnaud Colaitis 2011-01-05
    3 !
    4 ! Interface routine between physiq.F (driver for physics) and thermcell_main_mars (thermals model)
    5 ! Handles sub-timestep for the thermals model, outputs and diagnostics
    6 !
    7 ! The thermals model implemented in the Mars LMD-GCM is called after the vertical turbulent mixing scheme (Mellor and Yamada) 
    8 ! and surface layer scheme (Richardson-based surface layer with subgrid gustiness parameterization)
    9 ! Other routines called before the thermals model are Radiative transfer scheme (Madeleine et.al) 
    10 ! and gravity wave and subgrid scale topography drag.
    11 !
    12 ! Reference paper for the Martian thermals model is Colaitis et al. 2013 (JGR-planets)
    13 ! Mentions in the routines to a "paper" points to the above reference.
    14 !
    15 ! This model is based on the LMDZ thermals model from C. Rio and F. Hourdin
    16 ! ---------------------------------------------------------------------
    17       SUBROUTINE calltherm_interface (firstcall, &
     1! -----------------------------------------------------------------------
     2! CALLTHERM_INTERFACE
     3! -----------------------------------------------------------------------
     4! Main interface to the Martian thermal plume model
     5! This interface handles sub-timesteps for this model
     6! A call to this interface must be inserted in the main 'physics' routine
     7!   NB: for information:
     8!   In the Mars LMD-GCM, the thermal plume model is called after
     9!   the vertical turbulent mixing scheme (Mellor and Yamada) 
     10!   and the surface layer scheme (Richardson-based surface layer + subgrid gustiness)
     11!   Other routines called before the thermals model are
     12!   radiative transfer and (orographic) gravity wave drag.
     13! -----------------------------------------------------------------------
     14! ASSOCIATED FILES
     15! --> thermcell_dqup.F90
     16! --> thermcell_main_mars.F90
     17! --> comtherm_h.F90
     18! -----------------------------------------------------------------------
     19! Reference paper:
     20! A. Colaïtis, A. Spiga, F. Hourdin, C. Rio, F. Forget, and E. Millour.
     21! A thermal plume model for the Martian convective boundary layer.
     22! Journal of Geophysical Research (Planets), 118:1468-1487, July 2013.
     23! -----------------------------------------------------------------------
     24! Reference paper for terrestrial plume model:
     25! C. Rio and F. Hourdin.
     26! A thermal plume model for the convective boundary layer : Representation of cumulus clouds.
     27! Journal of the Atmospheric Sciences, 65:407-425, 2008.
     28! -----------------------------------------------------------------------
     29! Author : A. Colaitis 2011-01-05 (with updates 2011-2013)
     30! Institution : Laboratoire de Meteorologie Dynamique (LMD) Paris, France
     31! Corresponding author : A. Spiga aymeric.spiga_AT_upmc.fr
     32! -----------------------------------------------------------------------
     33
     34      SUBROUTINE calltherm_interface (ngrid,nlayer,nq, &
     35     & tracer,igcm_co2, &
    1836     & zzlev,zzlay, &
    1937     & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
    … …  
    2240     & pdhdif,hfmax,wstar,sensibFlux)
    2341
    24 
    25       USE ioipsl_getincom, only : getin ! to read options from external file "run.def"
     42      use comtherm_h
    2643      implicit none
    27 #include "callkeys.h" !contains logical values determining which subroutines and which options are used.
    28                       ! includes logical "tracer", which is .true. if tracers are present and to be transported
    29                       ! includes integer "iradia", frequency of calls to radiative scheme wrt calls to physics (and is typically==1)
    30 #include "dimensions.h" !contains global GCM grid dimension informations
    31 #include "dimphys.h" !similar to dimensions.h, and based on it; includes
    32                      ! "ngridmx": number of horizontal grid points
    33                      ! "nlayermx": number of atmospheric layers
     44
     45! SHARED VARIABLES. This needs adaptations in another climate model.
    3446#include "comcstfi.h" !contains physical constant values such as
    3547                      ! "g" : gravitational acceleration (m.s-2)
    3648                      ! "r" : recuced gas constant (J.K-1.mol-1)
    3749                      ! "cpp" : specific heat of the atmosphere (J.kg-1.K-1)
    38 #include "tracer.h" !contains tracer information such as
    39                     ! "nqmx": number of tracers
    40                     ! "noms()": tracer name
    4150
    4251!--------------------------------------------------------
    … …  
    4453!--------------------------------------------------------
    4554
     55      INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
     56      INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
     57      INTEGER, INTENT(IN) :: nq ! number of tracer species
    4658      REAL, INTENT(IN) :: ptimestep !timestep (s)
    47       REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1) !intermediate pressure levels (Pa)
    48       REAL, INTENT(IN) :: pplay(ngridmx,nlayermx) !Pressure at the middle of the layers (Pa)
    49       REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) !Geopotential at the middle of the layers (m2s-2)
    50       REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx) !u,v components of the wind (ms-1)
    51       REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)!temperature (K) and tracer concentration (kg/kg)
    52       REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx) ! altitude at the middle of the layers
    53       REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
    54       LOGICAL, INTENT(IN) :: firstcall ! =.true. if this is the first time the routine is called.
    55                                        ! Used to detect and store which tracer is co2
    56       REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx) ! wind velocity change from routines called
     59      REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) !intermediate pressure levels (Pa)
     60      REAL, INTENT(IN) :: pplay(ngrid,nlayer) !Pressure at the middle of the layers (Pa)
     61      REAL, INTENT(IN) :: pphi(ngrid,nlayer) !Geopotential at the middle of the layers (m2s-2)
     62      REAL, INTENT(IN) :: pu(ngrid,nlayer),pv(ngrid,nlayer) !u,v components of the wind (ms-1)
     63      REAL, INTENT(IN) :: pt(ngrid,nlayer),pq(ngrid,nlayer,nq)!temperature (K) and tracer concentration (kg/kg)
     64      REAL, INTENT(IN) :: zzlay(ngrid,nlayer) ! altitude at the middle of the layers
     65      REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries
     66      LOGICAL, INTENT(IN) :: tracer ! =.true. if tracers are present and to be transported
     67      INTEGER, INTENT(IN) :: igcm_co2 ! index of the CO2 tracer in mixing ratio array
     68                                      ! --> 0 if no tracer is CO2 (or no tracer at all)
     69                                      ! --> this prepares special treatment for polar night mixing
     70                                      !       (see thermcell_main_mars)
     71      REAL, INTENT(IN) :: pdu(ngrid,nlayer),pdv(ngrid,nlayer) ! wind velocity change from routines called
    5772                                                                      ! before thermals du/dt (m/s/s)
    58       REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx) ! tracer concentration change from routines called
     73      REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tracer concentration change from routines called
    5974                                                     ! before thermals dq/dt (kg/kg/s)
    60       REAL, INTENT(IN) :: pdt(ngridmx,nlayermx) ! temperature change from routines called before thermals dT/dt (K/s)
    61       REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1) ! turbulent kinetic energy
    62       REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx) ! ratio of pressure at middle of layer to surface pressure,
     75      REAL, INTENT(IN) :: pdt(ngrid,nlayer) ! temperature change from routines called before thermals dT/dt (K/s)
     76      REAL, INTENT(IN) :: q2(ngrid,nlayer+1) ! turbulent kinetic energy
     77      REAL, INTENT(IN) :: zpopsk(ngrid,nlayer) ! ratio of pressure at middle of layer to surface pressure,
    6378                                                   ! to the power r/cp, i.e. zpopsk=(pplay(ig,l)/pplev(ig,1))**rcp
    64       REAL, INTENT(IN) :: pdhdif(ngridmx,nlayermx) ! potential temperature change from turbulent diffusion scheme dT/dt (K/s)
    65       REAL, INTENT(IN) :: sensibFlux(ngridmx) ! sensible heat flux computed from surface layer scheme
     79      REAL, INTENT(IN) :: pdhdif(ngrid,nlayer) ! potential temperature change from turbulent diffusion scheme dT/dt (K/s)
     80      REAL, INTENT(IN) :: sensibFlux(ngrid) ! sensible heat flux computed from surface layer scheme
    6681
    6782!--------------------------------------------------------
    … …  
    6984!--------------------------------------------------------
    7085
    71       REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx) ! wind velocity change from thermals du/dt (m/s/s)
    72       REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx) ! wind velocity change from thermals dv/dt (m/s/s)
    73       REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx) ! temperature change from thermals dT/dt (K/s)
    74       REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx) ! tracer change from thermals dq/dt (kg/kg/s)
    75       INTEGER, INTENT(OUT) :: lmax(ngridmx) ! layer number reacher by thermals in grid point
    76       REAL, INTENT(OUT) :: zmaxth(ngridmx) ! equivalent to lmax, but in (m), interpolated
    77       REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1) ! turbulent kinetic energy change from thermals dtke/dt
    78       REAL, INTENT(OUT) :: wstar(ngridmx) ! free convection velocity (m/s)
    79 
    80 !--------------------------------------------------------
    81 ! Control parameters of the thermal model
    82 !--------------------------------------------------------
    83 ! such variables are:
    84 !  - outptherm to control internal diagnostics
    85 !  - qtransport_thermals to control tracer transport in thermals
    86 !  - dtke_thermals to control turbulent kinetic energy transport in thermals
    87      LOGICAL,SAVE :: outptherm,qtransport_thermals,dtke_thermals
     86      REAL, INTENT(OUT) :: pdu_th(ngrid,nlayer) ! wind velocity change from thermals du/dt (m/s/s)
     87      REAL, INTENT(OUT) :: pdv_th(ngrid,nlayer) ! wind velocity change from thermals dv/dt (m/s/s)
     88      REAL, INTENT(OUT) :: pdt_th(ngrid,nlayer) ! temperature change from thermals dT/dt (K/s)
     89      REAL, INTENT(OUT) :: pdq_th(ngrid,nlayer,nq) ! tracer change from thermals dq/dt (kg/kg/s)
     90      INTEGER, INTENT(OUT) :: lmax(ngrid) ! layer number reacher by thermals in grid point
     91      REAL, INTENT(OUT) :: zmaxth(ngrid) ! equivalent to lmax, but in (m), interpolated
     92      REAL, INTENT(OUT) :: pbl_dtke(ngrid,nlayer+1) ! turbulent kinetic energy change from thermals dtke/dt
     93      REAL, INTENT(OUT) :: wstar(ngrid) ! free convection velocity (m/s)
    8894
    8995!--------------------------------------------------------
    9096! Thermals local variables
    9197!--------------------------------------------------------
    92       REAL zu(ngridmx,nlayermx), zv(ngridmx,nlayermx)
    93       REAL zt(ngridmx,nlayermx)
    94       REAL d_t_ajs(ngridmx,nlayermx)
    95       REAL d_u_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
    96       REAL d_v_ajs(ngridmx,nlayermx)
    97       REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
    98       REAL detr_therm(ngridmx,nlayermx),detrmod(ngridmx,nlayermx)
    99       REAL zw2(ngridmx,nlayermx+1)
    100       REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
    101       REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
    102       REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
    103       REAL lmax_real(ngridmx)
    104       REAL masse(ngridmx,nlayermx)
     98      REAL zu(ngrid,nlayer), zv(ngrid,nlayer)
     99      REAL zt(ngrid,nlayer)
     100      REAL d_t_ajs(ngrid,nlayer)
     101      REAL d_u_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq)
     102      REAL d_v_ajs(ngrid,nlayer)
     103      REAL fm_therm(ngrid,nlayer+1), entr_therm(ngrid,nlayer)
     104      REAL detr_therm(ngrid,nlayer),detrmod(ngrid,nlayer)
     105      REAL zw2(ngrid,nlayer+1)
     106      REAL fraca(ngrid,nlayer+1),zfraca(ngrid,nlayer+1)
     107      REAL q_therm(ngrid,nlayer), pq_therm(ngrid,nlayer,nq)
     108      REAL q2_therm(ngrid,nlayer), dq2_therm(ngrid,nlayer)
     109      REAL lmax_real(ngrid)
     110      REAL masse(ngrid,nlayer)
    105111
    106112      INTEGER l,ig,iq,ii(1),k
    … …  
    111117!--------------------------------------------------------
    112118
    113       REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
     119      REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq)
    114120      INTEGER isplit
    115       INTEGER,SAVE :: nsplit_thermals
    116       REAL, SAVE :: r_aspect_thermals
    117121      REAL fact
    118       REAL zfm_therm(ngridmx,nlayermx+1),zdt
    119       REAL zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
    120       REAL zheatFlux(ngridmx,nlayermx)
    121       REAL zheatFlux_down(ngridmx,nlayermx)
    122       REAL zbuoyancyOut(ngridmx,nlayermx)
    123       REAL zbuoyancyEst(ngridmx,nlayermx)
    124       REAL zzw2(ngridmx,nlayermx+1)
    125       REAL zmax(ngridmx)
     122      REAL zfm_therm(ngrid,nlayer+1),zdt
     123      REAL zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer)
     124      REAL zheatFlux(ngrid,nlayer)
     125      REAL zheatFlux_down(ngrid,nlayer)
     126      REAL zbuoyancyOut(ngrid,nlayer)
     127      REAL zbuoyancyEst(ngrid,nlayer)
     128      REAL zzw2(ngrid,nlayer+1)
     129      REAL zmax(ngrid)
    126130      INTEGER ndt,zlmax
    127131
    … …  
    130134!--------------------------------------------------------
    131135
    132       REAL heatFlux(ngridmx,nlayermx)
    133       REAL heatFlux_down(ngridmx,nlayermx)
    134       REAL buoyancyOut(ngridmx,nlayermx)
    135       REAL buoyancyEst(ngridmx,nlayermx)
    136       REAL hfmax(ngridmx),wmax(ngridmx)
    137       REAL pbl_teta(ngridmx),dteta(ngridmx,nlayermx)
    138       REAL rpdhd(ngridmx,nlayermx)
    139       REAL wtdif(ngridmx,nlayermx),rho(ngridmx,nlayermx)
    140       REAL wtth(ngridmx,nlayermx)
    141 
    142 !--------------------------------------------------------
    143 ! Theta_m
    144 !--------------------------------------------------------
    145 
    146       INTEGER,SAVE :: ico2
    147 
    148 
    149       if(firstcall) then
    150         ico2=0
    151 ! **********************************************************************
    152 ! Polar night mixing : theta_m
    153 ! Get index of co2 tracer in tracer array
    154 ! **********************************************************************
    155         if (tracer) then !NB: tracer==.true. if tracers are present and to be transported
    156 !     Prepare Special treatment if one of the tracers is CO2 gas
    157            do iq=1,nqmx
    158              if (noms(iq).eq."co2") then
    159                 ico2=iq
    160              end if
    161            enddo
    162         endif
    163      
    164         ! control parameters of the thermal model:
    165         qtransport_thermals=.true. !! default setting: transport tracers in thermals
    166        
    167         dtke_thermals=.false. !! default setting
    168         !  switch to transport TKE in thermals. still experimental, for testing purposes only. not used in current thermal plume models both on Earth and Mars.
    169 
    170         outptherm=.false. !! default setting
    171         ! switch to control if internal quantities ofthe thermals must be output (typically set to .false., but very useful for fine diagnostics.
    172         ! get value (if set by user) of this parameter from run.def file
    173         call getin("outptherm",outptherm)
    174 
    175       endif !of if firstcall
     136      REAL heatFlux(ngrid,nlayer)
     137      REAL heatFlux_down(ngrid,nlayer)
     138      REAL buoyancyOut(ngrid,nlayer)
     139      REAL buoyancyEst(ngrid,nlayer)
     140      REAL hfmax(ngrid),wmax(ngrid)
     141      REAL pbl_teta(ngrid),dteta(ngrid,nlayer)
     142      REAL rpdhd(ngrid,nlayer)
     143      REAL wtdif(ngrid,nlayer),rho(ngrid,nlayer)
     144      REAL wtth(ngrid,nlayer)
    176145
    177146! **********************************************************************
    … …  
    226195
    227196       IF(dtke_thermals) THEN
    228           DO l=1,nlayermx
     197          DO l=1,nlayer
    229198              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
    230199          ENDDO
    231200       ENDIF
    232201
    233 
    234 ! **********************************************************************
    235 ! **********************************************************************
    236 ! **********************************************************************
    237 ! CALLTHERM
    238 ! **********************************************************************
    239 ! **********************************************************************
    240 ! **********************************************************************
    241 
    242 !         r_aspect_thermals     ! Mainly control the shape of the temperature profile
    243                                 ! in the surface layer. Decreasing it goes toward
    244                                 ! a convective-adjustment like profile.
    245 !         nsplit_thermals       ! Sub-timestep for the thermals. Very dependant on the
    246                                 ! chosen timestep for the radiative transfer.
    247                                 ! It is recommended to run with 96 timestep per day and
    248                                 ! call radiative transfer at each timestep,
    249                                 ! configuration in which thermals can run
    250                                 ! very well with a sub-timestep of 10.
    251          IF (firstcall) THEN
    252             r_aspect_thermals=1.  ! same value is OK for GCM and mesoscale
    253 ! Sub Timestep for the mesoscale model:
    254 #ifdef MESOSCALE
    255             !! valid for timesteps < 200s
    256             nsplit_thermals=4
    257 #else
    258 ! Sub Timestep for the GCM:
    259             ! If there is at least 96 timesteps per day in the gcm, subtimestep factor 10
    260             IF ((ptimestep .le. 3699.*24./96.) .and. (iradia .eq. 1)) THEN
    261                nsplit_thermals=10
    262             ELSE !if there is less, subtimestep factor 35
    263                nsplit_thermals=35
    264             ENDIF
    265 #endif
    266          ENDIF
    267 
    268 ! **********************************************************************
     202! **********************************************************************
     203! --> CALLTHERM
    269204! SUB-TIMESTEP LOOP
    270205! **********************************************************************
    … …  
    280215         lmax(:)=0
    281216
    282          if (nqmx .ne. 0 .and. ico2 .ne. 0) then !initialize co2 tracer tendancy
    283             d_q_the(:,:,ico2)=0.
     217         if (nq .ne. 0 .and. igcm_co2 .ne. 0) then !initialize co2 tracer tendancy
     218            d_q_the(:,:,igcm_co2)=0.
    284219         endif
    285220
    286221! CALL to main thermal routine
    287              CALL thermcell_main_mars(zdt  &
     222             CALL thermcell_main_mars(ngrid,nlayer,nq &
     223     &      ,tracer,igcm_co2 &
     224     &      ,zdt  &
    288225     &      ,pplay,pplev,pphi,zzlev,zzlay  &
    289226     &      ,zu,zv,zt,pq_therm,q2_therm  &
    290227     &      ,d_t_the,d_q_the  &
    291228     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
    292      &      ,r_aspect_thermals &
    293229     &      ,zzw2,fraca,zpopsk &
    294230     &      ,zheatFlux,zheatFlux_down &
    … …  
    299235! Update thermals tendancies
    300236            d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact !temperature
    301             if (ico2 .ne. 0) then
    302                d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact !co2 mass mixing ratio
     237            if (igcm_co2 .ne. 0) then
     238               d_q_the(:,:,igcm_co2)=d_q_the(:,:,igcm_co2)*ptimestep*fact !co2 mass mixing ratio
    303239            endif
    304240            zmaxth(:)=zmaxth(:)+zmax(:)*fact !thermals height
    … …  
    323259! Save tendancies
    324260           d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) !temperature tendancy (delta T)
    325             if (ico2 .ne. 0) then
    326                d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2) !tracer tendancy (delta q)
     261            if (igcm_co2 .ne. 0) then
     262               d_q_ajs(:,:,igcm_co2)=d_q_ajs(:,:,igcm_co2)+d_q_the(:,:,igcm_co2) !tracer tendancy (delta q)
    327263            endif
    328264
    329265!  Increment temperature and co2 concentration for next pass in subtimestep loop
    330266            zt(:,:) = zt(:,:) + d_t_the(:,:) !temperature
    331             if (ico2 .ne. 0) then
    332              pq_therm(:,:,ico2) = &
    333      &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2) !co2 tracer
     267            if (igcm_co2 .ne. 0) then
     268             pq_therm(:,:,igcm_co2) = &
     269     &          pq_therm(:,:,igcm_co2) + d_q_the(:,:,igcm_co2) !co2 tracer
    334270            endif
    335271
    … …  
    342278      lmax(:)=nint(lmax_real(:))
    343279      zlmax=MAXVAL(lmax(:))+2
    344       if (zlmax .ge. nlayermx) then
     280      if (zlmax .ge. nlayer) then
    345281        print*,'thermals have reached last layer of the model'
    346282        print*,'this is not good !'
    … …  
    356292
    357293! mass of cells
    358       do l=1,nlayermx
     294      do l=1,nlayer
    359295         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
    360296      enddo
    … …  
    364300      detrmod(:,:)=0.
    365301      do l=1,zlmax
    366          do ig=1,ngridmx
     302         do ig=1,ngrid
    367303            detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
    368304     &      +entr_therm(ig,l)
    … …  
    377313      ! result is a derivative (d_u_ajs in m/s/s)
    378314      ndt=10
    379       call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
     315      call thermcell_dqup(ngrid,nlayer,ptimestep                &
    380316     &      ,fm_therm,entr_therm,detrmod,  &
    381317     &     masse,zu,d_u_ajs,ndt,zlmax)
    … …  
    383319      ! v component of wind velocity advection in thermals
    384320      ! result is a derivative (d_v_ajs in m/s/s)
    385       call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
     321      call thermcell_dqup(ngrid,nlayer,ptimestep    &
    386322     &       ,fm_therm,entr_therm,detrmod,  &
    387323     &     masse,zv,d_v_ajs,ndt,zlmax)
    … …  
    390326      ! result is a derivative (d_q_ajs in kg/kg/s)
    391327
    392       if (nqmx .ne. 0.) then
    393       DO iq=1,nqmx
    394       if (iq .ne. ico2) then
    395       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
     328      if (nq .ne. 0.) then
     329      DO iq=1,nq
     330      if (iq .ne. igcm_co2) then
     331      call thermcell_dqup(ngrid,nlayer,ptimestep     &
    396332     &     ,fm_therm,entr_therm,detrmod,  &
    397333     &    masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),ndt,zlmax)
    … …  
    403339      ! result is a tendancy (d_u_ajs in J)
    404340      if (dtke_thermals) then
    405       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
     341      call thermcell_dqup(ngrid,nlayer,ptimestep     &
    406342     &     ,fm_therm,entr_therm,detrmod,  &
    407343     &    masse,q2_therm,dq2_therm,ndt,zlmax)
    … …  
    409345
    410346      ! compute wmax for diagnostics
    411       DO ig=1,ngridmx
     347      DO ig=1,ngrid
    412348         wmax(ig)=MAXVAL(zw2(ig,:))
    413349      ENDDO
    … …  
    434370           if(qtransport_thermals) then
    435371              if(tracer) then
    436                do iq=1,nqmx
    437                 if (iq .ne. ico2) then
     372               do iq=1,nq
     373                if (iq .ne. igcm_co2) then
    438374                  do l=1,zlmax
    439375                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq) !non-co2 tracers d_q_ajs are dq/dt (kg/kg/s)
    … …  
    450386           ! if tke is transported in thermals, update output variable, else this is 0.
    451387           IF(dtke_thermals) THEN
    452               DO l=2,nlayermx
     388              DO l=2,nlayer
    453389                 pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))
    454390              ENDDO
    455391 
    456392              pbl_dtke(:,1)=0.5*dq2_therm(:,1)
    457               pbl_dtke(:,nlayermx+1)=0.
     393              pbl_dtke(:,nlayer+1)=0.
    458394           ENDIF
    459395
    … …  
    475411! Potential temperature gradient
    476412
    477       dteta(:,nlayermx)=0.
    478       DO l=1,nlayermx-1
    479          DO ig=1, ngridmx
     413      dteta(:,nlayer)=0.
     414      DO l=1,nlayer-1
     415         DO ig=1, ngrid
    480416            dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l))          &
    481417     &              /(zzlay(ig,l+1)-zzlay(ig,l))
    … …  
    485421! Computation of the PBL mixed layer temperature
    486422
    487       DO ig=1, ngridmx
     423      DO ig=1, ngrid
    488424         ii=MINLOC(abs(dteta(ig,1:lmax(ig))))
    489425         pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1))
    … …  
    504440      rpdhd(:,:)=0.
    505441
    506       DO ig=1,ngridmx
     442      DO ig=1,ngrid
    507443       DO l=1,lmax(ig)
    508444        rpdhd(ig,l)=0.
    … …  
    527463! integrate -rho*pdhdif
    528464
    529       DO ig=1,ngridmx
     465      DO ig=1,ngrid
    530466       DO l=1,lmax(ig)
    531467        rpdhd(ig,l)=0.
    … …  
    538474       ENDDO
    539475      ENDDO
    540       rpdhd(:,nlayermx)=0.
     476      rpdhd(:,nlayer)=0.
    541477
    542478! compute w'teta' from thermals
    … …  
    547483! and compute the maximum
    548484
    549       DO ig=1,ngridmx
     485      DO ig=1,ngrid
    550486        hfmax(ig)=MAXVAL(wtth(ig,:)+wtdif(ig,:))
    551487      ENDDO
    … …  
    555491! ------------
    556492
    557       DO ig=1, ngridmx
     493      DO ig=1, ngrid
    558494         IF (zmax(ig) .gt. 0.) THEN
    559495            wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.)
    … …  
    563499      ENDDO
    564500
    565 ! **********************************************************************
    566 ! Diagnostics
    567 ! **********************************************************************
    568 ! THESE DIAGNOSTICS ARE WRITTEN IN THE LMD-GCM FORMAT
    569 ! THESE ARE LEFT AS A FURTHER INDICATION OF THE QUANTITIES DEFINED IN THIS ROUTINE
    570        
    571 ! outptherm is a logical in callkeys that controls if internal quantities of the thermals must be output.
    572 
    573         if (outptherm) then
    574         if (ngridmx .eq. 1) then !these are 1D outputs
    575         call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
    576      &                       'kg/m-2',1,entr_therm)
    577         call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
    578      &                       'kg/m-2',1,detr_therm)
    579         call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
    580      &                       'kg/m-2',1,fm_therm)
    581         call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
    582      &                       'm/s',1,zw2)
    583         call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
    584      &                       'SI',1,heatFlux)
    585        call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
    586      &                       'SI',1,heatFlux_down)
    587         call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
    588      &                       'percent',1,fraca)
    589         call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
    590      &                       'm.s-2',1,buoyancyOut)
    591         call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
    592      &                       'm.s-2',1,buoyancyEst)
    593         call WRITEDIAGFI(ngridmx,'d_t_th',  &
    594      &         'tendance temp TH','K',1,d_t_ajs)
    595         call WRITEDIAGFI(ngridmx,'d_q_th',  &
    596      &         'tendance traceur TH','kg/kg',1,d_q_ajs)
    597         call WRITEDIAGFI(ngridmx,'zmax',  &
    598      &         'pbl height','m',0,zmaxth)
    599         call WRITEDIAGFI(ngridmx,'d_u_th',  &
    600      &         'tendance moment','m/s',1,pdu_th)
    601         call WRITEDIAGFI(ngridmx,'wtdif',  &
    602      &         'heat flux from diffusion','K.m/s',1,wtdif)
    603         call WRITEDIAGFI(ngridmx,'wtth',  &
    604      &         'heat flux from thermals','K.m/s',1,wtth)
    605         call WRITEDIAGFI(ngridmx,'wttot',  &
    606      &         'heat flux PBL','K.m/s',1,wtdif(:,:)+wtth(:,:))
    607 
    608       else !these are 3D outputs
    609 
    610         call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
    611      &                       'kg/m-2',3,entr_therm)
    612         call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
    613      &                       'kg/m-2',3,detr_therm)
    614         call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
    615      &                       'kg/m-2',3,fm_therm)
    616         call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
    617      &                       'm/s',3,zw2)
    618         call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
    619      &                       'SI',3,heatFlux)
    620         call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
    621      &                       'SI',3,buoyancyOut)
    622         call WRITEDIAGFI(ngridmx,'d_t_th',  &
    623      &         'tendance temp TH','K',3,d_t_ajs)
    624 
    625       endif
    626       endif ! of if (outptherm)
    627 
    628501       END

  • trunk/LMDZ.MARS/libf/phymars/physiq.F

    r1013 r1032  
    896896      if(calltherm) then
    897897 
    898         call calltherm_interface(firstcall,
     898        call calltherm_interface(ngrid,nlayer,nq,
     899     $ tracer,igcm_co2,
    899900     $ zzlev,zzlay,
    900901     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
    … …  
    21062107     &                        'part/kg',3,ndust)
    21072108#ifdef MESOINI
    2108              call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
    2109      &                        'kg/kg',3,qdust)
    2110              call WRITEDIAGFI(ngridmx,'dustN','Dust number',
    2111      &                        'part/kg',3,ndust)
    2112              call WRITEDIAGFI(ngridmx,'ccn','Nuclei mass mr',
    2113      &                        'kg/kg',3,qccn)
    2114              call WRITEDIAGFI(ngridmx,'ccnN','Nuclei number',
    2115      &                        'part/kg',3,nccn)
     2109!     !!! to initialize mesoscale we need scaled variables
     2110!     !!! because this must correspond to starting point for tracers
     2111!             call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
     2112!     &           'kg/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_mass))
     2113!             call WRITEDIAGFI(ngridmx,'dustN','Dust number',
     2114!     &           'part/kg',3,pq(1:ngridmx,1:nlayermx,igcm_dust_number))
     2115!             call WRITEDIAGFI(ngridmx,'ccn','Nuclei mass mr',
     2116!     &           'kg/kg',3,pq(1:ngridmx,1:nlayermx,igcm_ccn_mass))
     2117!             call WRITEDIAGFI(ngridmx,'ccnN','Nuclei number',
     2118!     &           'part/kg',3,pq(1:ngridmx,1:nlayermx,igcm_ccn_number))
     2119              call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr',
     2120     &                        'kg/kg',3,pq(1,1,igcm_dust_mass))
     2121              call WRITEDIAGFI(ngridmx,'dustN','Dust number',
     2122     &                        'part/kg',3,pq(1,1,igcm_dust_number))
     2123              call WRITEDIAGFI(ngridmx,'ccn','Nuclei mass mr',
     2124     &                        'kg/kg',3,pq(1,1,igcm_ccn_mass))
     2125              call WRITEDIAGFI(ngridmx,'ccnN','Nuclei number',
     2126     &                        'part/kg',3,pq(1,1,igcm_ccn_number))
    21162127#endif
    21172128           else

  • trunk/LMDZ.MARS/libf/phymars/thermcell_dqup.F90

    r1028 r1032  
    1 ! ---------------------------------------------------------------------
    2 ! Arnaud Colaitis 2011-01-05
    3 ! --------------------------------------------------------------------
    41      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm,entr,detr,  &
    52     &    masse0,q_therm,dq_therm,ndt,zlmax)
    … …  
    1613!
    1714!=======================================================================
    18 
    19 #include "dimensions.h" !contains global GCM grid dimension informations
    20 #include "dimphys.h" !similar to dimensions.h, and based on it; includes
    21                      ! "ngridmx": number of horizontal grid points
    22                      ! "nlayermx": number of atmospheric layers
     15! Arnaud Colaitis 2011-01-05 (modified 2011-2013)
     16! SEE COMMENTS IN CALLTHERM_INTERFACE
     17!=======================================================================
    2318
    2419! ============================ INPUTS ============================
    … …  
    2621      INTEGER, INTENT(IN) :: ngrid,nlayer ! number of grid points and number of levels
    2722      REAL, INTENT(IN) :: ptimestep ! timestep (s)
    28       REAL, INTENT(IN) :: fm(ngridmx,nlayermx+1) ! upward mass flux
    29       REAL, INTENT(IN) :: entr(ngridmx,nlayermx) ! entrainment mass flux
    30       REAL, INTENT(IN) :: detr(ngridmx,nlayermx) ! detrainment mass flux
    31       REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx) ! initial profil of q
    32       REAL, INTENT(IN) :: masse0(ngridmx,nlayermx) ! mass of cells
     23      REAL, INTENT(IN) :: fm(ngrid,nlayer+1) ! upward mass flux
     24      REAL, INTENT(IN) :: entr(ngrid,nlayer) ! entrainment mass flux
     25      REAL, INTENT(IN) :: detr(ngrid,nlayer) ! detrainment mass flux
     26      REAL, INTENT(IN) :: q_therm(ngrid,nlayer) ! initial profil of q
     27      REAL, INTENT(IN) :: masse0(ngrid,nlayer) ! mass of cells
    3328      INTEGER, INTENT(IN) :: ndt ! number of subtimesteps
    3429      INTEGER, INTENT(IN) :: zlmax ! index of maximum layer
    … …  
    3631! ============================ OUTPUTS ===========================
    3732
    38       REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)  ! dq/dt -> derivative
     33      REAL, INTENT(OUT) :: dq_therm(ngrid,nlayer)  ! dq/dt -> derivative
    3934
    4035! ============================ LOCAL =============================
    4136
    42       REAL q(ngridmx,nlayermx)
    43       REAL qa(ngridmx,nlayermx)
     37      REAL q(ngrid,nlayer)
     38      REAL qa(ngrid,nlayer)
    4439      INTEGER ig,k,i
    45       REAL invflux0(ngridmx,nlayermx)
     40      REAL invflux0(ngrid,nlayer)
    4641      REAL ztimestep
    4742
    … …  
    6560
    6661        do k=2,zlmax
    67            do ig=1,ngridmx
     62           do ig=1,ngrid
    6863              if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
    6964     &        1.e-5*masse0(ig,k)) then

  • trunk/LMDZ.MARS/libf/phymars/thermcell_main_mars.F90

    r1028 r1032  
    11!=======================================================================
    22! THERMCELL_MAIN_MARS
    3 ! Author : A Colaitis
     3! Author : A. Colaitis after C. Rio and F. Hourdin
    44!
    55! This routine is called by calltherm_interface and is inside a sub-timestep
    … …  
    77! detrainment rate as well as the source profile.
    88! Mass flux are then computed and temperature and CO2 MMR are transported.
    9 !
    10 ! This routine is based on the terrestrial version thermcell_main of the
    11 ! LMDZ model, written by C. Rio and F. Hourdin
    129!=======================================================================
    13 !
    14 !
    15       SUBROUTINE thermcell_main_mars(ptimestep  &
     10! SEE COMMENTS IN CALLTHERM_INTERFACE
     11!=======================================================================
     12!
     13!
     14      SUBROUTINE thermcell_main_mars(ngrid,nlayer,nq &
     15     &                  ,tracer,igcm_co2 &
     16     &                  ,ptimestep  &
    1617     &                  ,pplay,pplev,pphi,zlev,zlay  &
    1718     &                  ,pu,pv,pt,pq,pq2  &
    1819     &                  ,pdtadj,pdqadj  &
    1920     &                  ,fm,entr,detr,lmax,zmax  &
    20      &                  ,r_aspect &
    2121     &                  ,zw2,fraca &
    2222     &                  ,zpopsk,heatFlux,heatFlux_down &
    2323     &                  ,buoyancyOut, buoyancyEst)
    2424
     25      USE comtherm_h
     26
    2527      IMPLICIT NONE
    2628
    2729!=======================================================================
    2830
    29 #include "callkeys.h" !contains logical values determining which subroutines and which options are used.
    30                       ! includes logical "tracer", which is .true. if tracers are present and to be transported
    31                       ! includes logical "lwrite", which is .true. if one wants more verbose outputs as GCM runs
    32 #include "dimensions.h" !contains global GCM grid dimension informations
    33 #include "dimphys.h" !similar to dimensions.h, and based on it; includes
    34                      ! "ngridmx": number of horizontal grid points
    35                      ! "nlayermx": number of atmospheric layers
     31! SHARED VARIABLES. This needs adaptations in another climate model.
    3632#include "comcstfi.h" !contains physical constant values such as
    3733                      ! "g" : gravitational acceleration (m.s-2)
    3834                      ! "r" : recuced gas constant (J.K-1.mol-1)
    39 #include "tracer.h" !contains tracer information such as
    40                     ! "nqmx": number of tracers
    41                     ! "noms()": tracer name
     35
    4236! ============== INPUTS ==============
    4337
     38      INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
     39      INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
     40      INTEGER, INTENT(IN) :: nq ! number of tracer species
     41      LOGICAL, INTENT(IN) :: tracer ! =.true. if tracers are present and to be transported
     42      INTEGER, INTENT(IN) :: igcm_co2 ! index of the CO2 tracer in mixing ratio array
     43                                      ! --> 0 if no tracer is CO2 (or no tracer at all)
     44                                      ! --> this prepares special treatment for polar night mixing
    4445      REAL, INTENT(IN) :: ptimestep !subtimestep (s)
    45       REAL, INTENT(IN) :: r_aspect !aspect ratio (see paragraph 45 of paper and appendix S4)
    46       REAL, INTENT(IN) :: pt(ngridmx,nlayermx) !temperature (K)
    47       REAL, INTENT(IN) :: pu(ngridmx,nlayermx) !u component of the wind (ms-1)
    48       REAL, INTENT(IN) :: pv(ngridmx,nlayermx) !v component of the wind (ms-1)
    49       REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx) !tracer concentration (kg/kg)
    50       REAL, INTENT(IN) :: pq2(ngridmx,nlayermx) ! Turbulent Kinetic Energy
    51       REAL, INTENT(IN) :: pplay(ngridmx,nlayermx) !Pressure at the middle of the layers (Pa)
    52       REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1) !intermediate pressure levels (Pa)
    53       REAL, INTENT(IN) :: pphi(ngridmx,nlayermx) !Geopotential at the middle of the layers (m2s-2)
    54       REAL, INTENT(IN) :: zlay(ngridmx,nlayermx) ! altitude at the middle of the layers
    55       REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
     46      REAL, INTENT(IN) :: pt(ngrid,nlayer) !temperature (K)
     47      REAL, INTENT(IN) :: pu(ngrid,nlayer) !u component of the wind (ms-1)
     48      REAL, INTENT(IN) :: pv(ngrid,nlayer) !v component of the wind (ms-1)
     49      REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) !tracer concentration (kg/kg)
     50      REAL, INTENT(IN) :: pq2(ngrid,nlayer) ! Turbulent Kinetic Energy
     51      REAL, INTENT(IN) :: pplay(ngrid,nlayer) !Pressure at the middle of the layers (Pa)
     52      REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) !intermediate pressure levels (Pa)
     53      REAL, INTENT(IN) :: pphi(ngrid,nlayer) !Geopotential at the middle of the layers (m2s-2)
     54      REAL, INTENT(IN) :: zlay(ngrid,nlayer) ! altitude at the middle of the layers
     55      REAL, INTENT(IN) :: zlev(ngrid,nlayer+1) ! altitude at layer boundaries
    5656
    5757! ============== OUTPUTS ==============
    5858
    5959! TEMPERATURE
    60       REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx) !temperature change from thermals dT/dt (K/s)
     60      REAL, INTENT(OUT) :: pdtadj(ngrid,nlayer) !temperature change from thermals dT/dt (K/s)
    6161
    6262! DIAGNOSTICS
    63       REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1) ! vertical velocity (m/s)
    64       REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
    65       REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
     63      REAL, INTENT(OUT) :: zw2(ngrid,nlayer+1) ! vertical velocity (m/s)
     64      REAL, INTENT(OUT) :: heatFlux(ngrid,nlayer)   ! interface heatflux
     65      REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlayer) ! interface heat flux from downdraft
    6666
    6767! ============== LOCAL ================
    68       REAL :: pdqadj(ngridmx,nlayermx,nqmx) !tracer change from thermals dq/dt, only for CO2 (the rest can be advected outside of the loop)
     68      REAL :: pdqadj(ngrid,nlayer,nq) !tracer change from thermals dq/dt, only for CO2 (the rest can be advected outside of the loop)
    6969
    7070! dummy variables when output not needed :
    7171
    72       REAL :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
    73       REAL :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
     72      REAL :: buoyancyOut(ngrid,nlayer)  ! interlayer buoyancy term
     73      REAL :: buoyancyEst(ngrid,nlayer)  ! interlayer estimated buoyancy term
    7474
    7575! ============== LOCAL ================
    7676
    7777      INTEGER ig,k,l,ll,iq
    78       INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
    79       REAL zmax(ngridmx)
    80       REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx)
    81       REAL zh(ngridmx,nlayermx)
    82       REAL zdthladj(ngridmx,nlayermx)
    83       REAL zdthladj_down(ngridmx,nlayermx)
    84       REAL ztvd(ngridmx,nlayermx)
    85       REAL ztv(ngridmx,nlayermx)
    86       REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx),zo(ngridmx,nlayermx)
    87       REAL zva(ngridmx,nlayermx)
    88       REAL zua(ngridmx,nlayermx)
    89 
    90       REAL zta(ngridmx,nlayermx)
    91       REAL fraca(ngridmx,nlayermx+1)
    92       REAL q2(ngridmx,nlayermx)
    93       REAL rho(ngridmx,nlayermx),rhobarz(ngridmx,nlayermx),masse(ngridmx,nlayermx)
    94       REAL zpopsk(ngridmx,nlayermx)
    95 
    96       REAL wmax(ngridmx)
    97       REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
    98 
    99       REAL fm_down(ngridmx,nlayermx+1)
    100 
    101       REAL ztla(ngridmx,nlayermx)
    102 
    103       REAL f_star(ngridmx,nlayermx+1),entr_star(ngridmx,nlayermx)
    104       REAL detr_star(ngridmx,nlayermx)
    105       REAL alim_star_tot(ngridmx)
    106       REAL alim_star(ngridmx,nlayermx)
    107       REAL alim_star_clos(ngridmx,nlayermx)
    108       REAL f(ngridmx)
    109 
    110       REAL detrmod(ngridmx,nlayermx)
    111 
    112       REAL teta_th_int(ngridmx,nlayermx)
    113       REAL teta_env_int(ngridmx,nlayermx)
    114       REAL teta_down_int(ngridmx,nlayermx)
     78      INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid)
     79      REAL zmax(ngrid)
     80      REAL ztva(ngrid,nlayer),zw_est(ngrid,nlayer+1),ztva_est(ngrid,nlayer)
     81      REAL zh(ngrid,nlayer)
     82      REAL zdthladj(ngrid,nlayer)
     83      REAL zdthladj_down(ngrid,nlayer)
     84      REAL ztvd(ngrid,nlayer)
     85      REAL ztv(ngrid,nlayer)
     86      REAL zu(ngrid,nlayer),zv(ngrid,nlayer),zo(ngrid,nlayer)
     87      REAL zva(ngrid,nlayer)
     88      REAL zua(ngrid,nlayer)
     89
     90      REAL zta(ngrid,nlayer)
     91      REAL fraca(ngrid,nlayer+1)
     92      REAL q2(ngrid,nlayer)
     93      REAL rho(ngrid,nlayer),rhobarz(ngrid,nlayer),masse(ngrid,nlayer)
     94      REAL zpopsk(ngrid,nlayer)
     95
     96      REAL wmax(ngrid)
     97      REAL fm(ngrid,nlayer+1),entr(ngrid,nlayer),detr(ngrid,nlayer)
     98
     99      REAL fm_down(ngrid,nlayer+1)
     100
     101      REAL ztla(ngrid,nlayer)
     102
     103      REAL f_star(ngrid,nlayer+1),entr_star(ngrid,nlayer)
     104      REAL detr_star(ngrid,nlayer)
     105      REAL alim_star_tot(ngrid)
     106      REAL alim_star(ngrid,nlayer)
     107      REAL alim_star_clos(ngrid,nlayer)
     108      REAL f(ngrid)
     109
     110      REAL detrmod(ngrid,nlayer)
     111
     112      REAL teta_th_int(ngrid,nlayer)
     113      REAL teta_env_int(ngrid,nlayer)
     114      REAL teta_down_int(ngrid,nlayer)
    115115
    116116      CHARACTER (LEN=80) :: abort_message
    … …  
    119119! ============= PLUME VARIABLES ============
    120120
    121       REAL w_est(ngridmx,nlayermx+1)
    122       REAL wa_moy(ngridmx,nlayermx+1)
    123       REAL wmaxa(ngridmx)
    124       REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
    125       LOGICAL activecell(ngridmx),activetmp(ngridmx)
     121      REAL w_est(ngrid,nlayer+1)
     122      REAL wa_moy(ngrid,nlayer+1)
     123      REAL wmaxa(ngrid)
     124      REAL zdz,zbuoy(ngrid,nlayer),zw2m
     125      LOGICAL activecell(ngrid),activetmp(ngrid)
    126126      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega
    127127      INTEGER tic
    … …  
    131131! ============= HEIGHT VARIABLES ===========
    132132
    133       REAL num(ngridmx)
    134       REAL denom(ngridmx)
    135       REAL zlevinter(ngridmx)
     133      REAL num(ngrid)
     134      REAL denom(ngrid)
     135      REAL zlevinter(ngrid)
    136136      INTEGER zlmax
    137137
    … …  
    140140! ============= CLOSURE VARIABLES =========
    141141
    142       REAL zdenom(ngridmx)
    143       REAL alim_star2(ngridmx)
    144       REAL alim_star_tot_clos(ngridmx)
     142      REAL zdenom(ngrid)
     143      REAL alim_star2(ngrid)
     144      REAL alim_star_tot_clos(ngrid)
    145145      INTEGER llmax
    146146
    … …  
    159159! ============== Theta_M Variables ========
    160160
    161       INTEGER ico2
    162       SAVE ico2
    163161      REAL m_co2, m_noco2, A , B
    164162      SAVE A, B
    165       LOGICAL firstcall
    166       save firstcall
    167       data firstcall/.true./
    168       REAL zhc(ngridmx,nlayermx)
    169       REAL ratiom(ngridmx,nlayermx)
     163      REAL zhc(ngrid,nlayer)
     164      REAL ratiom(ngrid,nlayer)
    170165
    171166! =========================================
    … …  
    181176      ndt=1
    182177
     178!.......................................................................
     179!  Special treatment for co2:
     180!.......................................................................
    183181! **********************************************************************
    184182! In order to take into account the effect of vertical molar mass
    … …  
    193191! This is especially important for modelling polar convection.
    194192! **********************************************************************
    195 
    196      
    197 !.......................................................................
    198 !  Special treatment for co2 at firstcall: compute coefficients and
    199 !  get index of co2 tracer
    200 !.......................................................................
    201 
    202       if(firstcall) then
    203         ico2=0
    204         if (tracer) then
    205 !     Prepare Special treatment if one of the tracers is CO2 gas
    206            do iq=1,nqmx
    207              if (noms(iq).eq."co2") then
    208                 ico2=iq
    209                 m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
    210                 m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
    211 !               Compute A and B coefficient use to compute
    212 !               mean molecular mass Mair defined by
    213 !               1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2
    214 !               1/Mair = A*q(ico2) + B
    215                 A =(1/m_co2 - 1/m_noco2)
    216                 B=1/m_noco2
    217              end if
    218            enddo
    219         endif
    220 
    221       firstcall=.false.
    222       endif !of if firstcall
    223 
    224 !.......................................................................
    225 !  Special treatment for co2
    226 !.......................................................................
    227 
    228       if (ico2.ne.0) then
     193      if (igcm_co2.ne.0) then
     194
     195         m_co2 = 44.01E-3  ! CO2 molecular mass (kg/mol)
     196         m_noco2 = 33.37E-3  ! Non condensible mol mass (kg/mol)
     197         ! Compute A and B coefficient use to compute
     198         ! mean molecular mass Mair defined by
     199         ! 1/Mair = q(igcm_co2)/m_co2 + (1-q(igcm_co2))/m_noco2
     200         ! 1/Mair = A*q(igcm_co2) + B
     201         A =(1/m_co2 - 1/m_noco2)
     202         B=1/m_noco2
     203
    229204!     Special case if one of the tracers is CO2 gas
    230          DO l=1,nlayermx
    231            DO ig=1,ngridmx
    232             ztv(ig,l) = zhc(ig,l)*(A*pq(ig,l,ico2)+B)
     205         DO l=1,nlayer
     206           DO ig=1,ngrid
     207            ztv(ig,l) = zhc(ig,l)*(A*pq(ig,l,igcm_co2)+B)
    233208           ENDDO
    234209         ENDDO
    … …  
    265240      rhobarz(:,1)=rho(:,1)
    266241
    267       do l=2,nlayermx
     242      do l=2,nlayer
    268243          rhobarz(:,l)=pplev(:,l)/(r*0.5*(pt(:,l)+pt(:,l-1)))
    269244      enddo
    270245
    271246! mass computation
    272       do l=1,nlayermx
     247      do l=1,nlayer
    273248         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
    274249      enddo
    … …  
    395370
    396371      activecell(:)=ztv(:,1)>ztv(:,2)
    397           do ig=1,ngridmx
     372          do ig=1,ngrid
    398373            if (ztv(ig,1)>=(ztv(ig,2))) then
    399374               alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
    … …  
    404379         enddo
    405380
    406        do l=2,nlayermx-1
    407          do ig=1,ngridmx
     381       do l=2,nlayer-1
     382         do ig=1,ngrid
    408383           if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) &
    409384     &             .and. (alim_star(ig,l-1).ne. 0.)) then
    … …  
    415390        enddo
    416391      enddo
    417       do l=1,nlayermx
    418          do ig=1,ngridmx
     392      do l=1,nlayer
     393         do ig=1,ngrid
    419394            if (alim_star_tot(ig) > 1.e-10 ) then
    420395               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    … …  
    428403! (f_star) in the first and second layer from the source profile.
    429404
    430       do ig=1,ngridmx
     405      do ig=1,ngrid
    431406          if (activecell(ig)) then
    432407          ztla(ig,1)=ztv(ig,1)
    … …  
    445420! LOOP ON VERTICAL LEVELS
    446421!==============================================================================
    447       do l=2,nlayermx-1
     422      do l=2,nlayer-1
    448423!==============================================================================
    449424!==============================================================================
    … …  
    452427
    453428! is the thermal plume still active ?
    454           do ig=1,ngridmx
     429          do ig=1,ngrid
    455430             activecell(ig)=activecell(ig) &
    456431      &                 .and. zw2(ig,l)>1.e-10 &
    … …  
    470445!---------------------------------------------------------------------------
    471446
    472           do ig=1,ngridmx
     447          do ig=1,ngrid
    473448             if(activecell(ig)) then
    474449                ztva_est(ig,l)=ztla(ig,l-1)
    … …  
    496471!-------------------------------------------------
    497472
    498       do ig=1,ngridmx
     473      do ig=1,ngrid
    499474        if (activecell(ig)) then
    500475
    … …  
    573548
    574549! If the cell is active, compute temperature inside updraft
    575       do ig=1,ngridmx
     550      do ig=1,ngrid
    576551       if (activetmp(ig)) then
    577552
    … …  
    587562
    588563! Compute new buoyancu and vertical velocity
    589       do ig=1,ngridmx
     564      do ig=1,ngrid
    590565      if (activetmp(ig)) then
    591566           ztva(ig,l) = ztla(ig,l)
    … …  
    606581         ! ND detrainment rate, see paragraph 44 of paper
    607582
    608       do ig=1,ngridmx
     583      do ig=1,ngrid
    609584        if (activetmp(ig)) then
    610585
    … …  
    664639!---------------------------------------------------------------------------
    665640
    666       do ig=1,ngridmx
     641      do ig=1,ngrid
    667642            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    668       IF (lwrite) THEN
     643      IF (thermverbose) THEN
    669644           print*,'thermcell_plume, particular case in velocity profile'
    670645      ENDIF
    … …  
    695670! alim_star.
    696671
    697        do ig=1,ngridmx
     672       do ig=1,ngrid
    698673          alim_star_tot(ig)=0.
    699674       enddo
    700        do ig=1,ngridmx
     675       do ig=1,ngrid
    701676          do l=1,lalim(ig)-1
    702677          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    … …  
    704679       enddo
    705680
    706       do l=1,nlayermx
    707          do ig=1,ngridmx
     681      do l=1,nlayer
     682         do ig=1,ngrid
    708683            if (alim_star_tot(ig) > 1.e-10 ) then
    709684               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    … …  
    727702
    728703! Index of the thermal plume height
    729       do ig=1,ngridmx
     704      do ig=1,ngrid
    730705         lmax(ig)=lalim(ig)
    731706      enddo
    732       do ig=1,ngridmx
    733          do l=nlayermx,lalim(ig)+1,-1
     707      do ig=1,ngrid
     708         do l=nlayer,lalim(ig)+1,-1
    734709            if (zw2(ig,l).le.1.e-10) then
    735710               lmax(ig)=l-1
    … …  
    739714
    740715! Particular case when the thermal reached the model top, which is not a good sign
    741       do ig=1,ngridmx
    742       if ( zw2(ig,nlayermx) > 1.e-10 ) then
     716      do ig=1,ngrid
     717      if ( zw2(ig,nlayer) > 1.e-10 ) then
    743718          print*,'WARNING !!!!! W2 non-zero in last layer'
    744           lmax(ig)=nlayermx
     719          lmax(ig)=nlayer
    745720      endif
    746721      enddo
    747722
    748723! Maximum vertical velocity zw2
    749       do ig=1,ngridmx
     724      do ig=1,ngrid
    750725         wmax(ig)=0.
    751726      enddo
    752727
    753       do l=1,nlayermx
    754          do ig=1,ngridmx
     728      do l=1,nlayer
     729         do ig=1,ngrid
    755730            if (l.le.lmax(ig)) then
    756731                if (zw2(ig,l).lt.0.)then
    … …  
    769744! zmax=Integral[z*w(z)*dz]/Integral[w(z)*dz]
    770745!
    771       do  ig=1,ngridmx
     746      do  ig=1,ngrid
    772747         zmax(ig)=0.
    773748         zlevinter(ig)=zlev(ig,1)
    … …  
    776751         num(:)=0.
    777752         denom(:)=0.
    778          do ig=1,ngridmx
    779           do l=1,nlayermx
     753         do ig=1,ngrid
     754          do l=1,nlayer
    780755             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    781756             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    782757          enddo
    783758       enddo
    784        do ig=1,ngridmx
     759       do ig=1,ngrid
    785760       if (denom(ig).gt.1.e-10) then
    786761          zmax(ig)=2.*num(ig)/denom(ig)
    … …  
    793768
    794769      zlmax=MAXVAL(lmax(:))+2
    795       if (zlmax .ge. nlayermx) then
     770      if (zlmax .ge. nlayer) then
    796771        print*,'thermals have reached last layer of the model'
    797772        print*,'this is not good !'
    … …  
    819794! llmax is the index of the heighest thermal in the simulation domain
    820795      llmax=1
    821       do ig=1,ngridmx
     796      do ig=1,ngrid
    822797         if (lalim(ig)>llmax) llmax=lalim(ig)
    823798      enddo
    … …  
    826801
    827802      do k=1,llmax-1
    828          do ig=1,ngridmx
     803         do ig=1,ngrid
    829804            if (k<lalim(ig)) then
    830805         alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)*alim_star_clos(ig,k)  &
    … …  
    837812! Closure mass flux, equation 13 of appendix 4.2 in paper
    838813
    839       do ig=1,ngridmx
     814      do ig=1,ngrid
    840815         if (alim_star2(ig)>1.e-10) then
    841816             f(ig)=wmax(ig)*alim_star_tot_clos(ig)/  &
    842       &     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
     817      &     (max(500.,zmax(ig))*r_aspect_thermals*alim_star2(ig))
    843818
    844819         endif
    … …  
    864839
    865840
    866 !    these variables allow to follow corrections made to the mass flux when lwrite=.true.
     841!    these variables allow to follow corrections made to the mass flux when thermverbose=.true.
    867842      ncorecfm1=0
    868843      ncorecfm2=0
    … …  
    887862
    888863      do l=1,zlmax
    889          do ig=1,ngridmx
     864         do ig=1,ngrid
    890865            if (l.lt.lmax(ig)) then
    891866               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
    … …  
    918893! Upward mass flux at level l+1
    919894
    920          do ig=1,ngridmx
     895         do ig=1,ngrid
    921896            if (l.lt.lmax(ig)) then
    922897               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
    … …  
    934909!-------------------------------------------------------------------------
    935910
    936          do ig=1,ngridmx
     911         do ig=1,ngrid
    937912
    938913            if (fm(ig,l+1).lt.0.) then
    … …  
    943918               ncorecfm2=ncorecfm2+1
    944919               else
    945                IF (lwrite) THEN
     920               IF (thermverbose) THEN
    946921          print*,'fm(l+1)<0 : ig, l+1,lmax :',ig,l+1,lmax(ig),fm(ig,l+1)
    947922               ENDIF
    … …  
    958933!-------------------------------------------------------------------------
    959934
    960          do ig=1,ngridmx
     935         do ig=1,ngrid
    961936            if (detr(ig,l).gt.fm(ig,l)) then
    962937               ncorecfm6=ncorecfm6+1
    … …  
    981956!-------------------------------------------------------------------------
    982957
    983          do ig=1,ngridmx
     958         do ig=1,ngrid
    984959            if (fm(ig,l+1).lt.0.) then
    985960               detr(ig,l)=detr(ig,l)+fm(ig,l+1)
    … …  
    993968!-----------------------------------------------------------------------
    994969
    995         do ig=1,ngridmx
     970        do ig=1,ngrid
    996971           if (zw2(ig,l+1).gt.1.e-10) then
    997972           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
    … …  
    1016991
    1017992      do l=1,zlmax
    1018          do ig=1,ngridmx
     993         do ig=1,ngrid
    1019994            eee0=entr(ig,l)
    1020995            ddd0=detr(ig,l)
    … …  
    10441019
    10451020! Check again that everything cancels at zmax
    1046       do ig=1,ngridmx
     1021      do ig=1,ngrid
    10471022         fm(ig,lmax(ig)+1)=0.
    10481023         entr(ig,lmax(ig))=0.
    … …  
    10511026
    10521027!-----------------------------------------------------------------------
    1053 ! Summary of the number of modifications that were necessary (if lwrite=.true.
     1028! Summary of the number of modifications that were necessary (if thermverbose=.true.
    10541029! and only if there were a lot of them)
    10551030!-----------------------------------------------------------------------
    10561031
    10571032!IM 090508 beg
    1058       IF (lwrite) THEN
    1059       if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngridmx/4. ) then
     1033      IF (thermverbose) THEN
     1034      if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngrid/4. ) then
    10601035         print*,'thermcell warning : large number of corrections'
    10611036         print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
    … …  
    10891064! Based on equation 14 in appendix 4.2
    10901065
    1091       do ig=1,ngridmx
     1066      do ig=1,ngrid
    10921067         if(lmax(ig) .gt. 1) then
    10931068         do k=1,lmax(ig)
    … …  
    10951070     &   fm(ig,k)*ztv(ig,k)+fm(ig,k)*ztva(ig,k)-fm(ig,k+1)*ztva(ig,k+1))
    10961071            if (ztv(ig,k) + ptimestep*zdthladj(ig,k) .le. 0.) then
    1097       IF (lwrite) THEN
     1072      IF (thermverbose) THEN
    10981073              print*,'Teta<0 in thermcell_dTeta up: qenv .. dq : ', ztv(ig,k),ptimestep*zdthladj(ig,k)
    10991074      ENDIF
    … …  
    11121087      ztvd(:,:)=ztv(:,:)
    11131088      fm_down(:,:)=0.
    1114       do ig=1,ngridmx
     1089      do ig=1,ngrid
    11151090         if (lmax(ig) .gt. 1) then
    11161091         do l=1,lmax(ig)
    … …  
    11391114       zdthladj_down(:,:)=0.
    11401115
    1141       do ig=1,ngridmx
     1116      do ig=1,ngrid
    11421117         if(lmax(ig) .gt. 1) then
    11431118! No downdraft in the very-near surface layer, we begin at k=3
    … …  
    11481123     & fm_down(ig,k)*ztv(ig,k)+fm_down(ig,k)*ztvd(ig,k)-fm_down(ig,k+1)*ztvd(ig,k+1))
    11491124            if (ztv(ig,k) + ptimestep*zdthladj_down(ig,k) .le. 0.) then
    1150       IF (lwrite) THEN
     1125      IF (thermverbose) THEN
    11511126              print*,'q<0 in thermcell_dTeta down: qenv .. dq : ', ztv(ig,k),ptimestep*zdthladj_down(ig,k)
    11521127      ENDIF
    … …  
    11641139      fraca(:,:)=0.
    11651140      do l=2,zlmax
    1166          do ig=1,ngridmx
     1141         do ig=1,ngrid
    11671142            if (zw2(ig,l).gt.1.e-10) then
    11681143            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
    … …  
    11821157      ratiom(:,:)=1.
    11831158
    1184       if (ico2.ne.0) then
     1159      if (igcm_co2.ne.0) then
    11851160      detrmod(:,:)=0.
    11861161      do k=1,zlmax
    1187          do ig=1,ngridmx
     1162         do ig=1,ngrid
    11881163            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
    11891164     &      +entr(ig,k)
    … …  
    11951170      enddo
    11961171
    1197       call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
     1172      call thermcell_dqup(ngrid,nlayer,ptimestep     &
    11981173     &     ,fm,entr,detrmod,  &
    1199      &    masse,pq(:,:,ico2),pdqadj(:,:,ico2),ndt,zlmax)
     1174     &    masse,pq(:,:,igcm_co2),pdqadj(:,:,igcm_co2),ndt,zlmax)
    12001175
    12011176! Compute the ratio between theta and theta_m
    12021177     
    12031178       do l=1,zlmax
    1204           do ig=1,ngridmx
    1205              ratiom(ig,l)=1./(A*(pq(ig,l,ico2)+pdqadj(ig,l,ico2)*ptimestep)+B)
     1179          do ig=1,ngrid
     1180             ratiom(ig,l)=1./(A*(pq(ig,l,igcm_co2)+pdqadj(ig,l,igcm_co2)*ptimestep)+B)
    12061181          enddo
    12071182       enddo
    … …  
    12151190      pdtadj(:,:)=0.
    12161191      do l=1,zlmax
    1217          do ig=1,ngridmx
     1192         do ig=1,ngrid
    12181193         pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
    12191194         enddo
    … …  
    12321207
    12331208   
    1234       do l=1,nlayermx-1
    1235        do ig=1,ngridmx
     1209      do l=1,nlayer-1
     1210       do ig=1,ngrid
    12361211        teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))*ratiom(ig,l)
    12371212        teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))*ratiom(ig,l)
    … …  
    12391214       enddo
    12401215      enddo
    1241       do ig=1,ngridmx
    1242        teta_th_int(ig,nlayermx)=teta_th_int(ig,nlayermx-1)
    1243        teta_env_int(ig,nlayermx)=teta_env_int(ig,nlayermx-1)
    1244        teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
     1216      do ig=1,ngrid
     1217       teta_th_int(ig,nlayer)=teta_th_int(ig,nlayer-1)
     1218       teta_env_int(ig,nlayer)=teta_env_int(ig,nlayer-1)
     1219       teta_down_int(ig,nlayer)=teta_down_int(ig,nlayer-1)
    12451220      enddo
    12461221        heatFlux(:,:)=0.
    … …  
    12491224        heatFlux_down(:,:)=0.
    12501225      do l=1,zlmax
    1251        do ig=1,ngridmx
     1226       do ig=1,ngrid
    12521227        heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
    12531228        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
Note: See TracChangeset for help on using the changeset viewer.