Changeset 1028


Ignore:
Timestamp:
Sep 5, 2013, 12:29:13 PM (11 years ago)
Author:
emillour
Message:

Mars GCM: Cleanup of the thermals routines (these changes in the files do not change GCM results).
AC+EM

Location:
trunk/LMDZ.MARS/libf/phymars
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/phymars/callkeys.h

    r833 r1028  
    1212     &   ,callg2d,linear,rayleigh,tracer,active,doubleq,submicron       &
    1313     &   ,lifting,callddevil,scavenging,sedimentation,activice,water    &
    14      &   ,tifeedback,microphys,caps,photochem,calltherm,outptherm       &
     14     &   ,tifeedback,microphys,caps,photochem,calltherm                 &
    1515     &   ,callrichsl,callslope,tituscap
    1616     
    … …  
    2626     &   ,callnirco2,callnlte,callthermos,callconduct,                  &
    2727     &    calleuv,callmolvis,callmoldiff,thermochem,thermoswater        &
    28      &   ,calltherm,outptherm,callrichsl,callslope,tituscap
     28     &   ,calltherm,callrichsl,callslope,tituscap
    2929
    3030

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

    r709 r1028  
     1! ---------------------------------------------------------------------
     2! Arnaud Colaitis 2011-01-05
    13!
    2 ! AC 2011-01-05
     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.
    311!
     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! ---------------------------------------------------------------------
    417      SUBROUTINE calltherm_interface (firstcall, &
    518     & zzlev,zzlay, &
    … …  
    922     & pdhdif,hfmax,wstar,sensibFlux)
    1023
    11        USE ioipsl_getincom
    12 
     24
     25      USE ioipsl_getincom, only : getin ! to read options from external file "run.def"
    1326      implicit none
    14 #include "callkeys.h"
    15 #include "dimensions.h"
    16 #include "dimphys.h"
    17 #include "comcstfi.h"
    18 #include "tracer.h"
     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
     34#include "comcstfi.h" !contains physical constant values such as
     35                      ! "g" : gravitational acceleration (m.s-2)
     36                      ! "r" : recuced gas constant (J.K-1.mol-1)
     37                      ! "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
    1941
    2042!--------------------------------------------------------
    … …  
    2244!--------------------------------------------------------
    2345
    24 !      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
    25       REAL, INTENT(IN) :: ptimestep
    26       REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
    27       REAL, INTENT(IN) :: pplay(ngridmx,nlayermx)
    28       REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
    29       REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx)
    30       REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)
    31       REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
    32       REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
    33       LOGICAL, INTENT(IN) :: firstcall
    34       REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx)
    35       REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx)
    36       REAL, INTENT(IN) :: pdt(ngridmx,nlayermx)
    37       REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
    38       REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
    39       REAL, INTENT(IN) :: pdhdif(ngridmx,nlayermx)
    40       REAL, INTENT(IN) :: sensibFlux(ngridmx)
     46      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
     57                                                                      ! before thermals du/dt (m/s/s)
     58      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx) ! tracer concentration change from routines called
     59                                                     ! 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,
     63                                                   ! 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
    4166
    4267!--------------------------------------------------------
    … …  
    4469!--------------------------------------------------------
    4570
    46       REAL, INTENT(OUT) :: pdu_th(ngridmx,nlayermx)
    47       REAL, INTENT(OUT) :: pdv_th(ngridmx,nlayermx)
    48       REAL, INTENT(OUT) :: pdt_th(ngridmx,nlayermx)
    49       REAL, INTENT(OUT) :: pdq_th(ngridmx,nlayermx,nqmx)
    50       INTEGER, INTENT(OUT) :: lmax(ngridmx)
    51       REAL, INTENT(OUT) :: zmaxth(ngridmx)
    52       REAL, INTENT(OUT) :: pbl_dtke(ngridmx,nlayermx+1)
    53       REAL, INTENT(OUT) :: wstar(ngridmx)
     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
    5488
    5589!--------------------------------------------------------
    … …  
    6599      REAL zw2(ngridmx,nlayermx+1)
    66100      REAL fraca(ngridmx,nlayermx+1),zfraca(ngridmx,nlayermx+1)
    67       REAL ztla(ngridmx,nlayermx)
    68101      REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
    69102      REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
    70103      REAL lmax_real(ngridmx)
    71104      REAL masse(ngridmx,nlayermx)
    72       LOGICAL qtransport_thermals,dtke_thermals
     105
    73106      INTEGER l,ig,iq,ii(1),k
    74107      CHARACTER (LEN=20) modname
    … …  
    79112
    80113      REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
    81       REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
    82       REAL dq2_the(ngridmx,nlayermx)
    83114      INTEGER isplit
    84115      INTEGER,SAVE :: nsplit_thermals
    … …  
    113144!--------------------------------------------------------
    114145
    115       INTEGER ico2
    116       SAVE ico2
    117 
    118 ! **********************************************************************
    119 ! Initialization
     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
     176
     177! **********************************************************************
     178! Initializations
    120179! **********************************************************************
    121180
    … …  
    128187      q2_therm(:,:)=0.
    129188      dq2_therm(:,:)=0.
    130       ztla(:,:)=0.
    131189      pbl_dtke(:,:)=0.
    132190      fm_therm(:,:)=0.
    … …  
    150208
    151209! **********************************************************************
    152 ! Preparing inputs for the thermals
    153 ! **********************************************************************
    154 
    155        zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep
    156        zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep
    157        zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep
    158 
    159        pq_therm(:,:,:)=0.
    160        qtransport_thermals=.true. !! default setting
    161        !call getin("qtransport_thermals",qtransport_thermals)
     210! Preparing inputs for the thermals: increment tendancies
     211! from other subroutines called before the thermals model
     212! **********************************************************************
     213
     214       zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep ! u-component of wind velocity
     215       zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep ! v-component of wind velocity
     216       zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep ! temperature
     217
     218       pq_therm(:,:,:)=0. ! tracer concentration
    162219
    163220       if(qtransport_thermals) then
    164           if(tracer) then
    165                 pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep
     221          if(tracer) then ! tracer is a logical that is true if tracer must be transported in the GCM physics
     222                pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep ! tracer concentration
    166223          endif
    167224       endif
    168225
    169        dtke_thermals=.false. !! default setting
    170        call getin("dtke_thermals",dtke_thermals)
     226
    171227       IF(dtke_thermals) THEN
    172228          DO l=1,nlayermx
    … …  
    174230          ENDDO
    175231       ENDIF
    176 
    177 ! **********************************************************************
    178 ! Polar night mixing : theta_m
    179 ! **********************************************************************
    180 
    181       if(firstcall) then
    182         ico2=0
    183         if (tracer) then
    184 !     Prepare Special treatment if one of the tracers is CO2 gas
    185            do iq=1,nqmx
    186              if (noms(iq).eq."co2") then
    187                 ico2=iq
    188              end if
    189            enddo
    190         endif
    191       endif !of if firstcall
    192232
    193233
    … …  
    206246                                ! chosen timestep for the radiative transfer.
    207247                                ! It is recommended to run with 96 timestep per day and
    208                                 ! iradia = 1., configuration in which thermals can run
     248                                ! call radiative transfer at each timestep,
     249                                ! configuration in which thermals can run
    209250                                ! very well with a sub-timestep of 10.
    210251         IF (firstcall) THEN
    211252            r_aspect_thermals=1.  ! same value is OK for GCM and mesoscale
     253! Sub Timestep for the mesoscale model:
    212254#ifdef MESOSCALE
    213255            !! valid for timesteps < 200s
    214256            nsplit_thermals=4
    215257#else
     258! Sub Timestep for the GCM:
     259            ! If there is at least 96 timesteps per day in the gcm, subtimestep factor 10
    216260            IF ((ptimestep .le. 3699.*24./96.) .and. (iradia .eq. 1)) THEN
    217261               nsplit_thermals=10
    218             ELSE
     262            ELSE !if there is less, subtimestep factor 35
    219263               nsplit_thermals=35
    220264            ENDIF
    221265#endif
    222             call getin("nsplit_thermals",nsplit_thermals)
    223             call getin("r_aspect_thermals",r_aspect_thermals)
    224266         ENDIF
    225267
    … …  
    228270! **********************************************************************
    229271
    230          zdt=ptimestep/REAL(nsplit_thermals)
    231 
    232          DO isplit=1,nsplit_thermals
     272         zdt=ptimestep/REAL(nsplit_thermals) !subtimestep
     273
     274         DO isplit=1,nsplit_thermals !temporal loop on the subtimestep
    233275
    234276! Initialization of intermediary variables
    235277
    236 !         zfm_therm(:,:)=0. !init is done inside
    237 !         zentr_therm(:,:)=0.
    238 !         zdetr_therm(:,:)=0.
    239 !         zheatFlux(:,:)=0. 
    240 !         zheatFlux_down(:,:)=0.
    241 !         zbuoyancyOut(:,:)=0.
    242 !         zbuoyancyEst(:,:)=0.
    243278         zzw2(:,:)=0.
    244279         zmax(:)=0.
    245280         lmax(:)=0
    246 !         d_t_the(:,:)=0. !init is done inside
    247 
    248 !         d_u_the(:,:)=0. !transported outside
    249 !         d_v_the(:,:)=0.
    250          dq2_the(:,:)=0.
    251 
    252          if (nqmx .ne. 0 .and. ico2 .ne. 0) then
     281
     282         if (nqmx .ne. 0 .and. ico2 .ne. 0) then !initialize co2 tracer tendancy
    253283            d_q_the(:,:,ico2)=0.
    254284         endif
    255285
     286! CALL to main thermal routine
    256287             CALL thermcell_main_mars(zdt  &
    257288     &      ,pplay,pplev,pphi,zzlev,zzlay  &
    258289     &      ,zu,zv,zt,pq_therm,q2_therm  &
    259      &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
     290     &      ,d_t_the,d_q_the  &
    260291     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
    261292     &      ,r_aspect_thermals &
    262293     &      ,zzw2,fraca,zpopsk &
    263      &      ,ztla,zheatFlux,zheatFlux_down &
     294     &      ,zheatFlux,zheatFlux_down &
    264295     &      ,zbuoyancyOut,zbuoyancyEst)
    265296
    266297      fact=1./REAL(nsplit_thermals)
    267298
    268             d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact
    269 !            d_u_the(:,:)=d_u_the(:,:)*ptimestep*fact
    270 !            d_v_the(:,:)=d_v_the(:,:)*ptimestep*fact
    271             dq2_the(:,:)=dq2_the(:,:)*fact
     299! Update thermals tendancies
     300            d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact !temperature
    272301            if (ico2 .ne. 0) then
    273                d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact
     302               d_q_the(:,:,ico2)=d_q_the(:,:,ico2)*ptimestep*fact !co2 mass mixing ratio
    274303            endif
    275 
    276             zmaxth(:)=zmaxth(:)+zmax(:)*fact
    277             lmax_real(:)=lmax_real(:)+float(lmax(:))*fact
    278             fm_therm(:,:)=fm_therm(:,:)  &
     304            zmaxth(:)=zmaxth(:)+zmax(:)*fact !thermals height
     305            lmax_real(:)=lmax_real(:)+float(lmax(:))*fact !thermals height index
     306            fm_therm(:,:)=fm_therm(:,:)  & !upward mass flux
    279307     &      +zfm_therm(:,:)*fact
    280             entr_therm(:,:)=entr_therm(:,:)  &
     308            entr_therm(:,:)=entr_therm(:,:)  & !entrainment mass flux
    281309     &       +zentr_therm(:,:)*fact
    282             detr_therm(:,:)=detr_therm(:,:)  &
     310            detr_therm(:,:)=detr_therm(:,:)  & !detrainment mass flux
    283311     &       +zdetr_therm(:,:)*fact
    284             zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact
    285 
    286             heatFlux(:,:)=heatFlux(:,:) &
     312            zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact !updraft fractional coverage
     313            heatFlux(:,:)=heatFlux(:,:) & !upward heat flux
    287314     &       +zheatFlux(:,:)*fact
    288             heatFlux_down(:,:)=heatFlux_down(:,:) &
     315            heatFlux_down(:,:)=heatFlux_down(:,:) & !downward heat flux
    289316     &       +zheatFlux_down(:,:)*fact
    290             buoyancyOut(:,:)=buoyancyOut(:,:) &
     317            buoyancyOut(:,:)=buoyancyOut(:,:) & !plume final buoyancy
    291318     &       +zbuoyancyOut(:,:)*fact
    292             buoyancyEst(:,:)=buoyancyEst(:,:) &
     319            buoyancyEst(:,:)=buoyancyEst(:,:) & !plume estimated buoyancy used for vertical velocity computation
    293320     &       +zbuoyancyEst(:,:)*fact
    294  
    295 
    296             zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact
    297 
    298 !  accumulation de la tendance
    299 
    300            d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:)
    301 !           d_u_ajs(:,:)=d_u_ajs(:,:)+d_u_the(:,:)
    302 !           d_v_ajs(:,:)=d_v_ajs(:,:)+d_v_the(:,:)
     321            zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact !vertical velocity
     322
     323! Save tendancies
     324           d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) !temperature tendancy (delta T)
    303325            if (ico2 .ne. 0) then
    304                d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2)
     326               d_q_ajs(:,:,ico2)=d_q_ajs(:,:,ico2)+d_q_the(:,:,ico2) !tracer tendancy (delta q)
    305327            endif
    306 !            dq2_therm(:,:)=dq2_therm(:,:)+dq2_the(:,:)
    307 !  incrementation des variables meteo
    308 
    309             zt(:,:) = zt(:,:) + d_t_the(:,:)
    310 !            zu(:,:) = zu(:,:) + d_u_the(:,:)
    311 !            zv(:,:) = zv(:,:) + d_v_the(:,:)
     328
     329!  Increment temperature and co2 concentration for next pass in subtimestep loop
     330            zt(:,:) = zt(:,:) + d_t_the(:,:) !temperature
    312331            if (ico2 .ne. 0) then
    313332             pq_therm(:,:,ico2) = &
    314      &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2)
     333     &          pq_therm(:,:,ico2) + d_q_the(:,:,ico2) !co2 tracer
    315334            endif
    316 !            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
    317335
    318336
    319337         ENDDO ! isplit
    320338!****************************************************************
     339
     340!     Do we have thermals that are too high ?
    321341
    322342      lmax(:)=nint(lmax_real(:))
    … …  
    329349
    330350! Now that we have computed total entrainment and detrainment, we can
    331 ! advect u, v, and q in thermals. (theta already advected). We can do
    332 ! that separatly because u,v,and q are not used in thermcell_main for
     351! advect u, v, and q in thermals. (potential temperature and co2 MMR
     352! have already been advected in thermcell_main because they are coupled
     353! to the determination of the thermals caracteristics). This is done
     354! separatly because u,v, and q are not used in thermcell_main for
    333355! any thermals-related computation : they are purely passive.
    334356
    … …  
    338360      enddo
    339361
     362! recompute detrainment mass flux from entrainment and upward mass flux
     363! this ensure mass flux conservation
    340364      detrmod(:,:)=0.
    341365      do l=1,zlmax
    … …  
    349373         enddo
    350374      enddo
     375
     376      ! u component of wind velocity advection in thermals
     377      ! result is a derivative (d_u_ajs in m/s/s)
    351378      ndt=10
    352379      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
    … …  
    354381     &     masse,zu,d_u_ajs,ndt,zlmax)
    355382
     383      ! v component of wind velocity advection in thermals
     384      ! result is a derivative (d_v_ajs in m/s/s)
    356385      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
    357386     &       ,fm_therm,entr_therm,detrmod,  &
    358387     &     masse,zv,d_v_ajs,ndt,zlmax)
     388
     389      ! non co2 tracers advection in thermals
     390      ! result is a derivative (d_q_ajs in kg/kg/s)
    359391
    360392      if (nqmx .ne. 0.) then
    … …  
    368400      endif
    369401
     402      ! tke advection in thermals
     403      ! result is a tendancy (d_u_ajs in J)
    370404      if (dtke_thermals) then
    371       detrmod(:,:)=0.
    372       ndt=10
    373       do l=1,zlmax
    374          do ig=1,ngridmx
    375             detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) &
    376      &      +entr_therm(ig,l)
    377             if (detrmod(ig,l).lt.0.) then
    378                entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l)
    379                detrmod(ig,l)=0.
    380             endif
    381          enddo
    382       enddo
    383405      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    384406     &     ,fm_therm,entr_therm,detrmod,  &
    … …  
    386408      endif
    387409
     410      ! compute wmax for diagnostics
    388411      DO ig=1,ngridmx
    389412         wmax(ig)=MAXVAL(zw2(ig,:))
    … …  
    408431      enddo
    409432
     433           ! if tracers are transported in thermals, update output variables, else these are 0.
    410434           if(qtransport_thermals) then
    411435              if(tracer) then
    … …  
    413437                if (iq .ne. ico2) then
    414438                  do l=1,zlmax
    415                      pdq_th(:,l,iq)=d_q_ajs(:,l,iq)
     439                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq) !non-co2 tracers d_q_ajs are dq/dt (kg/kg/s)
    416440                  enddo
    417441                else
    418442                  do l=1,zlmax
    419                      pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep
     443                     pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep !co2 tracer d_q_ajs is dq (kg/kg)
    420444                  enddo
    421445                endif
    … …  
    424448           endif
    425449
     450           ! if tke is transported in thermals, update output variable, else this is 0.
    426451           IF(dtke_thermals) THEN
    427452              DO l=2,nlayermx
    … …  
    433458           ENDIF
    434459
     460           ! update output variable for temperature. d_t_ajs is delta T in (K), pdt_th is dT/dt in (K/s)
    435461           do l=1,zlmax
    436462              pdt_th(:,l)=d_t_ajs(:,l)/ptimestep
    … …  
    439465
    440466! **********************************************************************
    441 ! Compute the free convection velocity scale for vdifc
    442 ! **********************************************************************
    443 
     467! SURFACE LAYER INTERFACE
     468! Compute the free convection velocity w* scale for surface layer gustiness
     469! speed parameterization. The computed value of w* will be used at the next
     470! timestep to modify surface-atmosphere exchange fluxes, because the surface
     471! layer scheme and diffusion are called BEFORE the thermals. (outside of these
     472! routines)
     473! **********************************************************************
    444474
    445475! Potential temperature gradient
    … …  
    453483      ENDDO
    454484
    455 ! Computation of the pbl mixed layer temperature
     485! Computation of the PBL mixed layer temperature
    456486
    457487      DO ig=1, ngridmx
    … …  
    460490      ENDDO
    461491
    462 ! we must add the heat flux from the diffusion scheme to hfmax
     492! In order to have an accurate w*, we must add the heat flux from the
     493! diffusion scheme to the computation of the maximum heat flux hfmax
     494! Here pdhdif is the potential temperature change from the diffusion
     495! scheme (Mellor and Yamada, see paper section 6, paragraph 57)
    463496
    464497! compute rho as it is after the diffusion
    … …  
    482515      ENDDO
    483516
    484 ! compute w'teta' from diffusion
     517! compute w'theta' (vertical turbulent flux of temperature) from
     518! the diffusion scheme
    485519
    486520      wtdif(:,:)=rpdhd(:,:)/rho(:,:)
    487521
     522! Now we compute the contribution of the thermals to w'theta':
    488523! compute rho as it is after the thermals
    489524
    490525      rho(:,:)=pplay(:,:)/(r*(zt(:,:)))
     526
    491527! integrate -rho*pdhdif
    492528
    … …  
    508544      wtth(:,:)=rpdhd(:,:)/rho(:,:)
    509545
    510 ! We get the max heat flux from thermals and add the contribution from the diffusion
     546! Add vertical turbulent heat fluxes from the thermals and the diffusion scheme
     547! and compute the maximum
    511548
    512549      DO ig=1,ngridmx
    513550        hfmax(ig)=MAXVAL(wtth(ig,:)+wtdif(ig,:))
    514551      ENDDO
     552
     553! Finally we can compute the free convection velocity scale
    515554! We follow Spiga et. al 2010 (QJRMS)
    516555! ------------
    … …  
    524563      ENDDO
    525564
    526 
    527 
    528565! **********************************************************************
    529566! Diagnostics
    530567! **********************************************************************
     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
    531570       
    532         if(outptherm) then
    533         if (ngridmx .eq. 1) then
     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
    534575        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
    535576     &                       'kg/m-2',1,entr_therm)
    … …  
    565606     &         'heat flux PBL','K.m/s',1,wtdif(:,:)+wtth(:,:))
    566607
    567       else
     608      else !these are 3D outputs
    568609
    569610        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
    … …  
    583624
    584625      endif
    585       endif
     626      endif ! of if (outptherm)
    586627
    587628       END

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

    r833 r1028  
    267267         write(*,*) " calltherm = ",calltherm
    268268
    269          write(*,*) "output thermal diagnostics ?"
    270          outptherm=.false. ! default value
    271          call getin("outptherm",outptherm)
    272          write(*,*) " outptherm = ",outptherm
    273 
    274269         write(*,*) "call convective adjustment ?"
    275270         calladj=.true. ! default value

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

    r628 r1028  
     1! ---------------------------------------------------------------------
     2! Arnaud Colaitis 2011-01-05
     3! --------------------------------------------------------------------
    14      subroutine thermcell_dqup(ngrid,nlayer,ptimestep,fm,entr,detr,  &
    25     &    masse0,q_therm,dq_therm,ndt,zlmax)
    … …  
    58!=======================================================================
    69!
    7 !   Calcul du transport verticale dans la couche limite en presence
    8 !   de "thermiques" explicitement representes
    9 !   calcul du dq/dt une fois qu'on connait les ascendances
    10 !   Version modifiee pour prendre les downdrafts a la place de la
    11 !   subsidence compensatoire
     10!   Compute the thermals contribution of explicit thermals
     11!   to vertical transport in the PBL.
     12!   dq is computed once upward, entrainment and detrainment mass fluxes
     13!   are known.
    1214!
    1315!   Version with sub-timestep for Martian thin layers
    … …  
    1517!=======================================================================
    1618
    17 #include "dimensions.h"
    18 #include "dimphys.h"
     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
    1923
    2024! ============================ INPUTS ============================
    2125
    22       INTEGER, INTENT(IN) :: ngrid,nlayer
    23       REAL, INTENT(IN) :: ptimestep
    24       REAL, INTENT(IN) :: fm(ngridmx,nlayermx+1)
    25       REAL, INTENT(IN) :: entr(ngridmx,nlayermx)
    26       REAL, INTENT(IN) :: detr(ngridmx,nlayermx)
    27       REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
    28       REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
    29       INTEGER, INTENT(IN) :: ndt
    30       INTEGER, INTENT(IN) :: zlmax
     26      INTEGER, INTENT(IN) :: ngrid,nlayer ! number of grid points and number of levels
     27      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
     33      INTEGER, INTENT(IN) :: ndt ! number of subtimesteps
     34      INTEGER, INTENT(IN) :: zlmax ! index of maximum layer
    3135
    3236! ============================ OUTPUTS ===========================
    … …  
    4448! =========== Init ==============================================
    4549
    46       qa(:,:)=q_therm(:,:)
    47       q(:,:)=q_therm(:,:)
     50      qa(:,:)=q_therm(:,:) !q profile in the updraft
     51      q(:,:)=q_therm(:,:) !mean q profile
    4852
    4953! ====== Computing q ============================================
     54! Based on equation 14 in appendix 4.2
    5055
    5156      dq_therm(:,:)=0.
    … …  
    5358      invflux0(:,:)=ztimestep/masse0(:,:)     
    5459
    55       do i=1,ndt
     60      do i=1,ndt !subtimestep loop
    5661
    57       do ig=1,ngrid
    58          qa(ig,1)=q(ig,1)
    59       enddo
     62        do ig=1,ngrid
     63           qa(ig,1)=q(ig,1)
     64       enddo
    6065
    61       do k=2,zlmax
    62          do ig=1,ngridmx
    63             if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
    64      &         1.e-5*masse0(ig,k)) then
    65          qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
    66      &     /(fm(ig,k+1)+detr(ig,k))
    67             else
    68                qa(ig,k)=q(ig,k)
    69             endif
    70          enddo
    71       enddo
     66        do k=2,zlmax
     67           do ig=1,ngridmx
     68              if ((fm(ig,k+1)+detr(ig,k))*ptimestep.gt.  &
     69     &        1.e-5*masse0(ig,k)) then
     70                 qa(ig,k)=(fm(ig,k)*qa(ig,k-1)+entr(ig,k)*q(ig,k))  &
     71     &           /(fm(ig,k+1)+detr(ig,k))
     72              else
     73                 qa(ig,k)=q(ig,k)
     74              endif
     75           enddo
     76        enddo
    7277
    73       do k=1,zlmax
    74            q(:,k)=q(:,k)+         &
    75      & (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
     78        do k=1,zlmax
     79          q(:,k)=q(:,k)+         &
     80     &    (detr(:,k)*qa(:,k)-entr(:,k)*q(:,k) &
    7681     &    -fm(:,k)*q(:,k)+fm(:,k+1)*q(:,k+1))  &
    77      &               *invflux0(:,k)
    78       enddo
     82     &    *invflux0(:,k)
     83        enddo
    7984
    8085      enddo !of do i=1,ndt
    8186
    8287! ====== Derivative ==============================================
    83 
    8488
    8589         do k=1,zlmax

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

    r682 r1028  
     1!=======================================================================
     2! THERMCELL_MAIN_MARS
     3! Author : A Colaitis
     4!
     5! This routine is called by calltherm_interface and is inside a sub-timestep
     6! loop. It computes thermals properties from parametrized entrainment and
     7! detrainment rate as well as the source profile.
     8! 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
     12!=======================================================================
    113!
    214!
    … …  
    416     &                  ,pplay,pplev,pphi,zlev,zlay  &
    517     &                  ,pu,pv,pt,pq,pq2  &
    6      &                  ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj  &
     18     &                  ,pdtadj,pdqadj  &
    719     &                  ,fm,entr,detr,lmax,zmax  &
    820     &                  ,r_aspect &
    921     &                  ,zw2,fraca &
    10      &                  ,zpopsk,ztla,heatFlux,heatFlux_down &
     22     &                  ,zpopsk,heatFlux,heatFlux_down &
    1123     &                  ,buoyancyOut, buoyancyEst)
    1224
    … …  
    1426
    1527!=======================================================================
    16 ! Mars version of thermcell_main. Author : A Colaitis
    17 !=======================================================================
    18 
    19 #include "dimensions.h"
    20 #include "dimphys.h"
    21 #include "comcstfi.h"
    22 #include "tracer.h"
    23 #include "callkeys.h"
    24 
     28
     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
     36#include "comcstfi.h" !contains physical constant values such as
     37                      ! "g" : gravitational acceleration (m.s-2)
     38                      ! "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
    2542! ============== INPUTS ==============
    2643
    27       REAL, INTENT(IN) :: ptimestep,r_aspect
    28       REAL, INTENT(IN) :: pt(ngridmx,nlayermx)
    29       REAL, INTENT(IN) :: pu(ngridmx,nlayermx)
    30       REAL, INTENT(IN) :: pv(ngridmx,nlayermx)
    31       REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx)
    32       REAL, INTENT(IN) :: pq2(ngridmx,nlayermx)
    33       REAL, INTENT(IN) :: pplay(ngridmx,nlayermx)
    34       REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
    35       REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
    36       REAL, INTENT(IN) :: zlay(ngridmx,nlayermx)
    37       REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1)
     44      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
    3856
    3957! ============== OUTPUTS ==============
    4058
    41       REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx)
    42       REAL :: pduadj(ngridmx,nlayermx)
    43       REAL :: pdvadj(ngridmx,nlayermx)
    44       REAL :: pdqadj(ngridmx,nlayermx,nqmx)
    45 !      REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx)
    46       REAL :: pdq2adj(ngridmx,nlayermx)
    47       REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1)
    48 
    49 ! Diagnostics
     59! TEMPERATURE
     60      REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx) !temperature change from thermals dT/dt (K/s)
     61
     62! DIAGNOSTICS
     63      REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1) ! vertical velocity (m/s)
    5064      REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
    51      REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
    52 !      REAL, INTENT(OUT) :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
    53 !      REAL, INTENT(OUT) :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
     65      REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
     66
     67! ============== 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)
    5469
    5570! dummy variables when output not needed :
    5671
    57 !      REAL :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
    58 !      REAL :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
    5972      REAL :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
    6073      REAL :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
    6174
    62 
    6375! ============== LOCAL ================
    6476
    6577      INTEGER ig,k,l,ll,iq
    6678      INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
    67       REAL linter(ngridmx)
    6879      REAL zmax(ngridmx)
    6980      REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx)
    70       REAL zmax_sec(ngridmx)
    7181      REAL zh(ngridmx,nlayermx)
    7282      REAL zdthladj(ngridmx,nlayermx)
    … …  
    8595
    8696      REAL wmax(ngridmx)
    87       REAL wmax_sec(ngridmx)
    8897      REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
    8998
    … …  
    115124      REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
    116125      LOGICAL activecell(ngridmx),activetmp(ngridmx)
    117       REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega,adalim
     126      REAL a1,b1,ae,be,ad,bd,fdfu,b1inv,a1inv,omega
    118127      INTEGER tic
    119128
    … …  
    145154      REAL f_old,ddd0,eee0,ddd,eee,zzz
    146155      REAL fomass_max,alphamax
    147 
    148 ! =========================================
    149 
    150 ! ============= DTETA VARIABLES ===========
    151 
    152 ! rien : on prends la divergence du flux turbulent
    153 
    154 ! =========================================
    155 
    156 ! ============= DV2 VARIABLES =============
    157 !               not used for now
    158 
    159       REAL qa(ngridmx,nlayermx),detr_dv2(ngridmx,nlayermx),zf,zf2
    160       REAL wvd(ngridmx,nlayermx+1),wud(ngridmx,nlayermx+1)
    161       REAL gamma0(ngridmx,nlayermx+1),gamma(ngridmx,nlayermx+1)
    162       REAL ue(ngridmx,nlayermx),ve(ngridmx,nlayermx)
    163       LOGICAL ltherm(ngridmx,nlayermx)
    164       REAL dua(ngridmx,nlayermx),dva(ngridmx,nlayermx)
    165       INTEGER iter
    166       INTEGER nlarga0
    167156
    168157! =========================================
    … …  
    183172
    184173!-----------------------------------------------------------------------
    185 !   initialisation:
     174!   initialization:
    186175!   ---------------
    187176
    188       entr(:,:)=0.
    189       detr(:,:)=0.
    190       fm(:,:)=0.
    191 !      zu(:,:)=pu(:,:)
    192 !      zv(:,:)=pv(:,:)
    193       zhc(:,:)=pt(:,:)/zpopsk(:,:)
     177      entr(:,:)=0. ! entrainment mass flux
     178      detr(:,:)=0. ! detrainment mass flux
     179      fm(:,:)=0. ! upward mass flux
     180      zhc(:,:)=pt(:,:)/zpopsk(:,:) ! potential temperature
    194181      ndt=1
    195182
    196183! **********************************************************************
    197 ! Taking into account vertical molar mass gradients
     184! In order to take into account the effect of vertical molar mass
     185! gradient on convection, we define modified theta that depends
     186! on the mass mixing ratio of Co2 in the cell.
     187! See for details:
     188!
     189! Forget, F. and Millour, E. et al. "Non condensable gas enrichment and depletion
     190! in the martian polar regions", third international workshop on the Mars Atmosphere:
     191! Modeling and Observations, 1447, 9106. year: 2008
     192!
     193! This is especially important for modelling polar convection.
    198194! **********************************************************************
     195
     196     
     197!.......................................................................
     198!  Special treatment for co2 at firstcall: compute coefficients and
     199!  get index of co2 tracer
     200!.......................................................................
    199201
    200202      if(firstcall) then
    … …  
    235237       end if
    236238
    237 
     239!------------------------------------------------------------------------
     240! where are the different quantities defined ?
    238241!------------------------------------------------------------------------
    239242!                       --------------------
    … …  
    255258
    256259!-----------------------------------------------------------------------
    257 !   Calcul des altitudes des couches
     260!   Densities at layer and layer interface (see above), mass:
    258261!-----------------------------------------------------------------------
    259262
    260 !      do l=2,nlayermx
    261 !         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g
    262 !      enddo
    263 !         zlev(:,1)=0.
    264 !         zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g
    265 
    266 !         zlay(:,:)=pphi(:,:)/g
    267 !-----------------------------------------------------------------------
    268 !   Calcul des densites
    269 !-----------------------------------------------------------------------
    270 
    271263      rho(:,:)=pplay(:,:)/(r*pt(:,:))
    272264
    … …  
    274266
    275267      do l=2,nlayermx
    276 !         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
    277268          rhobarz(:,l)=pplev(:,l)/(r*0.5*(pt(:,l)+pt(:,l-1)))
    278269      enddo
    279270
    280 !calcul de la masse
     271! mass computation
    281272      do l=1,nlayermx
    282273         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
    … …  
    284275
    285276
     277!-----------------------------------------------------------------
     278!   Schematic representation of an updraft:
    286279!------------------------------------------------------------------
    287280!
    … …  
    330323
    331324!=============================================================================
    332 !  Calculs initiaux ne faisant pas intervenir les changements de phase
     325! Mars version: no phase change is considered, we use a "dry" definition
     326! for the potential temperature.
    333327!=============================================================================
    334328
    335329!------------------------------------------------------------------
    336 !  1. alim_star est le profil vertical de l'alimentation a la base du
    337 !     panache thermique, calcule a partir de la flotabilite de l'air sec
    338 !  2. lmin et lalim sont les indices inferieurs et superieurs de alim_star
     330!  1. alim_star is the source layer vertical profile in the lowest layers
     331! of the thermal plume. Computed from the air buoyancy
     332!  2. lmin and lalim are the indices of begining and end of source profile
    339333!------------------------------------------------------------------
    340334!
    … …  
    343337
    344338!-----------------------------------------------------------------------------
    345 !  3. wmax_sec et zmax_sec sont les vitesses et altitudes maximum d'un
    346 !     panache sec conservatif (e=d=0) alimente selon alim_star
    347 !     Il s'agit d'un calcul de type CAPE
    348 !     zmax_sec est utilise pour determiner la geometrie du thermique.
    349 !------------------------------------------------------------------------------
    350 !---------------------------------------------------------------------------------
    351 !calcul du melange et des variables dans le thermique
    352 !--------------------------------------------------------------------------------
     339!  3. wmax and zmax are maximum vertical velocity and altitude of a
     340!     conservative plume (entrainment = detrainment = 0) using only
     341!     the source layer. This is a CAPE computation used for determining
     342!     the closure mass flux.
     343!-----------------------------------------------------------------------------
    353344
    354345! ===========================================================================
    … …  
    356347! ===========================================================================
    357348
    358 ! Initialisations des variables reeles
    359       ztva(:,:)=ztv(:,:)
    360       ztva_est(:,:)=ztva(:,:)
    361       ztla(:,:)=0.
    362       zdz=0.
    363       zbuoy(:,:)=0.
    364       w_est(:,:)=0.
    365       f_star(:,:)=0.
    366       wa_moy(:,:)=0.
    367       linter(:)=1.
     349! Initialization
     350      ztva(:,:)=ztv(:,:) ! temperature in the updraft = temperature of the env.
     351      ztva_est(:,:)=ztva(:,:) ! estimated temp. in the updraft
     352      ztla(:,:)=0. !intermediary variable
     353      zdz=0. !layer thickness
     354      zbuoy(:,:)=0. !buoyancy
     355      w_est(:,:)=0. !estimated vertical velocity
     356      f_star(:,:)=0. !non-dimensional upward mass flux f*
     357      wa_moy(:,:)=0. !vertical velocity
    368358
    369359! --------------------------------------------------------------------------
    370360! -------------- MAIN PARAMETERS FOR THERMALS MODEL ------------------------
    371 ! --------------  see thermiques.pro and getfit.py -------------------------
    372 
    373 !      a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6  ! svn baseline
    374 
    375 ! Using broad downdraft selection
    376 !      a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57
    377 !      ad = 0.0005114  ; bd = -0.662
    378 !      fdfu = -1.9
    379 
    380 ! Using conditional sampling downdraft selection
    381       a1=1.4716 ; b1=0.0005698 ; ae=0.03683 ; be = 0.57421
    382       ad = 0.00048088  ; bd = -0.6697
    383 !      fdfu = -1.3
    384       fdfu=-0.8
    385       a1inv=a1
    386       b1inv=b1
    387       omega=0.
    388       adalim=0.
    389 
    390 ! One good config for 34/35 levels
    391 !      a1inv=a1
    392 !      b1inv=b1
    393 !      be=1.1*be
    394 
    395 ! Best configuration for 222 levels:
    396 
    397 !      omega=0.06
    398 !      b1=0.
    399 !      a1=1.
    400 !      a1inv=0.25*a1
    401 !      b1inv=0.0002
    402 !!
    403 !!
    404 !!      ae=0.9*ae
    405 
    406 ! Best config for norad 222 levels:
    407 ! and more specifically to C.large :
    408 
    409 !       omega=0.06
    410 !       omega=0.
    411        omega=-0.03
    412 !       omega=0.
    413        a1=1.
    414 !       b1=0.
    415        b1=0.0001
    416        a1inv=a1
    417        be=1.1*be
    418        ad = 0.0004
    419 !       ad=0.0003
    420 !       adalim=0.
    421 
    422 !       b1inv=0.00025
    423        b1inv=b1
    424 
    425 !       b1inv = 0.0003
    426 
    427 !      omega=0.06
    428 ! Trying stuff :
    429 
    430 !      ad=0.00035
    431 !      ae=0.95*ae
    432 !      b1=0.00055
    433 !      omega=0.04
    434 !
    435 !      ad = 0.0003
    436 !      ae=0.9*ae
    437 
    438 !      omega=0.04
    439 !!      b1=0.
    440 !      a1=1.
    441 !      a1inv=a1
    442 !      b1inv=0.0005689
    443 !!      be=1.1*be
    444 !!      ae=0.96*ae
    445 
    446 
    447 !       omega=0.06
    448 !       a1=1.
    449 !       b1=0.
    450 !       be=be
    451 !       a1inv=0.25*a1
    452 !       b1inv=0.0002   
    453 !       ad=1.1*ad
    454 !       ae=1.*ae
     361
     362! Detrainment
     363      ad = 0.0004 !D_2 in paper, see paragraph 44
     364      bd = -0.6697 !D_1 in paper, see paragraph 44
     365
     366! Entrainment
     367      ae = 0.03683 !E_1 in paper, see paragraph 43
     368      be = 0.631631 !E_2 in paper, see paragraph 43
     369
     370! Downdraft
     371      fdfu=-0.8 !downdraft to updraft mass flux ratio, see paper paragraph 48
     372      omega=-0.03 !see paper paragraph 48
     373
     374! Vertical velocity equation
     375      a1=1. !a in paper, see paragraph 41
     376      b1=0.0001 !b in paper, see paragraph 41
     377
     378! Inversion layer
     379      a1inv=a1 !a1 in inversion layer
     380      b1inv=b1 !b1 in inversion layer
    455381! --------------------------------------------------------------------------
    456382! --------------------------------------------------------------------------
    457383! --------------------------------------------------------------------------
    458384
    459 ! Initialisation des variables entieres
     385! Some more initializations
    460386      wmaxa(:)=0.
    461387      lalim(:)=1
    462388
    463389!-------------------------------------------------------------------------
    464 ! On ne considere comme actif que les colonnes dont les deux premieres
    465 ! couches sont instables.
     390! We consider as an activecell columns where the two first layers are
     391! convectively unstable
     392! When it is the case, we compute the source layer profile (alim_star)
     393! see paper appendix 4.1 for details on the source layer
    466394!-------------------------------------------------------------------------
     395
    467396      activecell(:)=ztv(:,1)>ztv(:,2)
    468397          do ig=1,ngridmx
    469398            if (ztv(ig,1)>=(ztv(ig,2))) then
    470399               alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
    471 !     &                       *log(1.+zlev(ig,2))
    472400     &                       *sqrt(zlev(ig,2))
    473 !     &                       *sqrt(sqrt(zlev(ig,2)))
    474 !     &                       /sqrt(zlev(ig,2))
    475 !      &                       *zlev(ig,2)
    476 !      &                     *exp(-zlev(ig,2)/1000.)
    477401               lalim(ig)=2
    478402               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,1)
    … …  
    481405
    482406       do l=2,nlayermx-1
    483 !        do l=2,4
    484          do ig=1,ngridmx
    485            if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1).ne. 0.)) then ! .and. (zlev(ig,l+1) .lt. 1000.)) then
     407         do ig=1,ngridmx
     408           if (ztv(ig,l)>(ztv(ig,l+1)) .and. ztv(ig,1)>=ztv(ig,l) &
     409     &             .and. (alim_star(ig,l-1).ne. 0.)) then
    486410               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    487 !     &                       *log(1.+zlev(ig,l+1))
    488411     &                       *sqrt(zlev(ig,l+1))
    489 !     &                       *sqrt(sqrt(zlev(ig,l+1)))
    490 !     &                       /sqrt(zlev(ig,l+1))
    491 !      &                       *zlev(ig,l+1)
    492 !      &                     *exp(-zlev(ig,l+1)/1000.)
    493412                lalim(ig)=l+1
    494413               alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    … …  
    505424
    506425      alim_star_tot(:)=1.
    507 !      if(alim_star(1,1) .ne. 0.) then
    508 !      print*, 'lalim star' ,lalim(1)
    509 !      endif
    510 
    511 !------------------------------------------------------------------------------
    512 ! Calcul dans la premiere couche
    513 ! On decide dans cette version que le thermique n'est actif que si la premiere
    514 ! couche est instable.
    515 ! Pourrait etre change si on veut que le thermiques puisse se déclencher
    516 ! dans une couche l>1
    517 !------------------------------------------------------------------------------
    518 
    519       do ig=1,ngridmx
    520 ! Le panache va prendre au debut les caracteristiques de l'air contenu
    521 ! dans cette couche.
     426
     427! We compute the initial squared velocity (zw2) and non-dimensional upward mass flux
     428! (f_star) in the first and second layer from the source profile.
     429
     430      do ig=1,ngridmx
    522431          if (activecell(ig)) then
    523432          ztla(ig,1)=ztv(ig,1)
    524 !cr: attention, prise en compte de f*(1)=1 => AC : what ? f*(1) =0. ! (d'ou f*(2)=a*(1)
    525 ! dans un panache conservatif
    526433          f_star(ig,1)=0.
    527          
    528434          f_star(ig,2)=alim_star(ig,1)
    529 
    530435          zw2(ig,2)=2.*g*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    531436      &                     *(zlev(ig,2)-zlev(ig,1))  &
    532       &                    *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
     437      &                    *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1)) !0.4=von Karman constant
    533438          w_est(ig,2)=zw2(ig,2)
    534 
    535439          endif
    536440      enddo
    537441
    538 
    539442!==============================================================================
    540 !boucle de calcul de la vitesse verticale dans le thermique
     443!==============================================================================
     444!==============================================================================
     445! LOOP ON VERTICAL LEVELS
    541446!==============================================================================
    542447      do l=2,nlayermx-1
    543448!==============================================================================
    544 
    545 
    546 ! On decide si le thermique est encore actif ou non
    547 ! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
     449!==============================================================================
     450!==============================================================================
     451
     452
     453! is the thermal plume still active ?
    548454          do ig=1,ngridmx
    549455             activecell(ig)=activecell(ig) &
    … …  
    553459
    554460!---------------------------------------------------------------------------
    555 ! calcul des proprietes thermodynamiques et de la vitesse de la couche l
    556 ! sans tenir compte du detrainement et de l'entrainement dans cette
    557 ! couche
    558 ! C'est a dire qu'on suppose
    559 ! ztla(l)=ztla(l-1)
    560 ! Ici encore, on doit pouvoir ajouter entr_star (qui peut etre calculer
    561 ! avant) a l'alimentation pour avoir un calcul plus propre
     461!
     462! .I. INITIALIZATION
     463!
     464! Computations of the temperature and buoyancy properties in layer l,
     465! without accounting for entrainment and detrainment. We are therefore
     466! assuming constant temperature in the updraft
     467!
     468! This computation yields an estimation of the buoyancy (zbuoy) and thereforce
     469! an estimation of the velocity squared (w_est)
    562470!---------------------------------------------------------------------------
    563471
    564472          do ig=1,ngridmx
    565473             if(activecell(ig)) then
    566 !                if(l .lt. lalim(ig)) then
    567 !               ztva_est(ig,l)=(f_star(ig,l)*ztla(ig,l-1)+  &
    568 !     &            alim_star(ig,l)*ztv(ig,l))  &
    569 !     &            /(f_star(ig,l)+alim_star(ig,l))
    570 !                else
    571474                ztva_est(ig,l)=ztla(ig,l-1)
    572 !                endif
    573475
    574476                zdz=zlev(ig,l+1)-zlev(ig,l)
    575477                zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     478
     479                ! Estimated vertical velocity squared
     480                ! (discretized version of equation 12 in paragraph 40 of paper)
    576481
    577482                if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
    … …  
    588493
    589494!-------------------------------------------------
    590 !calcul des taux d'entrainement et de detrainement
     495! Compute corresponding non-dimensional (ND) entrainment and detrainment rates
    591496!-------------------------------------------------
    592497
    … …  
    597502
    598503          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
     504
     505         ! ND entrainment rate, see equation 16 of paper (paragraph 43)
     506
    599507          entr_star(ig,l)=f_star(ig,l)*zdz*  &
    600508        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
    601 !         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
    602509          else
    603510          entr_star(ig,l)=0.
    … …  
    607514             if(l .lt. lalim(ig)) then
    608515
    609 !                detr_star(ig,l)=0.
    610                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
    611             &  adalim
     516                detr_star(ig,l)=0.
    612517             else
    613518
    614 !                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    615 !              &     0.0105*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.7)
    616 !                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    617 !              &     0.0085*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.55)
    618 
    619 ! last baseline from direct les
    620 !                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    621 !               &     0.065*(2.5*(zbuoy(ig,l)/zw2m))**0.75
    622 
    623 ! new param from continuity eq with a fit on dfdz
     519         ! ND detrainment rate, see paragraph 44 of paper
     520
    624521                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    625522            &  ad
    626 !             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
    627 
    628 !               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
    629 !                &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0008)     
    630 
    631 !              &     0.014*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.35)
    632 !                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    633 !              &     ((zbuoy(ig,l)/zw2m)/2.222)! + 0.0002)
    634523
    635524             endif
    … …  
    637526          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    638527                &     MAX(ad,bd*zbuoy(ig,l)/zw2m)
    639 !             &  Max(0., 0.001 - 0.45*zbuoy(ig,l)/zw2m)
    640 !             &  Max(0., Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
    641 
    642 
    643 !               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)      !svn baseline
    644 
    645 !              &  *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m)
    646 !              &  *5.*(-afact*zbuoy(ig,l)/zw2m)
    647 
    648 ! last baseline from direct les
    649 !               &     0.065*(-2.5*(zbuoy(ig,l)/zw2m))**0.75
    650 
    651 ! new param from continuity eq with a fit on dfdz
    652 
    653528
    654529          endif
    655530
    656 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
    657 ! alim_star et 0 sinon
     531! If we are still in the source layer, we define the source layer entr. rate  (alim_star) as the
     532! maximum between the source entrainment rate and the estimated entrainment rate.
    658533
    659534       if (l.lt.lalim(ig)) then
    … …  
    662537       endif
    663538
    664 ! Calcul du flux montant normalise
     539! Compute the non-dimensional upward mass flux at layer l+1
     540! using equation 11 of appendix 4.2 in paper
    665541
    666542      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    … …  
    670546      enddo
    671547
    672 
    673 !----------------------------------------------------------------------------
    674 !calcul de la vitesse verticale en melangeant Tl et qt du thermique
    675 !---------------------------------------------------------------------------
    676 
     548! -----------------------------------------------------------------------------------
     549!
     550! .II. CONVERGENCE LOOP
     551!
     552! We have estimated a vertical velocity profile and refined the source layer profile
     553! We now conduct iterations to compute:
     554!
     555! - the temperature inside the updraft from the estimated entrainment/source, detrainment,
     556! and upward mass flux.
     557! - the buoyancy from the new temperature inside the updraft
     558! - the vertical velocity from the new buoyancy
     559! - the entr., detr. and upward mass flux from the new buoyancy and vertical velocity
     560!
     561! This loop (tic) converges quickly. We have hardcoded 6 iterations from empirical observations.
     562! Convergence occurs in 1 or 2 iterations in most cases.
     563! -----------------------------------------------------------------------------------
     564
     565! -----------------------------------------------------------------------------------
     566! -----------------------------------------------------------------------------------
    677567      DO tic=0,5  ! internal convergence loop
     568! -----------------------------------------------------------------------------------
     569! -----------------------------------------------------------------------------------
     570
     571! Is the cell still active ?
    678572      activetmp(:)=activecell(:) .and. f_star(:,l+1)>1.e-10
     573
     574! If the cell is active, compute temperature inside updraft
    679575      do ig=1,ngridmx
    680576       if (activetmp(ig)) then
    … …  
    687583      enddo
    688584
     585! Is the cell still active with respect to temperature variations ?
    689586      activetmp(:)=activetmp(:).and.(abs(ztla(:,l)-ztva(:,l)).gt.0.01)
    690587
     588! Compute new buoyancu and vertical velocity
    691589      do ig=1,ngridmx
    692590      if (activetmp(ig)) then
    … …  
    694592           zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    695593
     594          ! (discretized version of equation 12 in paragraph 40 of paper)
    696595           if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
    697596           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-         &
    … …  
    704603
    705604! ================ RECOMPUTE ENTR, DETR, and F FROM NEW W2 ===================
     605         ! ND entrainment rate, see equation 16 of paper (paragraph 43)
     606         ! ND detrainment rate, see paragraph 44 of paper
    706607
    707608      do ig=1,ngridmx
    … …  
    713614          entr_star(ig,l)=f_star(ig,l)*zdz*  &
    714615        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
    715 !         &  MAX(0.,log(1. + 0.03*sqrt(a1*(zbuoy(ig,l)/zw2m)-b1)))
    716616          else
    717617          entr_star(ig,l)=0.
    … …  
    721621             if(l .lt. lalim(ig)) then
    722622
    723 !                detr_star(ig,l)=0.
    724                  detr_star(ig,l) = f_star(ig,l)*zdz*              &
    725             &  adalim
     623                detr_star(ig,l)=0.
    726624
    727625             else
    728626                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    729627            &  ad
    730 !             &  Max(0., 0.0005 - 0.55*zbuoy(ig,l)/zw2m)
    731628
    732629             endif
    … …  
    734631          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    735632                &     MAX(ad,bd*zbuoy(ig,l)/zw2m)
    736 !             &  Max(0.,Min(0.001,0.0005 - 0.55*zbuoy(ig,l)/zw2m))
    737633
    738634          endif
    … …  
    742638          endif
    743639
    744 ! En dessous de lalim, on prend le max de alim_star et entr_star pour
    745 ! alim_star et 0 sinon
     640! If we are still in the source layer, we define the source layer entr. rate  (alim_star) as the
     641! maximum between the source entrainment rate and the estimated entrainment rate.
    746642
    747643        if (l.lt.lalim(ig)) then
    … …  
    750646        endif
    751647
    752 ! Calcul du flux montant normalise
     648! Compute the non-dimensional upward mass flux at layer l+1
     649! using equation 11 of appendix 4.2 in paper
    753650
    754651      f_star(ig,l+1)=f_star(ig,l)+alim_star(ig,l)+entr_star(ig,l)  &
    … …  
    757654      endif
    758655      enddo
    759  
    760       ENDDO   ! of tic
     656! -----------------------------------------------------------------------------------
     657! -----------------------------------------------------------------------------------
     658      ENDDO   ! of internal convergence loop
     659! -----------------------------------------------------------------------------------
     660! -----------------------------------------------------------------------------------
    761661
    762662!---------------------------------------------------------------------------
    763 !initialisations pour le calcul de la hauteur du thermique, de l'inversion et de la vitesse verticale max
     663! Miscellaneous computations for height
    764664!---------------------------------------------------------------------------
    765665
    … …  
    767667            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    768668      IF (lwrite) THEN
    769              print*,'On tombe sur le cas particulier de thermcell_plume'
     669           print*,'thermcell_plume, particular case in velocity profile'
    770670      ENDIF
    771671                zw2(ig,l+1)=0.
    772                 linter(ig)=l+1
    773672            endif
    774673
    775674        if (zw2(ig,l+1).lt.0.) then
    776            linter(ig)=(l*(zw2(ig,l+1)-zw2(ig,l))  &
    777      &               -zw2(ig,l))/(zw2(ig,l+1)-zw2(ig,l))
    778675           zw2(ig,l+1)=0.
    779676        endif
    … …  
    781678
    782679        if (wa_moy(ig,l+1).gt.wmaxa(ig)) then
    783 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
    784680            wmaxa(ig)=wa_moy(ig,l+1)
    785681        endif
    … …  
    787683
    788684!=========================================================================
    789 ! FIN DE LA BOUCLE VERTICALE
    790       enddo
    791685!=========================================================================
    792 
    793 !on recalcule alim_star_tot
     686!=========================================================================
     687! END OF THE LOOP ON VERTICAL LEVELS
     688      enddo
     689!=========================================================================
     690!=========================================================================
     691!=========================================================================
     692
     693! Recompute the source layer total entrainment alim_star_tot
     694! as alim_star may have been modified in the above loop. Renormalization of
     695! alim_star.
     696
    794697       do ig=1,ngridmx
    795698          alim_star_tot(ig)=0.
    … …  
    817720! ===========================================================================
    818721
    819 ! Attention, w2 est transforme en sa racine carree dans cette routine
     722! WARNING, W2 (squared velocity) IS TRANSFORMED IN ITS SQUARE ROOT HERE
    820723
    821724!-------------------------------------------------------------------------------
    822 ! Calcul des caracteristiques du thermique:zmax,wmax
     725! Computations of the thermal height zmax and maximum vertical velocity wmax
    823726!-------------------------------------------------------------------------------
    824727
    825 !calcul de la hauteur max du thermique
     728! Index of the thermal plume height
    826729      do ig=1,ngridmx
    827730         lmax(ig)=lalim(ig)
    … …  
    835738      enddo
    836739
    837 ! On traite le cas particulier qu'il faudrait éviter ou le thermique
    838 ! atteind le haut du modele ...
     740! Particular case when the thermal reached the model top, which is not a good sign
    839741      do ig=1,ngridmx
    840742      if ( zw2(ig,nlayermx) > 1.e-10 ) then
    841           print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
     743          print*,'WARNING !!!!! W2 non-zero in last layer'
    842744          lmax(ig)=nlayermx
    843745      endif
    844746      enddo
    845747
    846 ! pas de thermique si couche 1 stable
    847 !      do ig=1,ngridmx
    848 !         if (lmin(ig).gt.1) then
    849 !             lmax(ig)=1
    850 !             lmin(ig)=1
    851 !             lalim(ig)=1
    852 !         endif
    853 !      enddo
    854 !
    855 ! Determination de zw2 max
     748! Maximum vertical velocity zw2
    856749      do ig=1,ngridmx
    857750         wmax(ig)=0.
    … …  
    872765          enddo
    873766      enddo
    874 !   Longueur caracteristique correspondant a la hauteur des thermiques.
     767
     768! Height of the thermal plume, defined as the following:
     769! zmax=Integral[z*w(z)*dz]/Integral[w(z)*dz]
     770!
    875771      do  ig=1,ngridmx
    876772         zmax(ig)=0.
    … …  
    892788       enddo
    893789
    894 ! Attention, w2 est transforme en sa racine carree dans cette routine
    895 
    896790! ===========================================================================
    897791! ================= FIN HEIGHT ==============================================
    … …  
    904798      endif
    905799
    906 ! Choix de la fonction d'alimentation utilisee pour la fermeture.
     800! alim_star_clos is the source profile used for closure. It consists of the
     801! modified source profile in the source layers, and the entrainment profile
     802! above it.
    907803
    908804      alim_star_clos(:,:)=entr_star(:,:)+alim_star(:,:)
    … …  
    913809
    914810!-------------------------------------------------------------------------------
    915 ! Fermeture,determination de f
     811! Closure, determination of the upward mass flux
    916812!-------------------------------------------------------------------------------
    917 ! Appel avec la version seche
     813! Init.
    918814
    919815      alim_star2(:)=0.
    … …  
    921817      f(:)=0.
    922818
    923 ! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
     819! llmax is the index of the heighest thermal in the simulation domain
    924820      llmax=1
    925821      do ig=1,ngridmx
    … …  
    927823      enddo
    928824
    929 
    930 ! Calcul des integrales sur la verticale de alim_star et de
    931 !   alim_star^2/(rho dz)
     825! Integral of a**2/(rho* Delta z), see equation 13 of appendix 4.2 in paper
     826
    932827      do k=1,llmax-1
    933828         do ig=1,ngridmx
    … …  
    940835      enddo
    941836 
    942 ! WARNING : MARS MODIF : we have added 2. : ratio of wmax/vmoy
    943 ! True ratio is 3.5 but wetake into account the vmoy is the one alimentating
    944 ! the thermal, so there are vs=0 into the vmoy... the true vmoy is lower. (a la louche)
    945 ! And r_aspect has been changed from 2 to 1.5 from observations
     837! Closure mass flux, equation 13 of appendix 4.2 in paper
     838
    946839      do ig=1,ngridmx
    947840         if (alim_star2(ig)>1.e-10) then
    948 !            f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/  &
    949 !      &     (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
    950841             f(ig)=wmax(ig)*alim_star_tot_clos(ig)/  &
    951842      &     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
    … …  
    964855
    965856!-------------------------------------------------------------------------------
    966 !deduction des flux
     857! With the closure mass flux, we can compute the entrainment, detrainment and
     858! upward mass flux from the non-dimensional ones.
    967859!-------------------------------------------------------------------------------
    968860
    969       fomass_max=0.8
    970       alphamax=0.5
    971 
     861      fomass_max=0.8 !maximum mass fraction of a cell that can go upward in an
     862                     ! updraft
     863      alphamax=0.5 !maximum updraft coverage in a cell
     864
     865
     866!    these variables allow to follow corrections made to the mass flux when lwrite=.true.
    972867      ncorecfm1=0
    973868      ncorecfm2=0
    … …  
    981876
    982877!-------------------------------------------------------------------------
    983 ! Multiplication par le flux de masse issu de la femreture
     878! Multiply by the closure mass flux
    984879!-------------------------------------------------------------------------
    985880
    … …  
    988883         detr(:,l)=f(:)*detr_star(:,l)
    989884      enddo
     885
     886! Reconstruct the updraft mass flux everywhere
    990887
    991888      do l=1,zlmax
    … …  
    1002899      enddo
    1003900
    1004 ! Test provisoire : pour comprendre pourquoi on corrige plein de fois
    1005 ! le cas fm6, on commence par regarder une premiere fois avant les
    1006 ! autres corrections.
    1007 
    1008 !      do l=1,nlayermx
    1009 !         do ig=1,ngridmx
    1010 !            if (detr(ig,l).gt.fm(ig,l)) then
    1011 !               ncorecfm8=ncorecfm8+1
    1012 !            endif
    1013 !         enddo
    1014 !      enddo
    1015901
    1016902!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1017 ! FH Version en cours de test;
    1018 ! par rapport a thermcell_flux, on fait une grande boucle sur "l"
    1019 ! et on modifie le flux avec tous les contr�les appliques d'affilee
    1020 ! pour la meme couche
    1021 ! Momentanement, on duplique le calcule du flux pour pouvoir comparer
    1022 ! les flux avant et apres modif
     903!
     904! Now we will reconstruct once again the upward
     905! mass flux, but we will apply corrections
     906! in some cases. We can compare to the
     907! previously computed mass flux (above)
     908!
     909! This verification is done level by level
     910!
    1023911!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1024912
    1025       do l=1,zlmax
     913!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     914
     915      do l=1,zlmax !loop on the levels
     916!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     917
     918! Upward mass flux at level l+1
    1026919
    1027920         do ig=1,ngridmx
    … …  
    1038931
    1039932!-------------------------------------------------------------------------
    1040 ! Verification de la positivite des flux de masse
     933! Upward mass flux should be positive
    1041934!-------------------------------------------------------------------------
    1042935
    … …  
    1061954         enddo
    1062955
    1063 !  Les "optimisations" du flux sont desactivecelles : moins de bidouilles
    1064 !  je considere que celles ci ne sont pas justifiees ou trop delicates
    1065 !  pour MARS, d'apres des observations LES.
    1066956!-------------------------------------------------------------------------
    1067 !Test sur fraca croissant
    1068 !-------------------------------------------------------------------------
    1069 !      if (iflag_thermals_optflux==0) then
    1070 !         do ig=1,ngridmx
    1071 !          if (l.ge.lalim(ig).and.l.le.lmax(ig) &
    1072 !     &    .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
    1073 !!  zzz est le flux en l+1 a frac constant
    1074 !             zzz=fm(ig,l)*rhobarz(ig,l+1)*zw2(ig,l+1)  &
    1075 !     &                          /(rhobarz(ig,l)*zw2(ig,l))
    1076 !             if (fm(ig,l+1).gt.zzz) then
    1077 !                detr(ig,l)=detr(ig,l)+fm(ig,l+1)-zzz
    1078 !                fm(ig,l+1)=zzz
    1079 !                ncorecfm4=ncorecfm4+1
    1080 !             endif
    1081 !          endif
    1082 !        enddo
    1083 !      endif
    1084 !
    1085 !-------------------------------------------------------------------------
    1086 !test sur flux de masse croissant
    1087 !-------------------------------------------------------------------------
    1088 !      if (iflag_thermals_optflux==0) then
    1089 !         do ig=1,ngridmx
    1090 !            if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
    1091 !              f_old=fm(ig,l+1)
    1092 !              fm(ig,l+1)=fm(ig,l)
    1093 !              detr(ig,l)=detr(ig,l)+f_old-fm(ig,l+1)
    1094 !               ncorecfm5=ncorecfm5+1
    1095 !            endif
    1096 !         enddo
    1097 !      endif
    1098 !
    1099 !-------------------------------------------------------------------------
    1100 !detr ne peut pas etre superieur a fm
     957! Detrainment should be lower than upward mass flux
    1101958!-------------------------------------------------------------------------
    1102959
    … …  
    1107964               entr(ig,l)=fm(ig,l+1)
    1108965
    1109 ! Dans le cas ou on est au dessus de la couche d'alimentation et que le
    1110 ! detrainement est plus fort que le flux de masse, on stope le thermique.
    1111 !            endif
     966! When detrainment is stronger than upward mass flux, and we are above the
     967! thermal last level, the plume is stopped
    1112968
    1113969            if(l.gt.lmax(ig)) then
    1114 !            if(l.gt.lalim(ig)) then
    1115970               detr(ig,l)=0.
    1116971               fm(ig,l+1)=0.
    … …  
    1123978
    1124979!-------------------------------------------------------------------------
    1125 !fm ne peut pas etre negatif
     980! Check again for mass flux positivity
    1126981!-------------------------------------------------------------------------
    1127982
    … …  
    1135990
    1136991!-----------------------------------------------------------------------
    1137 !la fraction couverte ne peut pas etre superieure a 1
     992! Fractional coverage should be less than 1
    1138993!-----------------------------------------------------------------------
    1139 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1140 ! FH Partie a revisiter.
    1141 ! Il semble qu'etaient codees ici deux optiques dans le cas
    1142 ! F/ (rho *w) > 1
    1143 ! soit limiter la hauteur du thermique en considerant que c'est
    1144 ! la derniere chouche, soit limiter F a rho w.
    1145 ! Dans le second cas, il faut en fait limiter a un peu moins
    1146 ! que ca parce qu'on a des 1 / ( 1 -alpha) un peu plus loin
    1147 ! dans thermcell_main et qu'il semble de toutes facons deraisonable
    1148 ! d'avoir des fractions de 1..
    1149 ! Ci dessous, et dans l'etat actuel, le premier des  deux if est
    1150 ! sans doute inutile.
    1151 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1152994
    1153995        do ig=1,ngridmx
    … …  
    11641006        enddo
    11651007
    1166 ! Fin de la grande boucle sur les niveaux verticaux
    1167       enddo
     1008      enddo ! on vertical levels
    11681009
    11691010!-----------------------------------------------------------------------
    1170 ! On fait en sorte que la quantite totale d'air entraine dans le
    1171 ! panache ne soit pas trop grande comparee a la masse de la maille
     1011!
     1012! We limit the total mass going from one level to the next, compared to the
     1013! initial total mass fo the cell
     1014!
    11721015!-----------------------------------------------------------------------
    11731016
    … …  
    11821025                entr(ig,l)=entr(ig,l)-eee
    11831026                if ( ddd.gt.0.) then
    1184 !   l'entrainement est trop fort mais l'exces peut etre compense par une
    1185 !   diminution du detrainement)
     1027! The entrainment is too strong but we can compensate the excess by a detrainment decrease
    11861028                   detr(ig,l)=ddd
    11871029                else
    1188 !   l'entrainement est trop fort mais l'exces doit etre compense en partie
    1189 !   par un entrainement plus fort dans la couche superieure
     1030! The entrainment is too strong and we compensate the excess by a stronger entrainment
     1031! in the layer above
    11901032                   if(l.eq.lmax(ig)) then
    11911033                      detr(ig,l)=fm(ig,l)+entr(ig,l)
    … …  
    12001042         enddo
    12011043      enddo
    1202 !
    1203 !              ddd=detr(ig)-entre
    1204 !on s assure que tout s annule bien en zmax
     1044
     1045! Check again that everything cancels at zmax
    12051046      do ig=1,ngridmx
    12061047         fm(ig,lmax(ig)+1)=0.
    … …  
    12101051
    12111052!-----------------------------------------------------------------------
    1212 ! Impression du nombre de bidouilles qui ont ete necessaires
     1053! Summary of the number of modifications that were necessary (if lwrite=.true.
     1054! and only if there were a lot of them)
    12131055!-----------------------------------------------------------------------
    12141056
    … …  
    12371079
    12381080!------------------------------------------------------------------
    1239 !   calcul du transport vertical
     1081! vertical transport computation
    12401082!------------------------------------------------------------------
    12411083
    12421084! ------------------------------------------------------------------
    1243 ! Transport de teta dans l'updraft : (remplace thermcell_dq)
     1085! IN THE UPDRAFT
    12441086! ------------------------------------------------------------------
    12451087
    12461088      zdthladj(:,:)=0.
    1247 
    1248       if (1 .eq. 0) then
    1249 !      call thermcell_dqup(ngridmx,nlayermx,ptimestep     &
    1250 !     &     ,fm,entr,  &
    1251 !     &    masse,ztv,zdthladj)
    1252       else
    1253 
     1089! Based on equation 14 in appendix 4.2
    12541090
    12551091      do ig=1,ngridmx
    … …  
    12701106      enddo
    12711107
    1272       endif
    1273 
    12741108! ------------------------------------------------------------------
    1275 ! Prescription des proprietes du downdraft
     1109! DOWNDRAFT PARAMETERIZATION
    12761110! ------------------------------------------------------------------
    12771111
    … …  
    12811115         if (lmax(ig) .gt. 1) then
    12821116         do l=1,lmax(ig)
    1283 !              if(zlay(ig,l) .le. 0.8*zmax(ig)) then
    12841117              if(zlay(ig,l) .le. zmax(ig)) then
     1118
     1119! see equation 18 of paragraph 48 in paper
    12851120              fm_down(ig,l) =fm(ig,l)* &
    12861121     &      max(fdfu,-4*max(0.,(zlay(ig,l)/zmax(ig)))-0.6)
    12871122              endif
    12881123
    1289 !             if(zlay(ig,l) .le. 0.06*zmax(ig)) then
    1290 !          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(sqrt((zlay(ig,l)/zmax(ig))/0.122449) - 1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
    1291 !             elseif(zlay(ig,l) .le. 0.4*zmax(ig)) then
    1292 !          ztvd(ig,l)=ztv(ig,l)*max(0.,1.-0.25*(ztva(ig,l)/ztv(ig,l) - 1.))
    1293 !             elseif(zlay(ig,l) .le. 0.7*zmax(ig)) then
    1294 !          ztvd(ig,l)=ztv(ig,l)*max(0.,(1.+(((zlay(ig,l)/zmax(ig))-0.7)/1.)*(ztva(ig,l)/ztv(ig,l) - 1.)))
    1295 !             else
    1296 !          ztvd(ig,l)=ztv(ig,l)
    1297 !            endif
    1298 
    1299 !            if(zlay(ig,l) .le. 0.6*zmax(ig)) then
    1300 !          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/179.848 + 0.997832)
    1301 !            elseif(zlay(ig,l) .le. 0.8*zmax(ig)) then
    1302 !          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-171.74)/170.94)
    1303 !             else
    1304 !          ztvd(ig,l)=ztv(ig,l)
    1305 !            endif
    1306 
    1307 
    1308 !            if(zlay(ig,l) .le. 0.8*zmax(ig)) then
    1309 !          ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/224.81 + 0.997832)
    1310 !            elseif(zlay(ig,l) .le. zmax(ig)) then
    1311 !          ztvd(ig,l)=-ztv(ig,l)*(((zlay(ig,l)/zmax(ig))-144.685)/143.885)
    1312 !             else
    1313 !          ztvd(ig,l)=ztv(ig,l)
    1314 !            endif
    1315 
    1316 
    1317 !             if (zbuoy(ig,l) .gt. 0.) then
    1318 !               ztvd(ig,l)=ztva(ig,l)*0.9998
    1319 !!               ztvd(ig,l)=ztv(ig,l)*0.997832
    1320 !!             else
    1321 !!               if(zlay(ig,l) .le. zmax(ig)) then           
    1322 !!               ztvd(ig,l)=ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832)
    1323 !!               endif
    1324 !             endif
    1325 
    13261124            if(zlay(ig,l) .le. zmax(ig)) then           
     1125! see equation 19 of paragraph 49 in paper
    13271126          ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/400. + 0.997832))
    1328 !          ztvd(ig,l)=min(ztv(ig,l),ztv(ig,l)*((zlay(ig,l)/zmax(ig))/299.7 + 0.997832))
    13291127             else
    13301128          ztvd(ig,l)=ztv(ig,l)
    … …  
    13361134
    13371135! ------------------------------------------------------------------
    1338 ! Transport de teta dans le downdraft : (remplace thermcell_dq)
     1136! TRANSPORT IN DOWNDRAFT
    13391137! ------------------------------------------------------------------
    13401138
    … …  
    13441142         if(lmax(ig) .gt. 1) then
    13451143! No downdraft in the very-near surface layer, we begin at k=3
     1144! Based on equation 14 in appendix 4.2
    13461145 
    13471146         do k=3,lmax(ig)
    … …  
    13611160
    13621161!------------------------------------------------------------------
    1363 ! Calcul de la fraction de l'ascendance
     1162! Final fraction coverage with the clean upward mass flux, computed at interfaces
    13641163!------------------------------------------------------------------
    13651164      fraca(:,:)=0.
    … …  
    13741173      enddo
    13751174
    1376 
    1377 
    1378 ! ===========================================================================
    1379 ! ============= DV2 =========================================================
    1380 ! ===========================================================================
    1381 ! ROUTINE OVERIDE : ne prends pas en compte le downdraft
    1382 ! de plus, le gradient de pression horizontal semble tout deregler... A VOIR
    1383 
    1384       if (0 .eq. 1) then
    1385 
    13861175!------------------------------------------------------------------
    1387 !  calcul du transport vertical du moment horizontal
    1388 !------------------------------------------------------------------
    1389 
    1390 ! Calcul du transport de V tenant compte d'echange par gradient
    1391 ! de pression horizontal avec l'environnement
    1392 
    1393 !   calcul du detrainement
    1394 !---------------------------
    1395 
    1396       nlarga0=0.
    1397 
    1398       do k=1,nlayermx
    1399          do ig=1,ngridmx
    1400             detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    1401          enddo
    1402       enddo
    1403 
    1404 !   calcul de la valeur dans les ascendances
    1405       do ig=1,ngridmx
    1406          zua(ig,1)=zu(ig,1)
    1407          zva(ig,1)=zv(ig,1)
    1408          ue(ig,1)=zu(ig,1)
    1409          ve(ig,1)=zv(ig,1)
    1410       enddo
    1411 
    1412       gamma(1:ngridmx,1)=0.
    1413       do k=2,nlayermx
    1414          do ig=1,ngridmx
    1415             ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))*ptimestep > 1.e-5*masse(ig,k)
    1416             if(ltherm(ig,k).and.zmax(ig)>0.) then
    1417                gamma0(ig,k)=masse(ig,k)  &
    1418      &         *sqrt( 0.5*(fraca(ig,k+1)+fraca(ig,k)) )  &
    1419      &         *0.5/zmax(ig)  &
    1420      &         *1.
    1421             else
    1422                gamma0(ig,k)=0.
    1423             endif
    1424             if (ltherm(ig,k).and.zmax(ig)<=0.) nlarga0=nlarga0+1
    1425           enddo
    1426       enddo
    1427 
    1428       gamma(:,:)=0.
    1429 
    1430       do k=2,nlayermx
    1431 
    1432          do ig=1,ngridmx
    1433 
    1434             if (ltherm(ig,k)) then
    1435                dua(ig,k)=zua(ig,k-1)-zu(ig,k-1)
    1436                dva(ig,k)=zva(ig,k-1)-zv(ig,k-1)
    1437             else
    1438                zua(ig,k)=zu(ig,k)
    1439                zva(ig,k)=zv(ig,k)
    1440                ue(ig,k)=zu(ig,k)
    1441                ve(ig,k)=zv(ig,k)
    1442             endif
    1443          enddo
    1444 
    1445 
    1446 ! Debut des iterations
    1447 !----------------------
    1448 
    1449 ! AC WARNING : SALE !
    1450 
    1451       do iter=1,5
    1452          do ig=1,ngridmx
    1453 ! Pour memoire : calcul prenant en compte la fraction reelle
    1454 !              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
    1455 !              zf2=1./(1.-zf)
    1456 ! Calcul avec fraction infiniement petite
    1457                zf=0.
    1458                zf2=1.
    1459 
    1460 !  la première fois on multiplie le coefficient de freinage
    1461 !  par le module du vent dans la couche en dessous.
    1462 !  Mais pourquoi donc ???
    1463                if (ltherm(ig,k)) then
    1464 !   On choisit une relaxation lineaire.
    1465 !                 gamma(ig,k)=gamma0(ig,k)
    1466 !   On choisit une relaxation quadratique.
    1467                 gamma(ig,k)=gamma0(ig,k)*sqrt(dua(ig,k)**2+dva(ig,k)**2)
    1468                   zua(ig,k)=(fm(ig,k)*zua(ig,k-1)  &
    1469      &               +(zf2*entr(ig,k)+gamma(ig,k))*zu(ig,k))  &
    1470      &               /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2  &
    1471      &                 +gamma(ig,k))
    1472                   zva(ig,k)=(fm(ig,k)*zva(ig,k-1)  &
    1473      &               +(zf2*entr(ig,k)+gamma(ig,k))*zv(ig,k))  &
    1474      &               /(fm(ig,k+1)+detr_dv2(ig,k)+entr(ig,k)*zf*zf2  &
    1475      &                 +gamma(ig,k))
    1476 
    1477 !                  print*,' OUTPUT DV2 !!!!!!!!!!!!!!!!!!!!',k,zua(ig,k),zva(ig,k),zu(ig,k),zv(ig,k),dua(ig,k),dva(ig,k)
    1478                   dua(ig,k)=zua(ig,k)-zu(ig,k)
    1479                   dva(ig,k)=zva(ig,k)-zv(ig,k)
    1480                   ue(ig,k)=(zu(ig,k)-zf*zua(ig,k))*zf2
    1481                   ve(ig,k)=(zv(ig,k)-zf*zva(ig,k))*zf2
    1482                endif
    1483          enddo
    1484 ! Fin des iterations
    1485 !--------------------
    1486       enddo
    1487 
    1488       enddo ! k=2,nlayermx
    1489 
    1490 ! Calcul du flux vertical de moment dans l'environnement.
    1491 !---------------------------------------------------------
    1492       do k=2,nlayermx
    1493          do ig=1,ngridmx
    1494             wud(ig,k)=fm(ig,k)*ue(ig,k)
    1495             wvd(ig,k)=fm(ig,k)*ve(ig,k)
    1496          enddo
    1497       enddo
    1498       do ig=1,ngridmx
    1499          wud(ig,1)=0.
    1500          wud(ig,nlayermx+1)=0.
    1501          wvd(ig,1)=0.
    1502          wvd(ig,nlayermx+1)=0.
    1503       enddo
    1504 
    1505 ! calcul des tendances.
    1506 !-----------------------
    1507       do k=1,nlayermx
    1508          do ig=1,ngridmx
    1509             pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k)  &
    1510      &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
    1511      &               -wud(ig,k)+wud(ig,k+1))  &
    1512      &               /masse(ig,k)
    1513             pdvadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zva(ig,k)  &
    1514      &               -(entr(ig,k)+gamma(ig,k))*ve(ig,k)  &
    1515      &               -wvd(ig,k)+wvd(ig,k+1))  &
    1516      &               /masse(ig,k)
    1517          enddo
    1518       enddo
    1519 
    1520 
    1521 ! Sorties eventuelles.
    1522 !----------------------
    1523 
    1524 !      if (nlarga0>0) then
    1525 !          print*,'WARNING !!!!!! DANS THERMCELL_DV2 '
    1526 !      print*,nlarga0,' points pour lesquels laraga=0. dans un thermique'
    1527 !          print*,'Il faudrait decortiquer ces points'
    1528 !      endif
    1529 
    1530 ! ===========================================================================
    1531 ! ============= FIN DV2 =====================================================
    1532 ! ===========================================================================
    1533 
    1534       else
    1535 !      detrmod(:,:)=0.
    1536 !      do k=1,zlmax
    1537 !         do ig=1,ngridmx
    1538 !            detrmod(ig,k)=fm(ig,k)-fm(ig,k+1) &
    1539 !     &      +entr(ig,k)
    1540 !            if (detrmod(ig,k).lt.0.) then
    1541 !               entr(ig,k)=entr(ig,k)-detrmod(ig,k)
    1542 !               detrmod(ig,k)=0.
    1543 !            endif
    1544 !         enddo
    1545 !      enddo
    1546 !
    1547 !
    1548 !      call thermcell_dqup(ngridmx,nlayermx,ptimestep                &
    1549 !     &      ,fm,entr,detrmod,  &
    1550 !     &     masse,zu,pduadj,ndt,zlmax)
    1551 !
    1552 !      call thermcell_dqup(ngridmx,nlayermx,ptimestep    &
    1553 !     &       ,fm,entr,detrmod,  &
    1554 !     &     masse,zv,pdvadj,ndt,zlmax)
    1555 
    1556       endif
    1557 
    1558 !------------------------------------------------------------------
    1559 !  calcul du transport vertical de traceurs
     1176! Transport of C02 Tracer
    15601177!------------------------------------------------------------------
    15611178
    … …  
    16001217         do ig=1,ngridmx
    16011218         pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)*ratiom(ig,l)
    1602 !         pdtadj(ig,l)=zdthladj(ig,l)*zpopsk(ig,l)*ratiom(ig,l)
    1603          enddo
    1604       enddo
     1219         enddo
     1220      enddo
     1221
     1222! ===========================================================================
     1223! ============= FIN TRANSPORT ===============================================
     1224! ===========================================================================
     1225
    16051226
    16061227!------------------------------------------------------------------
    1607 !  calcul du transport vertical de la tke
     1228!   Diagnostics for outputs
    16081229!------------------------------------------------------------------
    1609 
    1610 !      modname='tke'
    1611 !      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
    1612 !      &      masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)
    1613 
    1614 ! ===========================================================================
    1615 ! ============= FIN TRANSPORT ===============================================
    1616 ! ===========================================================================
    1617 
    1618 
    1619 !------------------------------------------------------------------
    1620 !   Calculs de diagnostiques pour les sorties
    1621 !------------------------------------------------------------------
    1622 ! DIAGNOSTIQUE
    16231230! We compute interface values for teta env and th. The last interface
    16241231! value does not matter, as the mass flux is 0 there.
Note: See TracChangeset for help on using the changeset viewer.