Changeset 185


Ignore:
Timestamp:
Jul 1, 2011, 5:22:34 PM (14 years ago)
Author:
acolaitis
Message:

17/06/2011 == AC

  • Added new settings for the Martian thermals from new LES observations
  • Revamped thermcell's module variables to allow it's removal
  • Minor changes in physiq and meso_physiq for the call to thermals
  • Switched from dynamic to static memory allocation for all thermals variable

to gain computation speed

  • Thermals now output maximum speed, maximum heat flux, and maximum height in thermal plumes
Location:
trunk/LMDZ.MARS
Files:
1 deleted
7 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/README

    r174 r185  
    763763- correct bug (introduced previously) in lect_start_archive.F on loop
    764764  boundaries for soil temperature.
     765
     766== 17/06/2011 == AC
     767- Added new settings for the Martian thermals from new LES observations
     768- Revamped thermcell's module variables to allow it's removal
     769- Minor changes in physiq and meso_physiq for the call to thermals
     770- Switched from dynamic to static memory allocation for all thermals variable
     771to gain computation speed

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

    r173 r185  
    22! AC 2011-01-05
    33!
    4       SUBROUTINE calltherm_interface (ngrid,nlayer,firstcall, &
     4      SUBROUTINE calltherm_interface (firstcall, &
    55     & long,lati,zzlev,zzlay, &
    66     & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, &
    7      & pplay,pplev,pphi,nq,zpopsk, &
    8      & pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,pbl_dtke,hfmax,wmax)
     7     & pplay,pplev,pphi,zpopsk, &
     8     & pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,pbl_dtke,hfmax,wmax)
    99
    1010       USE ioipsl_getincom
    … …  
    1212      implicit none
    1313#include "callkeys.h"
     14#include "dimensions.h"
     15#include "dimphys.h"
     16
    1417!--------------------------------------------------------
    1518! Variables d'entree
    1619!--------------------------------------------------------
    1720
    18       INTEGER, INTENT(IN) :: ngrid,nlayer,nq
    1921      REAL, INTENT(IN) :: ptimestep
    20       REAL, INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
    21       REAL, INTENT(IN) :: pphi(ngrid,nlayer)
    22       REAL, INTENT(IN) :: pu(ngrid,nlayer),pv(ngrid,nlayer)
    23       REAL, INTENT(IN) :: pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
    24       REAL, INTENT(IN) :: zzlay(ngrid,nlayer)
    25       REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1)
     22      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1),pplay(ngridmx,nlayermx)
     23      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
     24      REAL, INTENT(IN) :: pu(ngridmx,nlayermx),pv(ngridmx,nlayermx)
     25      REAL, INTENT(IN) :: pt(ngridmx,nlayermx),pq(ngridmx,nlayermx,nqmx)
     26      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
     27      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
    2628      LOGICAL, INTENT(IN) :: firstcall
    27       REAL, INTENT(IN) :: pdu(ngrid,nlayer),pdv(ngrid,nlayer)
    28       REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq),pdt(ngrid,nlayer)
    29       REAL, INTENT(IN) :: q2(ngrid,nlayer+1)
    30       REAL, INTENT(IN) :: long(ngrid),lati(ngrid)
    31       REAL, INTENT(IN) :: zpopsk(ngrid,nlayer)
     29      REAL, INTENT(IN) :: pdu(ngridmx,nlayermx),pdv(ngridmx,nlayermx)
     30      REAL, INTENT(IN) :: pdq(ngridmx,nlayermx,nqmx),pdt(ngridmx,nlayermx)
     31      REAL, INTENT(IN) :: q2(ngridmx,nlayermx+1)
     32      REAL, INTENT(IN) :: long(ngridmx),lati(ngridmx)
     33      REAL, INTENT(IN) :: zpopsk(ngridmx,nlayermx)
    3234
    3335!--------------------------------------------------------
    … …  
    3537!--------------------------------------------------------
    3638
    37       REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
    38       REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
    39       INTEGER lmax_th(ngrid)
    40       REAL pbl_dtke(ngrid,nlayer+1)
     39      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
     40      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
     41      INTEGER lmax_th(ngridmx)
     42      REAL zmax_th(ngridmx)
     43      REAL pbl_dtke(ngridmx,nlayermx+1)
    4144
    4245!--------------------------------------------------------
    4346! Variables du thermique
    4447!--------------------------------------------------------
    45       REAL u_seri(ngrid,nlayer), v_seri(ngrid,nlayer)
    46       REAL t_seri(ngrid,nlayer)
    47       REAL d_t_ajs(ngrid,nlayer)
    48       REAL d_u_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq)
    49       REAL d_v_ajs(ngrid,nlayer)
    50       REAL fm_therm(ngrid,nlayer+1), entr_therm(ngrid,nlayer)
    51       REAL detr_therm(ngrid,nlayer)
    52       REAL zw2(ngrid,nlayer+1)
    53       REAL fraca(ngrid,nlayer+1)
    54       REAL ztla(ngrid,nlayer)
    55       REAL q_therm(ngrid,nlayer), pq_therm(ngrid,nlayer,nq)
    56       REAL dq_therm(ngrid,nlayer), dq_thermdown(ngrid,nlayer)
    57       REAL q2_therm(ngrid,nlayer), dq2_therm(ngrid,nlayer)
     48      REAL u_seri(ngridmx,nlayermx), v_seri(ngridmx,nlayermx)
     49      REAL t_seri(ngridmx,nlayermx)
     50      REAL d_t_ajs(ngridmx,nlayermx)
     51      REAL d_u_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
     52      REAL d_v_ajs(ngridmx,nlayermx)
     53      REAL fm_therm(ngridmx,nlayermx+1), entr_therm(ngridmx,nlayermx)
     54      REAL detr_therm(ngridmx,nlayermx)
     55      REAL zw2(ngridmx,nlayermx+1)
     56      REAL fraca(ngridmx,nlayermx+1)
     57      REAL ztla(ngridmx,nlayermx)
     58      REAL q_therm(ngridmx,nlayermx), pq_therm(ngridmx,nlayermx,nqmx)
     59      REAL dq_therm(ngridmx,nlayermx), dq_thermdown(ngridmx,nlayermx)
     60      REAL q2_therm(ngridmx,nlayermx), dq2_therm(ngridmx,nlayermx)
    5861
    5962      LOGICAL qtransport_thermals,dtke_thermals
    … …  
    6366! Variable de diagnostique : flux de chaleur vertical
    6467
    65       REAL heatFlux(ngrid,nlayer)
    66       REAL heatFlux_down(ngrid,nlayer)
    67       REAL buoyancyOut(ngrid,nlayer)
    68       REAL buoyancyEst(ngrid,nlayer)
    69       REAL hfmax(ngrid),wmax(ngrid)
     68      REAL heatFlux(ngridmx,nlayermx)
     69      REAL heatFlux_down(ngridmx,nlayermx)
     70      REAL buoyancyOut(ngridmx,nlayermx)
     71      REAL buoyancyEst(ngridmx,nlayermx)
     72      REAL hfmax(ngridmx),wmax(ngridmx)
     73
     74      REAL tstart,tstop
    7075
    7176!---------------------------------------------------------
    … …  
    125130         if(dtke_thermals) then
    126131
    127          DO l=1,nlayer
     132         DO l=1,nlayermx
    128133              q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1))
    129134         ENDDO
    130135         endif
    131136
    132          CALL calltherm_mars(ngrid,nlayer,ptimestep,nq,zzlev,zzlay &
     137         call cpu_time(tstart)
     138
     139         CALL calltherm_mars(ptimestep,zzlev,zzlay &
    133140     &      ,pplay,pplev,pphi &
    134141     &      ,u_seri,v_seri,t_seri,pq_therm, q2_therm &
    135142     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, dq2_therm &
    136143     &      ,fm_therm,entr_therm,detr_therm &
    137      &      ,lmax_th &
     144     &      ,lmax_th,zmax_th &
    138145     &      ,zw2,fraca &
    139146     &      ,zpopsk,ztla,heatFlux,heatFlux_down &
    140147     &      ,buoyancyOut,buoyancyEst,hfmax,wmax)
     148          call cpu_time(tstop)
     149         print*,'TOTAL elapsed time in thermals : ',tstop-tstart
    141150
    142151
    … …  
    156165
    157166
    158          DO l=2,nlayer
     167         DO l=2,nlayermx
    159168              pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l))/ptimestep
    160169         ENDDO
    161170
    162171         pbl_dtke(:,1)=0.5*dq2_therm(:,1)/ptimestep
    163          pbl_dtke(:,nlayer+1)=0.
     172         pbl_dtke(:,nlayermx+1)=0.
    164173!! DIAGNOSTICS
    165174       
    166175        if(outptherm) then
    167         if (ngrid .eq. 1) then
    168         call WRITEDIAGFI(ngrid,'entr_therm','entrainement thermique',&
     176        if (ngridmx .eq. 1) then
     177        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
    169178     &                       'kg/m-2',1,entr_therm)
    170         call WRITEDIAGFI(ngrid,'detr_therm','detrainement thermique',&
     179        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
    171180     &                       'kg/m-2',1,detr_therm)
    172         call WRITEDIAGFI(ngrid,'fm_therm','flux masse thermique',&
     181        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
    173182     &                       'kg/m-2',1,fm_therm)
    174         call WRITEDIAGFI(ngrid,'zw2','vitesse verticale thermique',&
     183        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
    175184     &                       'm/s',1,zw2)
    176         call WRITEDIAGFI(ngrid,'heatFlux_up','heatFlux_updraft',&
     185        call WRITEDIAGFI(ngridmx,'heatFlux_up','heatFlux_updraft',&
    177186     &                       'SI',1,heatFlux)
    178        call WRITEDIAGFI(ngrid,'heatFlux_down','heatFlux_downdraft',&
     187       call WRITEDIAGFI(ngridmx,'heatFlux_down','heatFlux_downdraft',&
    179188     &                       'SI',1,heatFlux_down)
    180         call WRITEDIAGFI(ngrid,'fraca','fraction coverage',&
     189        call WRITEDIAGFI(ngridmx,'fraca','fraction coverage',&
    181190     &                       'percent',1,fraca)
    182         call WRITEDIAGFI(ngrid,'buoyancyOut','buoyancyOut',&
     191        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
    183192     &                       'm.s-2',1,buoyancyOut)
    184         call WRITEDIAGFI(ngrid,'buoyancyEst','buoyancyEst',&
     193        call WRITEDIAGFI(ngridmx,'buoyancyEst','buoyancyEst',&
    185194     &                       'm.s-2',1,buoyancyEst)
    186         call WRITEDIAGFI(ngrid,'d_t_th',  &
     195        call WRITEDIAGFI(ngridmx,'d_t_th',  &
    187196     &         'tendance temp TH','K',1,d_t_ajs)
     197        call WRITEDIAGFI(ngridmx,'zmax',  &
     198     &         'pbl height','m',0,zmax_th)
    188199      else
    189200
    190         call WRITEDIAGFI(ngrid,'entr_therm','entrainement thermique',&
     201        call WRITEDIAGFI(ngridmx,'entr_therm','entrainement thermique',&
    191202     &                       'kg/m-2',3,entr_therm)
    192         call WRITEDIAGFI(ngrid,'detr_therm','detrainement thermique',&
     203        call WRITEDIAGFI(ngridmx,'detr_therm','detrainement thermique',&
    193204     &                       'kg/m-2',3,detr_therm)
    194         call WRITEDIAGFI(ngrid,'fm_therm','flux masse thermique',&
     205        call WRITEDIAGFI(ngridmx,'fm_therm','flux masse thermique',&
    195206     &                       'kg/m-2',3,fm_therm)
    196         call WRITEDIAGFI(ngrid,'zw2','vitesse verticale thermique',&
     207        call WRITEDIAGFI(ngridmx,'zw2','vitesse verticale thermique',&
    197208     &                       'm/s',3,zw2)
    198         call WRITEDIAGFI(ngrid,'heatFlux','heatFlux',&
     209        call WRITEDIAGFI(ngridmx,'heatFlux','heatFlux',&
    199210     &                       'SI',3,heatFlux)
    200         call WRITEDIAGFI(ngrid,'buoyancyOut','buoyancyOut',&
     211        call WRITEDIAGFI(ngridmx,'buoyancyOut','buoyancyOut',&
    201212     &                       'SI',3,buoyancyOut)
    202         call WRITEDIAGFI(ngrid,'d_t_th',  &
     213        call WRITEDIAGFI(ngridmx,'d_t_th',  &
    203214     &         'tendance temp TH','K',3,d_t_ajs)
    204215

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

    r173 r185  
    22! $Id: calltherm.F90 1428 2010-09-13 08:43:37Z fairhead $
    33!
    4       subroutine calltherm_mars(ngrid,nlayer,dtime,nq,zzlev,zzlay  &
     4      subroutine calltherm_mars(dtime,zzlev,zzlay  &
    55     &      ,pplay,paprs,pphi  &
    66     &      ,u_seri,v_seri,t_seri,pq_therm,q2_therm  &
    77     &      ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs,dq2_therm  &
    8      &      ,fm_therm,entr_therm,detr_therm,lmax,&
     8     &      ,fm_therm,entr_therm,detr_therm,lmax,zmax,&
    99     &   zw2,fraca,zpopsk,ztla,heatFlux,heatFlux_down,&
    1010     &     buoyancyOut,buoyancyEst,hfmax,wmax)
    1111
    12        USE thermcell, only : nsplit_thermals,r_aspect_thermals
    1312       USE ioipsl_getincom
    1413      implicit none
    1514
    16       INTEGER, INTENT(IN) :: ngrid,nlayer
     15#include "dimensions.h"
     16#include "dimphys.h"
     17
    1718      REAL dtime
    18       LOGICAL logexpr0, logexpr2(ngrid,nlayer), logexpr1(ngrid)
     19      LOGICAL logexpr0, logexpr2(ngridmx,nlayermx), logexpr1(ngridmx)
    1920      REAL fact
    20       INTEGER nbptspb,nq
    21 
    22       REAL, INTENT(IN) :: zzlay(ngrid,nlayer)
    23       REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1)
    24 
    25       REAL u_seri(ngrid,nlayer),v_seri(ngrid,nlayer)
    26       REAL t_seri(ngrid,nlayer),pq_therm(ngrid,nlayer,nq)
    27       REAL q2_therm(ngrid,nlayer)
    28       REAL paprs(ngrid,nlayer+1)
    29       REAL pplay(ngrid,nlayer)
    30       REAL pphi(ngrid,nlayer)
    31       real zlev(ngrid,nlayer+1)
     21      INTEGER nbptspb
     22
     23      REAL, INTENT(IN) :: zzlay(ngridmx,nlayermx)
     24      REAL, INTENT(IN) :: zzlev(ngridmx,nlayermx+1)
     25
     26      REAL u_seri(ngridmx,nlayermx),v_seri(ngridmx,nlayermx)
     27      REAL t_seri(ngridmx,nlayermx),pq_therm(ngridmx,nlayermx,nqmx)
     28      REAL q2_therm(ngridmx,nlayermx)
     29      REAL paprs(ngridmx,nlayermx+1)
     30      REAL pplay(ngridmx,nlayermx)
     31      REAL pphi(ngridmx,nlayermx)
     32      real zlev(ngridmx,nlayermx+1)
    3233!test: on sort lentr et a* pour alimenter KE
    33       REAL zw2(ngrid,nlayer+1),fraca(ngrid,nlayer+1)
    34       REAL zzw2(ngrid,nlayer+1)
     34      REAL zw2(ngridmx,nlayermx+1),fraca(ngridmx,nlayermx+1)
     35      REAL zzw2(ngridmx,nlayermx+1)
    3536
    3637!FH Update Thermiques
    37       REAL d_t_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq)
    38       REAL d_u_ajs(ngrid,nlayer),d_v_ajs(ngrid,nlayer)
    39       REAL dq2_therm(ngrid,nlayer), dq2_the(ngrid,nlayer)
    40       real fm_therm(ngrid,nlayer+1)
    41       real entr_therm(ngrid,nlayer),detr_therm(ngrid,nlayer)
     38      REAL d_t_ajs(ngridmx,nlayermx), d_q_ajs(ngridmx,nlayermx,nqmx)
     39      REAL d_u_ajs(ngridmx,nlayermx),d_v_ajs(ngridmx,nlayermx)
     40      REAL dq2_therm(ngridmx,nlayermx), dq2_the(ngridmx,nlayermx)
     41      real fm_therm(ngridmx,nlayermx+1)
     42      real entr_therm(ngridmx,nlayermx),detr_therm(ngridmx,nlayermx)
    4243
    4344!********************************************************
    4445!     declarations
    45       real zpopsk(ngrid,nlayer)
    46       real ztla(ngrid,nlayer)
    47       real wmax(ngrid)
    48       real hfmax(ngrid)
    49       integer lmax(ngrid)
     46      real zpopsk(ngridmx,nlayermx)
     47      real ztla(ngridmx,nlayermx)
     48      real wmax(ngridmx)
     49      real hfmax(ngridmx)
     50      integer lmax(ngridmx)
     51      real zmax(ngridmx)
    5052
    5153!nouvelles variables pour la convection
    … …  
    5658
    5759! variables locales
    58       REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq)
    59       REAL d_u_the(ngrid,nlayer),d_v_the(ngrid,nlayer)
     60      REAL d_t_the(ngridmx,nlayermx), d_q_the(ngridmx,nlayermx,nqmx)
     61      REAL d_u_the(ngridmx,nlayermx),d_v_the(ngridmx,nlayermx)
    6062!
    61       integer isplit
    62       real zfm_therm(ngrid,nlayer+1),zdt
    63       real zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer)
    64       real heatFlux(ngrid,nlayer)
    65       real heatFlux_down(ngrid,nlayer)
    66       real buoyancyOut(ngrid,nlayer)
    67       real buoyancyEst(ngrid,nlayer)
    68       real zheatFlux(ngrid,nlayer)
    69       real zheatFlux_down(ngrid,nlayer)
    70       real zbuoyancyOut(ngrid,nlayer)
    71       real zbuoyancyEst(ngrid,nlayer)
     63      integer isplit,nsplit_thermals
     64      real r_aspect_thermals
     65
     66      real zfm_therm(ngridmx,nlayermx+1),zdt
     67      real zentr_therm(ngridmx,nlayermx),zdetr_therm(ngridmx,nlayermx)
     68      real heatFlux(ngridmx,nlayermx)
     69      real heatFlux_down(ngridmx,nlayermx)
     70      real buoyancyOut(ngridmx,nlayermx)
     71      real buoyancyEst(ngridmx,nlayermx)
     72      real zheatFlux(ngridmx,nlayermx)
     73      real zheatFlux_down(ngridmx,nlayermx)
     74      real zbuoyancyOut(ngridmx,nlayermx)
     75      real zbuoyancyEst(ngridmx,nlayermx)
    7276
    7377      character (len=20) :: modname='calltherm'
    … …  
    7781      logical, save :: first=.true.
    7882
     83      REAL tstart,tstop
     84
     85
    7986!  Modele du thermique
    8087!  ===================
    81 
    82          nsplit_thermals=20
     88 
     89         r_aspect_thermals=3.
     90         nsplit_thermals=40
    8391         call getin("nsplit_thermals",nsplit_thermals)
    8492
    … …  
    98106         do isplit=1,nsplit_thermals
    99107
     108!         call cpu_time(tstart)
     109
     110
    100111! On reinitialise les flux de masse a zero pour le cumul en
    101112! cas de splitting
    … …  
    116127         d_v_the(:,:)=0.
    117128         dq2_the(:,:)=0.
    118          if (nq .ne. 0) then
     129         if (nqmx .ne. 0) then
    119130            d_q_the(:,:,:)=0.
    120131         endif
    121132
    122              CALL thermcell_main_mars(ngrid,nlayer,nq,zdt  &
     133             CALL thermcell_main_mars(zdt  &
    123134     &      ,pplay,paprs,pphi,zzlev,zzlay  &
    124135     &      ,u_seri,v_seri,t_seri,pq_therm,q2_therm  &
    125136     &      ,d_u_the,d_v_the,d_t_the,d_q_the,dq2_the  &
    126      &      ,zfm_therm,zentr_therm,zdetr_therm,lmax  &
     137     &      ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax  &
    127138     &      ,r_aspect_thermals &
    128139     &      ,zzw2,fraca,zpopsk &
    … …  
    138149            dq2_the(:,:)=dq2_the(:,:)*fact           
    139150
    140             if (nq .ne. 0) then
     151            if (nqmx .ne. 0) then
    141152               d_q_the(:,:,:)=d_q_the(:,:,:)*fact
    142153            endif
    … …  
    175186            q2_therm(:,:) = q2_therm(:,:) + dq2_therm(:,:)
    176187
     188
     189!         call cpu_time(tstop)
     190!         print*,'elapsed time in thermals : ',tstop-tstart
     191
    177192         enddo ! isplit
    178193
    … …  
    180195!****************************************************************
    181196
    182 !          do i=1,ngrid
    183 !             do k=1,nlayer
     197!          do i=1,ngridmx
     198!             do k=1,nlayermx
    184199!                if (ztla(i,k) .lt. 1.e-10) fraca(i,k) =0.
    185200!               print*,'youpi je sers a quelque chose !'
    … …  
    187202!          enddo
    188203       
    189           DO i=1,ngrid
     204          DO i=1,ngridmx
    190205            hfmax(i)=MAXVAL(heatFlux(i,:)+heatFlux_down(i,:))
    191206            wmax(i)=MAXVAL(zw2(i,:))

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

    r173 r185  
    361361c Variables from thermal
    362362
    363       REAL lmax_th_out(ngrid)
    364       REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
    365       REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
    366       INTEGER lmax_th(ngrid)
    367       REAL dtke_th(ngrid,nlayer+1)
    368       REAL wmax_th(ngrid),hfmax_th(ngrid)
     363      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
     364      REAL wmax_th(ngridmx)
     365      REAL ,SAVE :: hfmax_th(ngridmx)
     366      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
     367      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
     368      INTEGER lmax_th(ngridmx)
     369      REAL dtke_th(ngridmx,nlayermx+1)
     370      REAL dummycol(ngridmx)
     371
    369372c=======================================================================
    370373#ifdef MESOSCALE
    … …  
    403406      PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngridmx,nsoilmx)
    404407      mlayer(0:nsoilmx-1)=wdsoil(1,:)
    405       PRINT*,'check: layer ', mlayer
     408      PRINT*,'check: midlayer ', mlayer(:)
    406409            !!!!!!!!!!!!!!!!! DONE in soil_setting.F
    407410            ! 1.5 Build layer(); following the same law as mlayer()
    … …  
    413416              layer(iloop)=lay1*(alpha**(iloop-1))
    414417            enddo
     418
     419      PRINT*,'check: layer ', layer(:)
     420
    415421            !!!!!!!!!!!!!!!!! DONE in soil_setting.F
    416422      tnom(:)=wtnom(:)   !! est rempli dans advtrac.h
    … …  
    10091015      if(calltherm) then
    10101016 
    1011         call calltherm_interface(ngrid,nlayer,firstcall,
     1017        call calltherm_interface(firstcall,
    10121018     $ long,lati,zzlev,zzlay,
    10131019     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
    1014      $ pplay,pplev,pphi,nq,zpopsk,
    1015      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th,hfmax_th,wmax_th)
    1016  
     1020     $ pplay,pplev,pphi,zpopsk,
     1021     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
     1022     $ dtke_th,hfmax_th,wmax_th)
     1023
    10171024         DO l=1,nlayer
    10181025           DO ig=1,ngrid

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

    r173 r185  
    206206
    207207      CHARACTER*80 fichier
    208       INTEGER l,ig,ierr,igout,iq,i, tapphys
     208      INTEGER l,ig,ierr,igout,iq,i, tapphys,k
     209      LOGICAL end_of_file
    209210
    210211      REAL fluxsurf_lw(ngridmx)      !incident LW (IR) surface flux (W.m-2)
    … …  
    302303c Variables from thermal
    303304
    304       REAL lmax_th_out(ngrid)
    305       REAL wmax_th(ngrid),hfmax_th(ngrid)
    306       REAL pdu_th(ngrid,nlayer),pdv_th(ngrid,nlayer)
    307       REAL pdt_th(ngrid,nlayer),pdq_th(ngrid,nlayer,nq)
    308       INTEGER lmax_th(ngrid)
    309       REAL dtke_th(ngrid,nlayer+1)
    310       REAL dummycol(ngrid)
    311       REAL hfx(ngrid)
     305      REAL lmax_th_out(ngridmx),zmax_th(ngridmx)
     306      REAL wmax_th(ngridmx)
     307      REAL ,SAVE :: hfmax_th(ngridmx)
     308      REAL pdu_th(ngridmx,nlayermx),pdv_th(ngridmx,nlayermx)
     309      REAL pdt_th(ngridmx,nlayermx),pdq_th(ngridmx,nlayermx,nqmx)
     310      INTEGER lmax_th(ngridmx)
     311      REAL dtke_th(ngridmx,nlayermx+1)
     312      REAL dummycol(ngridmx)
    312313
    313314c=======================================================================
    … …  
    562563           ENDIF
    563564
    564 
    565 
    566565        ENDIF ! of if(mod(icount-1,iradia).eq.0)
    567566
    … …  
    578577               fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig)
    579578           ENDDO
    580 
    581579
    582580         DO l=1,nlayer
    … …  
    624622            enddo
    625623         enddo
    626          
     624
    627625c        Calling vdif (Martian version WITH CO2 condensation)
    628626         CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk,
    … …  
    634632     &        zdqdif,zdqsdif)
    635633
    636           hfx(:) = zflubid(:)-capcal(:)*zdtsdif(:)
    637 
    638634         DO l=1,nlayer
    639635            DO ig=1,ngrid
    … …  
    647643         ENDDO
    648644
    649          DO ig=1,ngrid
    650             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
    651          ENDDO
     645          DO ig=1,ngrid
     646             zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)
     647          ENDDO
    652648
    653649         if (tracer) then
    … …  
    680676      if(calltherm) then
    681677 
    682         call calltherm_interface(ngrid,nlayer,firstcall,
     678        call calltherm_interface(firstcall,
    683679     $ long,lati,zzlev,zzlay,
    684680     $ ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2,
    685      $ pplay,pplev,pphi,nq,zpopsk,
    686      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,dtke_th,hfmax_th,wmax_th)
     681     $ pplay,pplev,pphi,zpopsk,
     682     $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
     683     $ dtke_th,hfmax_th,wmax_th)
    687684 
    688685         DO l=1,nlayer
    … …  
    767764     $              co2ice,albedo,emis,
    768765     $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
    769      $              fluxsurf_sw,zls)
     766     $              fluxsurf_sw,zls)
    770767
    771768         DO l=1,nlayer
    … …  
    14821479        lmax_th_out(:)=real(lmax_th(:))
    14831480
    1484         call WRITEDIAGFI(ngridmx,'hfx',
    1485      &              'sensible heat flux','W.m-2',
    1486      &                         2,hfx)
    1487 
    14881481        call WRITEDIAGFI(ngridmx,'lmax_th',
    14891482     &              'hauteur du thermique','K',
    … …  
    15341527! THERMALS STUFF 1D
    15351528
    1536         call WRITEDIAGFI(ngridmx,'hfx',
    1537      &              'sensible heat flux','W.m-2',
    1538      &                         0,hfx)
    15391529         if(calltherm) then
    15401530
    … …  
    15801570     &                  tsurf)
    15811571
    1582            call WRITEDIAGFI(ngrid,"fluxsurf_sw",
    1583      &                "Solar radiative flux to surface","W.m-2",0,
    1584      &                fluxsurf_sw_tot)
    1585 
    15861572         call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay)
    15871573         call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev)

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

    r161 r185  
    1       subroutine thermcell_dqupdown(ngrid,nlay,ptimestep,fm0,entr0,  &
     1      subroutine thermcell_dqupdown(ngrid,nlayer,ptimestep,fm0,entr0,  &
    22     &    detr0,masse0,q_therm,dq_therm,ztvd,fm_down,ztv,charvar,lmax)
    33      implicit none
    … …  
    1313!=======================================================================
    1414
     15#include "dimensions.h"
     16#include "dimphys.h"
     17
    1518! ============================ INPUTS ============================
    1619
    17       INTEGER, INTENT(IN) :: ngrid,nlay
     20      INTEGER, INTENT(IN) :: ngrid,nlayer
    1821      REAL, INTENT(IN) :: ptimestep
    19       REAL, INTENT(IN) :: fm0(ngrid,nlay+1)
    20       REAL, INTENT(IN) :: entr0(ngrid,nlay),detr0(ngrid,nlay)
    21       REAL, INTENT(IN) :: q_therm(ngrid,nlay)
    22       REAL, INTENT(IN) :: fm_down(ngrid,nlay+1)
    23       REAL, INTENT(IN) :: ztvd(ngrid,nlay)
    24       REAL, INTENT(IN) :: ztv(ngrid,nlay)
     22      REAL, INTENT(IN) :: fm0(ngridmx,nlayermx+1)
     23      REAL, INTENT(IN) :: entr0(ngridmx,nlayermx),detr0(ngridmx,nlayermx)
     24      REAL, INTENT(IN) :: q_therm(ngridmx,nlayermx)
     25      REAL, INTENT(IN) :: fm_down(ngridmx,nlayermx+1)
     26      REAL, INTENT(IN) :: ztvd(ngridmx,nlayermx)
     27      REAL, INTENT(IN) :: ztv(ngridmx,nlayermx)
    2528      CHARACTER (LEN=20), INTENT(IN) :: charvar
    26       REAL, INTENT(IN) :: masse0(ngrid,nlay)
    27       INTEGER, INTENT(IN) :: lmax(ngrid)
     29      REAL, INTENT(IN) :: masse0(ngridmx,nlayermx)
     30      INTEGER, INTENT(IN) :: lmax(ngridmx)
    2831
    2932! ============================ OUTPUTS ===========================
    3033
    31       REAL, INTENT(OUT) :: dq_therm(ngrid,nlay)
     34      REAL, INTENT(OUT) :: dq_therm(ngridmx,nlayermx)
    3235
    3336! ============================ LOCAL =============================
    3437
    35 !      REAL detr0(ngrid,nlay)
    36       REAL detrd(ngrid,nlay)
    37       REAL entrd(ngrid,nlay)     
    38       REAL fmd(ngrid,nlay+1)
    39       REAL q(ngrid,nlay)
    40       REAL qa(ngrid,nlay)
    41       REAL qd(ngrid,nlay)
     38!      REAL detr0(ngridmx,nlayermx)
     39      REAL detrd(ngridmx,nlayermx)
     40      REAL entrd(ngridmx,nlayermx)     
     41      REAL fmd(ngridmx,nlayermx+1)
     42      REAL q(ngridmx,nlayermx)
     43      REAL qa(ngridmx,nlayermx)
     44      REAL qd(ngridmx,nlayermx)
    4245      INTEGER ig,k
    43       LOGICAL active(ngrid,nlay)
    44       INTEGER lmax_down(ngrid),lmin_down(ngrid)
     46      LOGICAL active(ngridmx,nlayermx)
     47      INTEGER lmax_down(ngridmx),lmin_down(ngridmx)
    4548      INTEGER ncorec
    4649
    … …  
    6467!! ========== DOWNDRAFT TRANSPORT DISABLED FOR NOW ===============
    6568!
    66 !      do ig=1,ngrid
    67 !          if (ztv(ig,nlay)-ztvd(ig,nlay) .gt. 0.5) then
     69!      do ig=1,ngridmx
     70!          if (ztv(ig,nlayermx)-ztvd(ig,nlayermx) .gt. 0.5) then
    6871!            print*,"downdraft non nul derniere couche !!! (dqupdown)"
    6972!          endif
    70 !          detrd(ig,nlay)=0.
    71 !          entrd(ig,nlay)=0.
    72 !      enddo
    73 !
    74 !      do k=nlay-1,1,-1
    75 !          do ig=1,ngrid
     73!          detrd(ig,nlayermx)=0.
     74!          entrd(ig,nlayermx)=0.
     75!      enddo
     76!
     77!      do k=nlayermx-1,1,-1
     78!          do ig=1,ngridmx
    7679!
    7780!          if (ztv(ig,k)-ztvd(ig,k) .gt. 0.0001) then
    … …  
    9396!      lmin_down(:)=1
    9497!
    95 !      do k=1,nlay
    96 !        do ig=1,ngrid
     98!      do k=1,nlayermx
     99!        do ig=1,ngridmx
    97100!        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
    98101!!         if (entrd(ig,k).gt.detrd(ig,k)) then
    … …  
    101104!        enddo
    102105!      enddo
    103 !      do k=nlay,1,-1
    104 !        do ig=1,ngrid
     106!      do k=nlayermx,1,-1
     107!        do ig=1,ngridmx
    105108!        if ((entrd(ig,k).gt.0.) .or. (detrd(ig,k).gt.0.)) then
    106109!!         if (detrd(ig,k).gt.entrd(ig,k)) then
    … …  
    113116!      fmd(:,:)=0.
    114117!       
    115 !      do ig=1,ngrid
     118!      do ig=1,ngridmx
    116119!          if ((lmax_down(ig).gt.1) .and. ((lmax_down(ig)-lmin_down(ig)).gt.1)) then
    117120!!          fmd(ig,lmax_down(ig))=0.
    … …  
    137140!      enddo
    138141!         ncorec=0
    139 !         do k=nlay,2,-1
    140 !         do ig=1,ngrid
     142!         do k=nlayermx,2,-1
     143!         do ig=1,ngridmx
    141144!            if (fmd(ig,k).lt.0.) then
    142145!!               detrd(ig,k)=max(0.,detrd(ig,k)+fmd(ig,k-1))
    … …  
    160163!         print*, lmin_down(:),lmax_down(:)
    161164!
    162 !      do k=2,nlay
    163 !      do ig=1,ngrid
     165!      do k=2,nlayermx
     166!      do ig=1,ngridmx
    164167!          active(ig,k)=(k.ge.lmin_down(ig)).and.(k.le.lmax_down(ig)) &
    165168!     & .and.(((fmd(ig,k)+detrd(ig,k))*ptimestep).gt.1.e-6*masse0(ig,k))
    … …  
    168171!
    169172!
    170 !      do ig=1,ngrid
     173!      do ig=1,ngridmx
    171174!      do k=lmin_down(ig),lmax_down(ig)
    172175!          if(.not.active(ig,k)) then
    … …  
    180183!      endif
    181184!
    182 !!      do ig=1,ngrid
     185!!      do ig=1,ngridmx
    183186!!         active(ig,lmax_down(ig))=(((fmd(ig,lmax_down(ig))+detrd(ig,lmax_down(ig)))* &
    184187!!     &       ptimestep).gt.1.e-8*masse0(ig,lmax_down(ig)))
    185188!!      enddo
    186189!!
    187 !!      do ig=1,ngrid
     190!!      do ig=1,ngridmx
    188191!!          if (lmax_down(ig).gt.1) then
    189192!!            do k=lmax_down(ig)-1,lmin_down(ig),-1
    … …  
    199202!! ========== qa : q in updraft ==================================
    200203!
    201       do k=2,nlay
    202          do ig=1,ngrid
     204      do k=2,nlayermx
     205         do ig=1,ngridmx
    203206            if ((fm0(ig,k+1)+detr0(ig,k))*ptimestep.gt.  &
    204207     &         1.e-5*masse0(ig,k)) then
    … …  
    219222! ========== qd : q in downdraft =================================
    220223!
    221 !      do k=nlay-1,1,-1
    222 !         do ig=1,ngrid
     224!      do k=nlayermx-1,1,-1
     225!         do ig=1,ngridmx
    223226!            if (active(ig,k)) then
    224227!         qd(ig,k)=(fmd(ig,k+1)*qd(ig,k+1)+entrd(ig,k)*q(ig,k))  &
    … …  
    241244! ====== dq ======================================================
    242245
    243       do ig=1,ngrid
     246      do ig=1,ngridmx
    244247         if(active(ig,1)) then
    245248
    … …  
    258261       enddo
    259262     
    260        do k=2,nlay-1
    261          do ig=1, ngrid
     263       do k=2,nlayermx-1
     264         do ig=1, ngridmx
    262265
    263266         if(active(ig,k)) then
    … …  
    281284      enddo
    282285
    283          do ig=1, ngrid
    284 
    285          if(active(ig,nlay)) then
    286 
    287          dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay)+detrd(ig,nlay)*qd(ig,nlay) &
    288       &               +fmd(ig,nlay)*q(ig,nlay-1)   &
    289       &         -entr0(ig,nlay)*q(ig,nlay)-entrd(ig,nlay)*q(ig,nlay)   &
    290       &               -fm0(ig,nlay)*q(ig,nlay)) &
    291       &               *ptimestep/masse0(ig,nlay)
     286         do ig=1, ngridmx
     287
     288         if(active(ig,nlayermx)) then
     289
     290         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx)+detrd(ig,nlayermx)*qd(ig,nlayermx) &
     291      &               +fmd(ig,nlayermx)*q(ig,nlayermx-1)   &
     292      &         -entr0(ig,nlayermx)*q(ig,nlayermx)-entrd(ig,nlayermx)*q(ig,nlayermx)   &
     293      &               -fm0(ig,nlayermx)*q(ig,nlayermx)) &
     294      &               *ptimestep/masse0(ig,nlayermx)
    292295
    293296         else
    294          dq_therm(ig,nlay)=(detr0(ig,nlay)*qa(ig,nlay) &
    295       &             -entr0(ig,nlay)*q(ig,nlay)-fm0(ig,nlay)*q(ig,nlay)) &
    296       &               *ptimestep/masse0(ig,nlay)
     297         dq_therm(ig,nlayermx)=(detr0(ig,nlayermx)*qa(ig,nlayermx) &
     298      &             -entr0(ig,nlayermx)*q(ig,nlayermx)-fm0(ig,nlayermx)*q(ig,nlayermx)) &
     299      &               *ptimestep/masse0(ig,nlayermx)
    297300         endif
    298301         

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

    r173 r185  
    11!
    22!
    3       SUBROUTINE thermcell_main_mars(ngrid,nlay,nq,ptimestep  &
     3      SUBROUTINE thermcell_main_mars(ptimestep  &
    44     &                  ,pplay,pplev,pphi,zlev,zlay  &
    55     &                  ,pu,pv,pt,pq,pq2  &
    66     &                  ,pduadj,pdvadj,pdtadj,pdqadj,pdq2adj  &
    7      &                  ,fm,entr,detr,lmax  &
     7     &                  ,fm,entr,detr,lmax,zmax  &
    88     &                  ,r_aspect &
    99     &                  ,zw2,fraca &
    … …  
    1111     &                  ,buoyancyOut, buoyancyEst)
    1212
    13       USE thermcell,ONLY:RG,RD,RKAPPA
    1413      IMPLICIT NONE
    1514
    … …  
    1716! Mars version of thermcell_main. Author : A Colaitis
    1817!=======================================================================
     18
     19#include "dimensions.h"
     20#include "dimphys.h"
     21#include "comcstfi.h"
     22
    1923! ============== INPUTS ==============
    2024
    21       INTEGER, INTENT(IN) :: ngrid,nlay,nq
    2225      REAL, INTENT(IN) :: ptimestep,r_aspect
    23       REAL, INTENT(IN) :: pt(ngrid,nlay)
    24       REAL, INTENT(IN) :: pu(ngrid,nlay)
    25       REAL, INTENT(IN) :: pv(ngrid,nlay)
    26       REAL, INTENT(IN) :: pq(ngrid,nlay,nq)
    27       REAL, INTENT(IN) :: pq2(ngrid,nlay)
    28       REAL, INTENT(IN) :: pplay(ngrid,nlay)
    29       REAL, INTENT(IN) :: pplev(ngrid,nlay+1)
    30       REAL, INTENT(IN) :: pphi(ngrid,nlay)
    31       REAL, INTENT(IN) :: zlay(ngrid,nlay)
    32       REAL, INTENT(IN) :: zlev(ngrid,nlay+1)
     26      REAL, INTENT(IN) :: pt(ngridmx,nlayermx)
     27      REAL, INTENT(IN) :: pu(ngridmx,nlayermx)
     28      REAL, INTENT(IN) :: pv(ngridmx,nlayermx)
     29      REAL, INTENT(IN) :: pq(ngridmx,nlayermx,nqmx)
     30      REAL, INTENT(IN) :: pq2(ngridmx,nlayermx)
     31      REAL, INTENT(IN) :: pplay(ngridmx,nlayermx)
     32      REAL, INTENT(IN) :: pplev(ngridmx,nlayermx+1)
     33      REAL, INTENT(IN) :: pphi(ngridmx,nlayermx)
     34      REAL, INTENT(IN) :: zlay(ngridmx,nlayermx)
     35      REAL, INTENT(IN) :: zlev(ngridmx,nlayermx+1)
    3336
    3437! ============== OUTPUTS ==============
    3538
    36       REAL, INTENT(OUT) :: pdtadj(ngrid,nlay)
    37       REAL, INTENT(OUT) :: pduadj(ngrid,nlay)
    38       REAL, INTENT(OUT) :: pdvadj(ngrid,nlay)
    39       REAL, INTENT(OUT) :: pdqadj(ngrid,nlay,nq)
    40       REAL, INTENT(OUT) :: pdq2adj(ngrid,nlay)
    41       REAL, INTENT(OUT) :: zw2(ngrid,nlay+1)
     39      REAL, INTENT(OUT) :: pdtadj(ngridmx,nlayermx)
     40      REAL, INTENT(OUT) :: pduadj(ngridmx,nlayermx)
     41      REAL, INTENT(OUT) :: pdvadj(ngridmx,nlayermx)
     42      REAL, INTENT(OUT) :: pdqadj(ngridmx,nlayermx,nqmx)
     43      REAL, INTENT(OUT) :: pdq2adj(ngridmx,nlayermx)
     44      REAL, INTENT(OUT) :: zw2(ngridmx,nlayermx+1)
    4245
    4346! Diagnostics
    44       REAL, INTENT(OUT) :: heatFlux(ngrid,nlay)   ! interface heatflux
    45       REAL, INTENT(OUT) :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft
    46       REAL, INTENT(OUT) :: buoyancyOut(ngrid,nlay)  ! interlayer buoyancy term
    47       REAL, INTENT(OUT) :: buoyancyEst(ngrid,nlay)  ! interlayer estimated buoyancy term
     47      REAL, INTENT(OUT) :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
     48      REAL, INTENT(OUT) :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
     49      REAL, INTENT(OUT) :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
     50      REAL, INTENT(OUT) :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
    4851
    4952! dummy variables when output not needed :
    5053
    51 !      REAL :: heatFlux(ngrid,nlay)   ! interface heatflux
    52 !      REAL :: heatFlux_down(ngrid,nlay) ! interface heat flux from downdraft
    53 !      REAL :: buoyancyOut(ngrid,nlay)  ! interlayer buoyancy term
    54 !      REAL :: buoyancyEst(ngrid,nlay)  ! interlayer estimated buoyancy term
     54!      REAL :: heatFlux(ngridmx,nlayermx)   ! interface heatflux
     55!      REAL :: heatFlux_down(ngridmx,nlayermx) ! interface heat flux from downdraft
     56!      REAL :: buoyancyOut(ngridmx,nlayermx)  ! interlayer buoyancy term
     57!      REAL :: buoyancyEst(ngridmx,nlayermx)  ! interlayer estimated buoyancy term
    5558
    5659
    … …  
    5861
    5962      INTEGER ig,k,l,ll,iq
    60       INTEGER lmax(ngrid),lmin(ngrid),lalim(ngrid)
    61       INTEGER lmix(ngrid)
    62       INTEGER lmix_bis(ngrid)
    63       REAL linter(ngrid)
    64       REAL zmix(ngrid)
    65       REAL zmax(ngrid)
    66       REAL ztva(ngrid,nlay),zw_est(ngrid,nlay+1),ztva_est(ngrid,nlay)
    67       REAL zmax_sec(ngrid)
    68       REAL zh(ngrid,nlay)
    69       REAL zdthladj(ngrid,nlay)
    70       REAL zdthladj_down(ngrid,nlay)
    71       REAL ztvd(ngrid,nlay)
    72       REAL ztv(ngrid,nlay)
    73       REAL zu(ngrid,nlay),zv(ngrid,nlay),zo(ngrid,nlay)
    74       REAL zva(ngrid,nlay)
    75       REAL zua(ngrid,nlay)
    76 
    77       REAL zta(ngrid,nlay)
    78       REAL fraca(ngrid,nlay+1)
    79       REAL q2(ngrid,nlay)
    80       REAL rho(ngrid,nlay),rhobarz(ngrid,nlay),masse(ngrid,nlay)
    81       REAL zpopsk(ngrid,nlay)
    82 
    83       REAL wmax(ngrid)
    84       REAL wmax_sec(ngrid)
    85       REAL fm(ngrid,nlay+1),entr(ngrid,nlay),detr(ngrid,nlay)
    86 
    87       REAL fm_down(ngrid,nlay+1)
    88 
    89       REAL ztla(ngrid,nlay)
    90 
    91       REAL f_star(ngrid,nlay+1),entr_star(ngrid,nlay)
    92       REAL detr_star(ngrid,nlay)
    93       REAL alim_star_tot(ngrid)
    94       REAL alim_star(ngrid,nlay)
    95       REAL alim_star_clos(ngrid,nlay)
    96       REAL f(ngrid)
    97 
    98       REAL teta_th_int(ngrid,nlay)
    99       REAL teta_env_int(ngrid,nlay)
    100       REAL teta_down_int(ngrid,nlay)
     63      INTEGER lmax(ngridmx),lmin(ngridmx),lalim(ngridmx)
     64      INTEGER lmix(ngridmx)
     65      INTEGER lmix_bis(ngridmx)
     66      REAL linter(ngridmx)
     67      REAL zmix(ngridmx)
     68      REAL zmax(ngridmx)
     69      REAL ztva(ngridmx,nlayermx),zw_est(ngridmx,nlayermx+1),ztva_est(ngridmx,nlayermx)
     70      REAL zmax_sec(ngridmx)
     71      REAL zh(ngridmx,nlayermx)
     72      REAL zdthladj(ngridmx,nlayermx)
     73      REAL zdthladj_down(ngridmx,nlayermx)
     74      REAL ztvd(ngridmx,nlayermx)
     75      REAL ztv(ngridmx,nlayermx)
     76      REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx),zo(ngridmx,nlayermx)
     77      REAL zva(ngridmx,nlayermx)
     78      REAL zua(ngridmx,nlayermx)
     79
     80      REAL zta(ngridmx,nlayermx)
     81      REAL fraca(ngridmx,nlayermx+1)
     82      REAL q2(ngridmx,nlayermx)
     83      REAL rho(ngridmx,nlayermx),rhobarz(ngridmx,nlayermx),masse(ngridmx,nlayermx)
     84      REAL zpopsk(ngridmx,nlayermx)
     85
     86      REAL wmax(ngridmx)
     87      REAL wmax_sec(ngridmx)
     88      REAL fm(ngridmx,nlayermx+1),entr(ngridmx,nlayermx),detr(ngridmx,nlayermx)
     89
     90      REAL fm_down(ngridmx,nlayermx+1)
     91
     92      REAL ztla(ngridmx,nlayermx)
     93
     94      REAL f_star(ngridmx,nlayermx+1),entr_star(ngridmx,nlayermx)
     95      REAL detr_star(ngridmx,nlayermx)
     96      REAL alim_star_tot(ngridmx)
     97      REAL alim_star(ngridmx,nlayermx)
     98      REAL alim_star_clos(ngridmx,nlayermx)
     99      REAL f(ngridmx)
     100
     101      REAL teta_th_int(ngridmx,nlayermx)
     102      REAL teta_env_int(ngridmx,nlayermx)
     103      REAL teta_down_int(ngridmx,nlayermx)
    101104
    102105      CHARACTER (LEN=20) :: modname
    … …  
    105108! ============= PLUME VARIABLES ============
    106109
    107       REAL w_est(ngrid,nlay+1)
    108       REAL wa_moy(ngrid,nlay+1)
    109       REAL wmaxa(ngrid)
    110       REAL zdz,zbuoy(ngrid,nlay),zw2m
    111       LOGICAL active(ngrid),activetmp(ngrid)
     110      REAL w_est(ngridmx,nlayermx+1)
     111      REAL wa_moy(ngridmx,nlayermx+1)
     112      REAL wmaxa(ngridmx)
     113      REAL zdz,zbuoy(ngridmx,nlayermx),zw2m
     114      LOGICAL active(ngridmx),activetmp(ngridmx)
     115      REAL a1,b1,ae,be,ad,bd
    112116
    113117! ==========================================
    … …  
    115119! ============= HEIGHT VARIABLES ===========
    116120
    117       REAL num(ngrid)
    118       REAL denom(ngrid)
    119       REAL zlevinter(ngrid)
     121      REAL num(ngridmx)
     122      REAL denom(ngridmx)
     123      REAL zlevinter(ngridmx)
    120124
    121125! =========================================
    … …  
    123127! ============= DRY VARIABLES =============
    124128
    125        REAL zw2_dry(ngrid,nlay+1)
    126        REAL f_star_dry(ngrid,nlay+1)
    127        REAL ztva_dry(ngrid,nlay+1)
    128        REAL wmaxa_dry(ngrid)
    129        REAL wa_moy_dry(ngrid,nlay+1)
    130        REAL linter_dry(ngrid),zlevinter_dry(ngrid)
    131        INTEGER lmix_dry(ngrid),lmax_dry(ngrid)
     129       REAL zw2_dry(ngridmx,nlayermx+1)
     130       REAL f_star_dry(ngridmx,nlayermx+1)
     131       REAL ztva_dry(ngridmx,nlayermx+1)
     132       REAL wmaxa_dry(ngridmx)
     133       REAL wa_moy_dry(ngridmx,nlayermx+1)
     134       REAL linter_dry(ngridmx),zlevinter_dry(ngridmx)
     135       INTEGER lmix_dry(ngridmx),lmax_dry(ngridmx)
    132136
    133137! =========================================
    … …  
    135139! ============= CLOSURE VARIABLES =========
    136140
    137       REAL zdenom(ngrid)
    138       REAL alim_star2(ngrid)
    139       REAL alim_star_tot_clos(ngrid)
     141      REAL zdenom(ngridmx)
     142      REAL alim_star2(ngridmx)
     143      REAL alim_star_tot_clos(ngridmx)
    140144      INTEGER llmax
    141145
    … …  
    161165!               not used for now
    162166
    163       REAL qa(ngrid,nlay),detr_dv2(ngrid,nlay),zf,zf2
    164       REAL wvd(ngrid,nlay+1),wud(ngrid,nlay+1)
    165       REAL gamma0(ngrid,nlay+1),gamma(ngrid,nlay+1)
    166       REAL ue(ngrid,nlay),ve(ngrid,nlay)
    167       LOGICAL ltherm(ngrid,nlay)
    168       REAL dua(ngrid,nlay),dva(ngrid,nlay)
     167      REAL qa(ngridmx,nlayermx),detr_dv2(ngridmx,nlayermx),zf,zf2
     168      REAL wvd(ngridmx,nlayermx+1),wud(ngridmx,nlayermx+1)
     169      REAL gamma0(ngridmx,nlayermx+1),gamma(ngridmx,nlayermx+1)
     170      REAL ue(ngridmx,nlayermx),ve(ngridmx,nlayermx)
     171      LOGICAL ltherm(ngridmx,nlayermx)
     172      REAL dua(ngridmx,nlayermx),dva(ngridmx,nlayermx)
    169173      INTEGER iter
    170174      INTEGER nlarga0
    … …  
    202206!-----------------------------------------------------------------------
    203207
    204 !      do l=2,nlay
    205 !         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/RG
     208!      do l=2,nlayermx
     209!         zlev(:,l)=0.5*(pphi(:,l)+pphi(:,l-1))/g
    206210!      enddo
    207211!         zlev(:,1)=0.
    208 !         zlev(:,nlay+1)=(2.*pphi(:,nlay)-pphi(:,nlay-1))/RG
    209 
    210 !         zlay(:,:)=pphi(:,:)/RG
     212!         zlev(:,nlayermx+1)=(2.*pphi(:,nlayermx)-pphi(:,nlayermx-1))/g
     213
     214!         zlay(:,:)=pphi(:,:)/g
    211215!-----------------------------------------------------------------------
    212216!   Calcul des densites
    213217!-----------------------------------------------------------------------
    214218
    215       rho(:,:)=pplay(:,:)/(RD*pt(:,:))
     219      rho(:,:)=pplay(:,:)/(r*pt(:,:))
    216220
    217221      rhobarz(:,1)=rho(:,1)
    218222
    219       do l=2,nlay
     223      do l=2,nlayermx
    220224         rhobarz(:,l)=0.5*(rho(:,l)+rho(:,l-1))
    221225      enddo
    222226
    223227!calcul de la masse
    224       do l=1,nlay
    225          masse(:,l)=(pplev(:,l)-pplev(:,l+1))/RG
     228      do l=1,nlayermx
     229         masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g
    226230      enddo
    227231
    … …  
    309313      wa_moy(:,:)=0.
    310314      linter(:)=1.
     315
     316!      a1=2.5 ; b1=0.0015 ; ae=0.045 ; be = 0.6  ! svn baseline
     317
     318      a1=1.60226 ; b1=0.0006 ; ae=0.0454 ; be = 0.57 !improved fits
     319      ad = 0.0005114  ; bd = -0.662
     320
    311321! Initialisation des variables entieres
    312322      lmix(:)=1
    … …  
    320330!-------------------------------------------------------------------------
    321331      active(:)=ztv(:,1)>ztv(:,2)
    322           do ig=1,ngrid
     332          do ig=1,ngridmx
    323333            if (ztv(ig,1)>=(ztv(ig,2))) then
    324334               alim_star(ig,1)=MAX((ztv(ig,1)-ztv(ig,2)),0.)  &
    … …  
    330340         enddo
    331341
    332        do l=2,nlay-1
    333          do ig=1,ngrid
    334             if (ztv(ig,l)>(ztv(ig,l+1)+0.5) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then
     342       do l=2,nlayermx-1
     343!        do l=2,4
     344         do ig=1,ngridmx
     345           if (ztv(ig,l)>(ztv(ig,l+1)+0.) .and. ztv(ig,1)>=ztv(ig,l) .and. (alim_star(ig,l-1) .ne. 0.)) then
    335346               alim_star(ig,l)=MAX((ztv(ig,l)-ztv(ig,l+1)),0.)  &
    336347     &                       *sqrt(zlev(ig,l+1))
    … …  
    341352        enddo
    342353      enddo
    343       do l=1,nlay
    344          do ig=1,ngrid
     354      do l=1,nlayermx
     355         do ig=1,ngridmx
    345356            if (alim_star_tot(ig) > 1.e-10 ) then
    346357               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    … …  
    350361
    351362      alim_star_tot(:)=1.
    352       if(alim_star(1,1) .ne. 0.) then
    353       print*, alim_star(:,:)
    354       endif
     363!      if(alim_star(1,1) .ne. 0.) then
     364!      print*, 'lalim star' ,lalim(1)
     365!      endif
    355366
    356367!------------------------------------------------------------------------------
    … …  
    362373!------------------------------------------------------------------------------
    363374
    364       do ig=1,ngrid
     375      do ig=1,ngridmx
    365376! Le panache va prendre au debut les caracteristiques de l'air contenu
    366377! dans cette couche.
    … …  
    371382          f_star(ig,1)=0.
    372383          f_star(ig,2)=alim_star(ig,1)
    373           zw2(ig,2)=2.*RG*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
     384          zw2(ig,2)=2.*g*(ztv(ig,1)-ztv(ig,2))/ztv(ig,2)  &
    374385      &                     *(zlev(ig,2)-zlev(ig,1))  &
    375386      &                     *0.4*pphi(ig,1)/(pphi(ig,2)-pphi(ig,1))
    … …  
    383394!boucle de calcul de la vitesse verticale dans le thermique
    384395!==============================================================================
    385       do l=2,nlay-1
     396      do l=2,nlayermx-1
    386397!==============================================================================
    387398
    … …  
    389400! On decide si le thermique est encore actif ou non
    390401! AFaire : Il faut sans doute ajouter entr_star a alim_star dans ce test
    391           do ig=1,ngrid
     402          do ig=1,ngridmx
    392403             active(ig)=active(ig) &
    393404      &                 .and. zw2(ig,l)>1.e-10 &
    … …  
    405416!---------------------------------------------------------------------------
    406417
    407           do ig=1,ngrid
     418          do ig=1,ngridmx
    408419             if(active(ig)) then
    409420
    … …  
    417428
    418429                zdz=zlev(ig,l+1)-zlev(ig,l)
    419                 zbuoy(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    420                 if (((2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
    421                 w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015 &
    422      & -2.*zdz*w_est(ig,l)*0.045*(2.5*zbuoy(ig,l)/w_est(ig,l)-0.0015)**0.6)
     430                zbuoy(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     431
     432                if (((a1*zbuoy(ig,l)/w_est(ig,l)-b1) .gt. 0.) .and. (w_est(ig,l) .ne. 0.)) then
     433                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1 &
     434     & -2.*zdz*w_est(ig,l)*ae*(a1*zbuoy(ig,l)/w_est(ig,l)-b1)**be)
    423435                else
    424                 w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*0.0015)
     436                w_est(ig,l+1)=Max(0.0001,w_est(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*w_est(ig,l)*b1)
    425437                endif
    426438                if (w_est(ig,l+1).lt.0.) then
    … …  
    434446!-------------------------------------------------
    435447
    436       do ig=1,ngrid
     448      do ig=1,ngridmx
    437449        if (active(ig)) then
    438450
    439451          zw2m=w_est(ig,l+1)
    440           if((2.5*(zbuoy(ig,l)/zw2m)-0.0015) .gt. 0.) then
     452
     453          if((a1*(zbuoy(ig,l)/zw2m)-b1) .gt. 0.) then
    441454          entr_star(ig,l)=f_star(ig,l)*zdz*  &
    442 !       &   0.0118*((zbuoy(ig,l)/zw2m)/0.043)**(1./1.65)
    443 !       &   0.0090*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.90)
    444 !       &   0.0120*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.6)
    445         &   MAX(0.,0.045*(2.5*(zbuoy(ig,l)/zw2m)-0.0015)**0.6)
     455        &   MAX(0.,ae*(a1*(zbuoy(ig,l)/zw2m)-b1)**be)
    446456          else
    447457          entr_star(ig,l)=0.
    448458          endif
     459
    449460          if(zbuoy(ig,l) .gt. 0.) then
    450461             if(l .lt. lalim(ig)) then
    451                  detr_star(ig,l)=0.
     462                detr_star(ig,l)=0.
    452463             else
     464
    453465!                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    454466!              &     0.0105*((zbuoy(ig,l)/zw2m)/0.048)**(1./1.7)
    … …  
    462474! new param from continuity eq with a fit on dfdz
    463475                 detr_star(ig,l) = f_star(ig,l)*zdz*              &
    464                &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)
     476            &  ad
     477
     478!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)   !svn baseline
     479!                &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0008)     
    465480
    466481!              &     0.014*((zbuoy(ig,l)/zw2m)/0.05)**(1./1.35)
    … …  
    471486          else
    472487          detr_star(ig,l)=f_star(ig,l)*zdz*                        &
    473                &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)
     488                &     bd*zbuoy(ig,l)/zw2m
     489
     490!               &     MAX(0.,-0.38*zbuoy(ig,l)/zw2m+0.0005)      !svn baseline
    474491
    475492!              &  *5.*(-afact*zbetalpha*zbuoy(ig,l)/zw2m)
    … …  
    506523
    507524      activetmp(:)=active(:) .and. f_star(:,l+1)>1.e-10
    508       do ig=1,ngrid
     525      do ig=1,ngridmx
    509526       if (activetmp(ig)) then
    510527
    … …  
    516533      enddo
    517534
    518       do ig=1,ngrid
     535      do ig=1,ngridmx
    519536      if (activetmp(ig)) then
    520537           ztva(ig,l) = ztla(ig,l)
    521            zbuoy(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    522 
    523            if (((2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
    524            zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-         &
    525      & 2.*zdz*zw2(ig,l)*0.0015-2.*zdz*zw2(ig,l)*0.045*(2.5*zbuoy(ig,l)/zw2(ig,l)-0.0015)**0.6)
     538           zbuoy(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     539
     540           if (((a1*zbuoy(ig,l)/zw2(ig,l)-b1) .gt. 0.) .and. (zw2(ig,l) .ne. 0.) ) then
     541           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-         &
     542     & 2.*zdz*zw2(ig,l)*b1-2.*zdz*zw2(ig,l)*ae*(a1*zbuoy(ig,l)/zw2(ig,l)-b1)**be)
    526543           else
    527            zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*2.5*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*0.0015)
     544           zw2(ig,l+1)=Max(0.,zw2(ig,l)+2.*zdz*a1*zbuoy(ig,l)-2.*zdz*zw2(ig,l)*b1)
    528545           endif
    529546      endif
    … …  
    534551!---------------------------------------------------------------------------
    535552
    536       do ig=1,ngrid
     553      do ig=1,ngridmx
    537554            if (zw2(ig,l+1)>0. .and. zw2(ig,l+1).lt.1.e-10) then
    538555             print*,'On tombe sur le cas particulier de thermcell_plume'
    … …  
    562579
    563580!on recalcule alim_star_tot
    564        do ig=1,ngrid
     581       do ig=1,ngridmx
    565582          alim_star_tot(ig)=0.
    566583       enddo
    567        do ig=1,ngrid
     584       do ig=1,ngridmx
    568585          do l=1,lalim(ig)-1
    569586          alim_star_tot(ig)=alim_star_tot(ig)+alim_star(ig,l)
    … …  
    571588       enddo
    572589
    573       do l=1,nlay
    574          do ig=1,ngrid
     590      do l=1,nlayermx
     591         do ig=1,ngridmx
    575592            if (alim_star_tot(ig) > 1.e-10 ) then
    576593               alim_star(ig,l)=alim_star(ig,l)/alim_star_tot(ig)
    … …  
    595612
    596613!calcul de la hauteur max du thermique
    597       do ig=1,ngrid
     614      do ig=1,ngridmx
    598615         lmax(ig)=lalim(ig)
    599616      enddo
    600       do ig=1,ngrid
    601          do l=nlay,lalim(ig)+1,-1
     617      do ig=1,ngridmx
     618         do l=nlayermx,lalim(ig)+1,-1
    602619            if (zw2(ig,l).le.1.e-10) then
    603620               lmax(ig)=l-1
    … …  
    608625! On traite le cas particulier qu'il faudrait éviter ou le thermique
    609626! atteind le haut du modele ...
    610       do ig=1,ngrid
    611       if ( zw2(ig,nlay) > 1.e-10 ) then
     627      do ig=1,ngridmx
     628      if ( zw2(ig,nlayermx) > 1.e-10 ) then
    612629          print*,'WARNING !!!!! W2 thermiques non nul derniere couche '
    613           lmax(ig)=nlay
     630          lmax(ig)=nlayermx
    614631      endif
    615632      enddo
    616633
    617634! pas de thermique si couche 1 stable
    618 !      do ig=1,ngrid
     635!      do ig=1,ngridmx
    619636!         if (lmin(ig).gt.1) then
    620637!             lmax(ig)=1
    … …  
    625642!
    626643! Determination de zw2 max
    627       do ig=1,ngrid
     644      do ig=1,ngridmx
    628645         wmax(ig)=0.
    629646      enddo
    630647
    631       do l=1,nlay
    632          do ig=1,ngrid
     648      do l=1,nlayermx
     649         do ig=1,ngridmx
    633650            if (l.le.lmax(ig)) then
    634651                if (zw2(ig,l).lt.0.)then
    … …  
    643660      enddo
    644661!   Longueur caracteristique correspondant a la hauteur des thermiques.
    645       do  ig=1,ngrid
     662      do  ig=1,ngridmx
    646663         zmax(ig)=0.
    647664         zlevinter(ig)=zlev(ig,1)
    … …  
    650667         num(:)=0.
    651668         denom(:)=0.
    652          do ig=1,ngrid
    653           do l=1,nlay
     669         do ig=1,ngridmx
     670          do l=1,nlayermx
    654671             num(ig)=num(ig)+zw2(ig,l)*zlev(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    655672             denom(ig)=denom(ig)+zw2(ig,l)*(zlev(ig,l+1)-zlev(ig,l))
    656673          enddo
    657674       enddo
    658        do ig=1,ngrid
     675       do ig=1,ngridmx
    659676       if (denom(ig).gt.1.e-10) then
    660677          zmax(ig)=2.*num(ig)/denom(ig)
    … …  
    663680
    664681! def de  zmix continu (profil parabolique des vitesses)
    665       do ig=1,ngrid
     682      do ig=1,ngridmx
    666683           if (lmix(ig).gt.1) then
    667684! test
    … …  
    673690!
    674691            zmix(ig)=((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
    675      &        *((zlev(ig,lmix(ig)))**2-(zlev(ig,lmix(ig)+1))**2)  &
     692     &        *((zlev(ig,lmix(ig))*zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)*zlev(ig,lmix(ig)+1)))  &
    676693     &        -(zw2(ig,lmix(ig))-zw2(ig,lmix(ig)+1))  &
    677      &        *((zlev(ig,lmix(ig)-1))**2-(zlev(ig,lmix(ig)))**2))  &
     694     &        *((zlev(ig,lmix(ig)-1)*zlev(ig,lmix(ig)-1))-(zlev(ig,lmix(ig))*zlev(ig,lmix(ig)))))  &
    678695     &        /(2.*((zw2(ig,lmix(ig)-1)-zw2(ig,lmix(ig)))  &
    679696     &        *((zlev(ig,lmix(ig)))-(zlev(ig,lmix(ig)+1)))  &
    … …  
    694711!
    695712! calcul du nouveau lmix correspondant
    696       do ig=1,ngrid
    697          do l=1,nlay
     713      do ig=1,ngridmx
     714         do l=1,nlayermx
    698715            if (zmix(ig).ge.zlev(ig,l).and.  &
    699716     &          zmix(ig).lt.zlev(ig,l+1)) then
    … …  
    710727! ===========================================================================
    711728
    712 
    713 ! ===========================================================================
    714 ! ================= DRY =====================================================
    715 ! ===========================================================================
    716 
    717 !initialisations
    718        do ig=1,ngrid
    719           do l=1,nlay+1
    720              zw2_dry(ig,l)=0.
    721              wa_moy_dry(ig,l)=0.
    722           enddo
    723        enddo
    724        do ig=1,ngrid
    725           do l=1,nlay
    726              ztva_dry(ig,l)=ztv(ig,l)
    727           enddo
    728        enddo
    729        do ig=1,ngrid
    730           wmax_sec(ig)=0.
    731           wmaxa_dry(ig)=0.
    732        enddo
    733 !calcul de la vitesse a partir de la CAPE en melangeant thetav
    734 
    735 
    736 ! Calcul des F^*, integrale verticale de E^*
    737        f_star_dry(:,1)=0.
    738        do l=1,nlay
    739           f_star_dry(:,l+1)=f_star_dry(:,l)+alim_star(:,l)
    740        enddo
    741 
    742 ! niveau (reel) auquel zw2 s'annule FH :n'etait pas initialise
    743        linter_dry(:)=0.
    744 
    745 ! couche la plus haute concernee par le thermique.
    746        lmax_dry(:)=1
    747 
    748 ! Le niveau linter est une variable continue qui se trouve dans la couche
    749 ! lmax
    750 
    751        do l=1,nlay-2
    752          do ig=1,ngrid
    753             if (l.eq.lmin(ig).and.lalim(ig).gt.1) then
    754 
    755 !------------------------------------------------------------------------
    756 !  Calcul de la vitesse en haut de la premiere couche instable.
    757 !  Premiere couche du panache thermique
    758 !------------------------------------------------------------------------
    759 
    760              zw2_dry(ig,l+1)=2.*RG*(ztv(ig,l)-ztv(ig,l+1))/ztv(ig,l+1)  &
    761      &                     *(zlev(ig,l+1)-zlev(ig,l))  &
    762      &                     *0.4*pphi(ig,l)/(pphi(ig,l+1)-pphi(ig,l))
    763 !------------------------------------------------------------------------
    764 ! Tant que la vitesse en bas de la couche et la somme du flux de masse
    765 ! et de l'entrainement (c'est a dire le flux de masse en haut) sont
    766 ! positifs, on calcul
    767 ! 1. le flux de masse en haut  f_star(ig,l+1)
    768 ! 2. la temperature potentielle virtuelle dans la couche ztva(ig,l)
    769 ! 3. la vitesse au carré en haut zw2(ig,l+1)
    770 !------------------------------------------------------------------------
    771 
    772             else if (zw2_dry(ig,l).ge.1e-10) then
    773 
    774        ztva_dry(ig,l)=(f_star_dry(ig,l)*ztva_dry(ig,l-1)+alim_star(ig,l) &
    775      &                    *ztv(ig,l))/f_star_dry(ig,l+1)
    776       zw2_dry(ig,l+1)=zw2_dry(ig,l)*(f_star_dry(ig,l)/f_star_dry(ig,l+1))**2+  &
    777      &                     2.*RG*(ztva_dry(ig,l)-ztv(ig,l))/ztv(ig,l)  &
    778      &                     *(zlev(ig,l+1)-zlev(ig,l))
    779             endif
    780 ! determination de zmax continu par interpolation lineaire
    781 !------------------------------------------------------------------------
    782 
    783             if (zw2_dry(ig,l+1)>0. .and. zw2_dry(ig,l+1).lt.1.e-10) then
    784 !               stop'On tombe sur le cas particulier de thermcell_dry'
    785 !               print*,'On tombe sur le cas particulier de thermcell_dry'
    786                 zw2_dry(ig,l+1)=0.
    787                 linter_dry(ig)=l+1
    788                 lmax_dry(ig)=l
    789             endif
    790 
    791             if (zw2_dry(ig,l+1).lt.0.) then
    792                linter_dry(ig)=(l*(zw2_dry(ig,l+1)-zw2_dry(ig,l))  &
    793      &           -zw2_dry(ig,l))/(zw2_dry(ig,l+1)-zw2_dry(ig,l))
    794                zw2_dry(ig,l+1)=0.
    795                lmax_dry(ig)=l
    796             endif
    797 
    798                wa_moy_dry(ig,l+1)=sqrt(zw2_dry(ig,l+1))
    799 
    800             if (wa_moy_dry(ig,l+1).gt.wmaxa_dry(ig)) then
    801 !   lmix est le niveau de la couche ou w (wa_moy) est maximum
    802                lmix_dry(ig)=l+1
    803                wmaxa_dry(ig)=wa_moy_dry(ig,l+1)
    804             endif
    805          enddo
    806       enddo
    807 !
    808 ! Determination de zw2 max
    809       do ig=1,ngrid
    810          wmax_sec(ig)=0.
    811       enddo
    812       do l=1,nlay
    813          do ig=1,ngrid
    814             if (l.le.lmax_dry(ig)) then
    815                 zw2_dry(ig,l)=sqrt(zw2_dry(ig,l))
    816                 wmax_sec(ig)=max(wmax_sec(ig),zw2_dry(ig,l))
    817             else
    818                  zw2_dry(ig,l)=0.
    819             endif
    820           enddo
    821       enddo
    822 
    823 !   Longueur caracteristique correspondant a la hauteur des thermiques.
    824       do  ig=1,ngrid
    825          zmax_sec(ig)=0.
    826          zlevinter_dry(ig)=zlev(ig,1)
    827       enddo
    828       do  ig=1,ngrid
    829 ! calcul de zlevinter
    830           zlevinter_dry(ig)=zlev(ig,lmax_dry(ig)) + &
    831      &    (linter_dry(ig)-lmax_dry(ig))*(zlev(ig,lmax_dry(ig)+1)-zlev(ig,lmax_dry(ig)))
    832       zmax_sec(ig)=max(zmax_sec(ig),zlevinter_dry(ig)-zlev(ig,lmin(ig)))
    833       enddo
    834 
    835 ! ===========================================================================
    836 ! ============= FIN DRY =====================================================
    837 ! ===========================================================================
    838 
    839 
    840729! Choix de la fonction d'alimentation utilisee pour la fermeture.
    841730
    … …  
    857746! Indice vertical max (max de lalim) atteint par les thermiques sur le domaine
    858747      llmax=1
    859       do ig=1,ngrid
     748      do ig=1,ngridmx
    860749         if (lalim(ig)>llmax) llmax=lalim(ig)
    861750      enddo
    … …  
    865754!   alim_star^2/(rho dz)
    866755      do k=1,llmax-1
    867          do ig=1,ngrid
     756         do ig=1,ngridmx
    868757            if (k<lalim(ig)) then
    869                alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)**2  &
     758         alim_star2(ig)=alim_star2(ig)+alim_star_clos(ig,k)*alim_star_clos(ig,k)  &
    870759      &                    /(rho(ig,k)*(zlev(ig,k+1)-zlev(ig,k)))
    871760         alim_star_tot_clos(ig)=alim_star_tot_clos(ig)+alim_star_clos(ig,k)
    … …  
    873762         enddo
    874763      enddo
    875 
     764 
    876765! WARNING : MARS MODIF : we have added 2. : ratio of wmax/vmoy
    877766! True ratio is 3.5 but wetake into account the vmoy is the one alimentating
    878767! the thermal, so there are vs=0 into the vmoy... the true vmoy is lower. (a la louche)
    879768! And r_aspect has been changed from 2 to 1.5 from observations
    880       do ig=1,ngrid
     769      do ig=1,ngridmx
    881770         if (alim_star2(ig)>1.e-10) then
    882             f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/  &
    883       &     (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
     771!            f(ig)=wmax_sec(ig)*alim_star_tot_clos(ig)/  &
     772!      &     (max(500.,zmax_sec(ig))*r_aspect*alim_star2(ig))
     773             f(ig)=wmax(ig)*alim_star_tot_clos(ig)/  &
     774      &     (max(500.,zmax(ig))*r_aspect*alim_star2(ig))
     775
    884776         endif
    885777      enddo
    … …  
    915807!-------------------------------------------------------------------------
    916808
    917       do l=1,nlay
     809      do l=1,nlayermx
    918810         entr(:,l)=f(:)*(entr_star(:,l)+alim_star(:,l))
    919811         detr(:,l)=f(:)*detr_star(:,l)
    920812      enddo
    921813
    922       do l=1,nlay
    923          do ig=1,ngrid
     814      do l=1,nlayermx
     815         do ig=1,ngridmx
    924816            if (l.lt.lmax(ig)) then
    925817               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
    … …  
    937829! autres corrections.
    938830
    939       do l=1,nlay
    940          do ig=1,ngrid
     831      do l=1,nlayermx
     832         do ig=1,ngridmx
    941833            if (detr(ig,l).gt.fm(ig,l)) then
    942834               ncorecfm8=ncorecfm8+1
    … …  
    954846!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    955847
    956       do l=1,nlay
    957 
    958          do ig=1,ngrid
     848      do l=1,nlayermx
     849
     850         do ig=1,ngridmx
    959851            if (l.lt.lmax(ig)) then
    960852               fm(ig,l+1)=fm(ig,l)+entr(ig,l)-detr(ig,l)
    … …  
    972864!-------------------------------------------------------------------------
    973865
    974          do ig=1,ngrid
     866         do ig=1,ngridmx
    975867            if (fm(ig,l+1).lt.0.) then
    976868              print*,'fm1<0',l+1,lmax(ig),fm(ig,l+1)
    … …  
    988880!-------------------------------------------------------------------------
    989881!      if (iflag_thermals_optflux==0) then
    990 !         do ig=1,ngrid
     882!         do ig=1,ngridmx
    991883!          if (l.ge.lalim(ig).and.l.le.lmax(ig) &
    992884!     &    .and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) then
    … …  
    1007899!-------------------------------------------------------------------------
    1008900!      if (iflag_thermals_optflux==0) then
    1009 !         do ig=1,ngrid
     901!         do ig=1,ngridmx
    1010902!            if ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
    1011903!              f_old=fm(ig,l+1)
    … …  
    1021913!-------------------------------------------------------------------------
    1022914
    1023          do ig=1,ngrid
     915         do ig=1,ngridmx
    1024916            if (detr(ig,l).gt.fm(ig,l)) then
    1025917               ncorecfm6=ncorecfm6+1
    … …  
    1042934!-------------------------------------------------------------------------
    1043935
    1044          do ig=1,ngrid
     936         do ig=1,ngridmx
    1045937            if (fm(ig,l+1).lt.0.) then
    1046938               detr(ig,l)=detr(ig,l)+fm(ig,l+1)
    … …  
    1067959!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1068960
    1069         do ig=1,ngrid
     961        do ig=1,ngridmx
    1070962           if (zw2(ig,l+1).gt.1.e-10) then
    1071963           zfm=rhobarz(ig,l+1)*zw2(ig,l+1)*alphamax
    … …  
    1088980!-----------------------------------------------------------------------
    1089981
    1090       do l=1,nlay-1
    1091          do ig=1,ngrid
     982      do l=1,nlayermx-1
     983         do ig=1,ngridmx
    1092984            eee0=entr(ig,l)
    1093985            ddd0=detr(ig,l)
    … …  
    11191011!              ddd=detr(ig)-entre
    11201012!on s assure que tout s annule bien en zmax
    1121       do ig=1,ngrid
     1013      do ig=1,ngridmx
    11221014         fm(ig,lmax(ig)+1)=0.
    11231015         entr(ig,lmax(ig))=0.
    … …  
    11301022
    11311023!IM 090508 beg
    1132       if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngrid/4. ) then
     1024      if (ncorecfm1+ncorecfm2+ncorecfm3+ncorecfm4+ncorecfm5+ncorecalpha > ngridmx/4. ) then
    11331025         print*,'thermcell warning : large number of corrections'
    11341026         print*,'PB thermcell : on a du coriger ',ncorecfm1,'x fm1',&
    … …  
    11601052      zdthladj(:,:)=0.
    11611053
    1162       do ig=1,ngrid
     1054      do ig=1,ngridmx
    11631055         if(lmax(ig) .gt. 1) then
    11641056         do k=1,lmax(ig)
    … …  
    11811073      ztvd(:,:)=ztv(:,:)
    11821074      fm_down(:,:)=0.
    1183       do ig=1,ngrid
     1075      do ig=1,ngridmx
    11841076         if (lmax(ig) .gt. 1) then
    11851077         do l=1,lmax(ig)
    … …  
    12191111       zdthladj_down(:,:)=0.
    12201112
    1221       do ig=1,ngrid
     1113      do ig=1,ngridmx
    12221114         if(lmax(ig) .gt. 1) then
    12231115         do k=1,lmax(ig)
    … …  
    12371129! Calcul de la fraction de l'ascendance
    12381130!------------------------------------------------------------------
    1239       do ig=1,ngrid
     1131      do ig=1,ngridmx
    12401132         fraca(ig,1)=0.
    1241          fraca(ig,nlay+1)=0.
    1242       enddo
    1243       do l=2,nlay
    1244          do ig=1,ngrid
     1133         fraca(ig,nlayermx+1)=0.
     1134      enddo
     1135      do l=2,nlayermx
     1136         do ig=1,ngridmx
    12451137            if (zw2(ig,l).gt.1.e-10) then
    12461138            fraca(ig,l)=fm(ig,l)/(rhobarz(ig,l)*zw2(ig,l))
    … …  
    12731165      nlarga0=0.
    12741166
    1275       do k=1,nlay
    1276          do ig=1,ngrid
     1167      do k=1,nlayermx
     1168         do ig=1,ngridmx
    12771169            detr_dv2(ig,k)=fm(ig,k)-fm(ig,k+1)+entr(ig,k)
    12781170         enddo
    … …  
    12801172
    12811173!   calcul de la valeur dans les ascendances
    1282       do ig=1,ngrid
     1174      do ig=1,ngridmx
    12831175         zua(ig,1)=zu(ig,1)
    12841176         zva(ig,1)=zv(ig,1)
    … …  
    12871179      enddo
    12881180
    1289       gamma(1:ngrid,1)=0.
    1290       do k=2,nlay
    1291          do ig=1,ngrid
     1181      gamma(1:ngridmx,1)=0.
     1182      do k=2,nlayermx
     1183         do ig=1,ngridmx
    12921184            ltherm(ig,k)=(fm(ig,k+1)+detr_dv2(ig,k))*ptimestep > 1.e-5*masse(ig,k)
    12931185            if(ltherm(ig,k).and.zmax(ig)>0.) then
    … …  
    13051197      gamma(:,:)=0.
    13061198
    1307       do k=2,nlay
    1308 
    1309          do ig=1,ngrid
     1199      do k=2,nlayermx
     1200
     1201         do ig=1,ngridmx
    13101202
    13111203            if (ltherm(ig,k)) then
    … …  
    13271219
    13281220      do iter=1,5
    1329          do ig=1,ngrid
     1221         do ig=1,ngridmx
    13301222! Pour memoire : calcul prenant en compte la fraction reelle
    13311223!              zf=0.5*(fraca(ig,k)+fraca(ig,k+1))
    … …  
    13631255      enddo
    13641256
    1365       enddo ! k=2,nlay
     1257      enddo ! k=2,nlayermx
    13661258
    13671259! Calcul du flux vertical de moment dans l'environnement.
    13681260!---------------------------------------------------------
    1369       do k=2,nlay
    1370          do ig=1,ngrid
     1261      do k=2,nlayermx
     1262         do ig=1,ngridmx
    13711263            wud(ig,k)=fm(ig,k)*ue(ig,k)
    13721264            wvd(ig,k)=fm(ig,k)*ve(ig,k)
    13731265         enddo
    13741266      enddo
    1375       do ig=1,ngrid
     1267      do ig=1,ngridmx
    13761268         wud(ig,1)=0.
    1377          wud(ig,nlay+1)=0.
     1269         wud(ig,nlayermx+1)=0.
    13781270         wvd(ig,1)=0.
    1379          wvd(ig,nlay+1)=0.
     1271         wvd(ig,nlayermx+1)=0.
    13801272      enddo
    13811273
    13821274! calcul des tendances.
    13831275!-----------------------
    1384       do k=1,nlay
    1385          do ig=1,ngrid
     1276      do k=1,nlayermx
     1277         do ig=1,ngridmx
    13861278            pduadj(ig,k)=((detr_dv2(ig,k)+gamma(ig,k))*zua(ig,k)  &
    13871279     &               -(entr(ig,k)+gamma(ig,k))*ue(ig,k)  &
    … …  
    14121304
    14131305      modname='momentum'
    1414       call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
     1306      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
    14151307     &     masse,zu,pduadj,ztvd,fm_down,ztv,modname,lmax)
    14161308
    1417       call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
     1309      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
    14181310     &     masse,zv,pdvadj,ztvd,fm_down,ztv,modname,lmax)
    14191311
    … …  
    14241316!------------------------------------------------------------------
    14251317
    1426       do l=1,nlay
    1427          do ig=1,ngrid
     1318      do l=1,nlayermx
     1319         do ig=1,ngridmx
    14281320           pdtadj(ig,l)=(zdthladj(ig,l)+zdthladj_down(ig,l))*zpopsk(ig,l)
    14291321         enddo
    … …  
    14341326!------------------------------------------------------------------
    14351327
    1436       if (nq .ne. 0.) then
     1328      if (nqmx .ne. 0.) then
    14371329      modname='tracer'
    1438       DO iq=1,nq
    1439           call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
     1330      DO iq=1,nqmx
     1331          call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
    14401332          &      masse,pq(:,:,iq),pdqadj(:,:,iq),ztvd,fm_down,ztv,modname,lmax)
    14411333
    … …  
    14481340
    14491341      modname='tke'
    1450       call thermcell_dqupdown(ngrid,nlay,ptimestep,fm,entr,detr,  &
     1342      call thermcell_dqupdown(ngridmx,nlayermx,ptimestep,fm,entr,detr,  &
    14511343      &      masse,pq2,pdq2adj,ztvd,fm_down,ztv,modname,lmax)
    14521344
    … …  
    14641356
    14651357   
    1466       do l=1,nlay-1
    1467        do ig=1,ngrid
     1358      do l=1,nlayermx-1
     1359       do ig=1,ngridmx
    14681360        teta_th_int(ig,l)=0.5*(ztva(ig,l+1)+ztva(ig,l))
    14691361        teta_down_int(ig,l) = 0.5*(ztvd(ig,l+1)+ztvd(ig,l))
    … …  
    14711363       enddo
    14721364      enddo
    1473       do ig=1,ngrid
    1474        teta_th_int(ig,nlay)=teta_th_int(ig,nlay-1)
    1475        teta_env_int(ig,nlay)=teta_env_int(ig,nlay-1)
    1476        teta_down_int(ig,nlay)=teta_down_int(ig,nlay-1)
    1477       enddo
    1478       do l=1,nlay
    1479        do ig=1,ngrid
     1365      do ig=1,ngridmx
     1366       teta_th_int(ig,nlayermx)=teta_th_int(ig,nlayermx-1)
     1367       teta_env_int(ig,nlayermx)=teta_env_int(ig,nlayermx-1)
     1368       teta_down_int(ig,nlayermx)=teta_down_int(ig,nlayermx-1)
     1369      enddo
     1370      do l=1,nlayermx
     1371       do ig=1,ngridmx
    14801372        heatFlux(ig,l)=fm(ig,l)*(teta_th_int(ig,l)-teta_env_int(ig,l))/(rhobarz(ig,l))
    1481         buoyancyOut(ig,l)=RG*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
    1482         buoyancyEst(ig,l)=RG*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
     1373        buoyancyOut(ig,l)=g*(ztva(ig,l)-ztv(ig,l))/ztv(ig,l)
     1374        buoyancyEst(ig,l)=g*(ztva_est(ig,l)-ztv(ig,l))/ztv(ig,l)
    14831375        heatFlux_down(ig,l)=fm_down(ig,l)*(teta_down_int(ig,l)-teta_env_int(ig,l))/rhobarz(ig,l)
    14841376       enddo
Note: See TracChangeset for help on using the changeset viewer.