Ignore:
Timestamp:
Sep 19, 2012, 8:26:55 PM (12 years ago)
Author:
aslmd
Message:

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

  • nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this
Location:
trunk/LMDZ.GENERIC/libf
Files:
6 added
8 deleted
49 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/dyn3d/calfis.F

    r751 r787  
    7272#include "comgeom2.h"
    7373#include "control.h"
     74
     75#include "advtrac.h"
     76!! this is to get tnom (tracers name)
    7477
    7578c    Arguments :
    … …  
    422425
    423426      CALL physiq (ngridmx,llm,nq,
     427     .     tnom,
    424428     ,     debut,lafin,
    425429     ,     rday_ecri,heure,dtphys,

  • trunk/LMDZ.GENERIC/libf/dyn3d/ini_archive.F

    r711 r787  
    3434
    3535c=======================================================================
     36
     37      USE comsoil_h
    3638 
    3739      implicit none
    … …  
    4951#include "serre.h"
    5052#include "control.h"
    51 #include"comsoil.h"
    5253
    5354#include "netcdf.inc"

  • trunk/LMDZ.GENERIC/libf/dyn3d/lect_start_archive.F

    r588 r787  
    33     &     q,qsurf,surfith,nid)
    44
     5      USE surfdat_h
     6      USE comsoil_h
     7      USE tracer_h
     8
    59c=======================================================================
    610c
    … …  
    2226#include "dimensions.h"
    2327#include "dimphys.h"
    24 #include "surfdat.h"
    25 #include "comsoil.h"
    2628#include "planete.h"
    2729#include "paramet.h"
    … …  
    3638#include "lmdstd.h"
    3739#include "netcdf.inc"
    38 #include "tracer.h"
    3940#include"advtrac.h"
    4041c=======================================================================

  • trunk/LMDZ.GENERIC/libf/dyn3d/newstart.F

    r699 r787  
    1515c=======================================================================
    1616
     17      USE tracer_h
     18      USE comsoil_h
     19      USE surfdat_h
     20      USE comgeomfi_h
     21
    1722      implicit none
    1823
    1924#include "dimensions.h"
    2025#include "dimphys.h"
    21 #include "surfdat.h"
    22 #include "comsoil.h"
    2326#include "planete.h"
    2427#include "paramet.h"
    … …  
    3740#include "netcdf.inc"
    3841#include "advtrac.h"
    39 #include "tracer.h"
    4042c=======================================================================
    4143c   Declarations
    … …  
    243245      ! tnom(:) now contains tracer names
    244246! JL12 we will need the tracer names to read start in dyneta0
     247
     248! Initialize global tracer indexes (stored in tracer.h)
     249! ... this has to be done before phyetat0
     250      call initracer(ngridmx,nqmx,tnom)
     251
    245252
    246253c-----------------------------------------------------------------------
    … …  
    293300     .       ps,phis,time)
    294301
     302
     303        ! ALLOCATE ARRAYS IN comgeomfi_h (done in inifis usually)
     304        IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngridmx))
     305        IF (.not. ALLOCATED(long)) ALLOCATE(long(ngridmx))
     306        IF (.not. ALLOCATED(area)) ALLOCATE(area(ngridmx))
     307        cursor = 1 ! added by AS in dimphys. 1 for sequential runs.
     308
    295309        write(*,*) 'Reading file STARTFI'
    296310        fichnom = 'startfi.nc'
    297         CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,nqmx,
     311        CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqmx,
    298312     .        day_ini,time,
    299313     .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
    … …  
    352366c-----------------------------------------------------------------------
    353367      if (choix_1.eq.0) then
    354          call tabfi (nid,Lmodif,tab0,day_ini,lllm,p_rad,
     368         call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
    355369     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    356370      else if (choix_1.eq.1) then
    357371         Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
    358          call tabfi (nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
     372         call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
    359373     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    360374      endif
    … …  
    369383      kappa = 8.314*1000./(p_mugaz*p_cpp) ! added by RDW
    370384
    371 
    372385c=======================================================================
    373386c  INITIALISATIONS DIVERSES
    374387c=======================================================================
    375 ! Initialize global tracer indexes (stored in tracer.h)
    376       call initracer()
    377388! Load parameters from run.def file
    378389      CALL defrun_new( 99, .TRUE. )
    … …  
    14131424
    14141425
    1415       call physdem1("restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
     1426      call physdem1(ngridmx,"restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
    14161427     .              dtphys,float(day_ini),
    14171428     .              time,tsurf,tsoil,emis,q2,qsurf,
    14181429     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
    1419      .              cloudfrac,totalfrac,hice)
     1430     .              cloudfrac,totalfrac,hice,tnom)
    14201431
    14211432c=======================================================================

  • trunk/LMDZ.GENERIC/libf/dyn3d/start2archive.F

    r711 r787  
    1818c
    1919c=======================================================================
     20
     21      USE comsoil_h
    2022
    2123      implicit none
    … …  
    3436
    3537#include "dimphys.h"
    36 #include "comsoil.h"
    3738#include "planete.h"
    3839#include"advtrac.h"
    … …  
    162163      Lmodif=0
    163164
    164       CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,timefi,
     165      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,
     166     .      timefi,
    165167     .      tsurf,tsoil,emis,q2,qsurf,
    166168!       change FF 05/2011

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

    r741 r787  
    44       use radinc_h, only : L_TAUMAX,naerkind
    55       use aerosol_mod
     6       USE comgeomfi_h
     7       USE tracer_h
    68                 
    79       implicit none
    … …  
    2830!     pq                Aerosol mixing ratio
    2931!     reffrad(ngrid,nlayer,naerkind)         Aerosol effective radius
    30 !     QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
    31 !     QREFir3d(ngridmx,nlayermx,naerkind)  / at reference wavelengths
     32!     QREFvis3d(ngrid,nlayermx,naerkind) \ 3d extinction coefficients
     33!     QREFir3d(ngrid,nlayermx,naerkind)  / at reference wavelengths
    3234!
    3335!     Output
    … …  
    4244#include "callkeys.h"
    4345#include "comcstfi.h"
    44 #include "comgeomfi.h"
    45 #include "tracer.h"
    4646#include "comvert.h"
    4747
    … …  
    5353      REAL aerosol(ngrid,nlayer,naerkind)
    5454      REAL reffrad(ngrid,nlayer,naerkind)
    55       REAL QREFvis3d(ngridmx,nlayermx,naerkind)
    56       REAL QREFir3d(ngridmx,nlayermx,naerkind)
     55      REAL QREFvis3d(ngrid,nlayermx,naerkind)
     56      REAL QREFir3d(ngrid,nlayermx,naerkind)
    5757
    5858      REAL tau_col(ngrid)
    5959!      REAL tauref(ngrid), tau_col(ngrid)
    6060
    61       real cloudfrac(ngridmx,nlayermx)
     61      real cloudfrac(ngrid,nlayermx)
    6262      real aerosol0
    6363
    … …  
    7777      ! for fixed dust profiles
    7878      real topdust, expfactor, zp
    79       REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling
    80       REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling
     79      REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
     80      REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
    8181      ! BENJAMIN MODIFS
    8282      real CLFtot
    83       real totcloudfrac(ngridmx)
     83      real totcloudfrac(ngrid)
    8484      logical clearsky
    8585
    … …  
    8787      IF (firstcall) THEN
    8888
    89          ! are these tests of any real use ?
    90         IF(ngrid.NE.ngridmx) THEN
    91             PRINT*,'STOP in aeropacity!'
    92             PRINT*,'problem of dimensions:'
    93             PRINT*,'ngrid  =',ngrid
    94             PRINT*,'ngridmx  =',ngridmx
    95             STOP
    96         ENDIF
    97 
    98         if (nq.gt.nqmx) then
    99            write(*,*) 'STOP in aeropacity: (nq .gt. nqmx)!'
    100            write(*,*) 'nq=',nq,' nqmx=',nqmx
    101            stop
    102         endif
    103 
    10489        write(*,*) "Tracers found in aeropacity:"
    105         do iq=1,nqmx
     90        do iq=1,nq
    10691          tracername=noms(iq)
    10792          if (tracername.eq."co2_ice") then
    … …  
    150135           iaer=iaero_co2
    151136!       1. Initialization
    152             aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
     137            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
    153138!       2. Opacity calculation
    154139            if (noaero) then ! aerosol set to a very low value
    155              aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
     140             aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
    156141            elseif (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
    157                aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
    158                !aerosol(1:ngridmx,12,iaer)=4.0 ! single cloud layer option
     142               aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
     143               !aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
    159144            else
    160145               DO ig=1, ngrid
    … …  
    187172           iaer=iaero_h2o
    188173!       1. Initialization
    189             aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
     174            aerosol(1:ngrid,1:nlayermx,iaer)=0.0
    190175!       2. Opacity calculation
    191176            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
    192                aerosol(1:ngridmx,1:nlayermx,iaer) =1.e-9
     177               aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9
    193178
    194179               ! put cloud at cloudlvl
    … …  
    241226
    242227               if(CLFvarying)then
    243                   call totalcloudfrac(cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngridmx,1:nlayermx,iaer))
     228                  call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer))
    244229                  do ig=1, ngrid
    245230                     do l=1,nlayer-1 ! to stop the rad tran bug
    … …  
    268253          iaer=iaero_dust
    269254!         1. Initialization
    270           aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
     255          aerosol(1:ngrid,1:nlayermx,iaer)=0.0
    271256         
    272257          topdust=30.0 ! km  (used to be 10.0 km) LK
    … …  
    322307
    323308!       1. Initialization
    324          aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
     309         aerosol(1:ngrid,1:nlayermx,iaer)=0.0
    325310
    326311
    … …  
    386371      end do
    387372
    388       do ig=1, ngrid
     373      do ig=1,ngrid
    389374         do l=1,nlayer
    390375            do iaer = 1, naerkind
    … …  
    399384      end do
    400385
    401       do ig=1, ngrid
     386      do ig=1,ngrid
    402387         if(tau_col(ig).gt.1.e3)then
    403388            print*,'WARNING: tau_col=',tau_col(ig)

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

    r726 r787  
    6060      INTEGER :: grid_i,grid_j
    6161!     Intermediate variable
    62       REAL :: var_tmp,var3d_tmp(ngridmx,nlayermx)
     62      REAL :: var_tmp,var3d_tmp(ngrid,nlayermx)
    6363!     Bilinear interpolation factors
    6464      REAL :: kx,ky,k1,k2,k3,k4
    … …  
    161161      INTEGER :: ngrid,nlayer
    162162!     Aerosol effective radius used for radiative transfer (meter)
    163       REAL :: reffrad(ngridmx,nlayermx,naerkind)
     163      REAL :: reffrad(ngrid,nlayermx,naerkind)
    164164!     Aerosol effective variance used for radiative transfer (n.u.)
    165       REAL :: nueffrad(ngridmx,nlayermx,naerkind)
     165      REAL :: nueffrad(ngrid,nlayermx,naerkind)
    166166
    167167!     Outputs
    168168!     -------
    169169
    170       REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
    171       REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
    172       REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
    173 
    174       REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
    175       REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
    176       REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
    177 
    178       REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
    179       REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
    180 
    181       REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
    182       REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)
     170      REAL :: QVISsQREF3d(ngrid,nlayermx,L_NSPECTV,naerkind)
     171      REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
     172      REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
     173
     174      REAL :: QIRsQREF3d(ngrid,nlayermx,L_NSPECTI,naerkind)
     175      REAL :: omegaIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
     176      REAL :: gIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
     177
     178      REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
     179      REAL :: QREFir3d(ngrid,nlayermx,naerkind)
     180
     181      REAL :: omegaREFvis3d(ngrid,nlayermx,naerkind)
     182      REAL :: omegaREFir3d(ngrid,nlayermx,naerkind)
    183183
    184184      DO iaer = 1, naerkind ! Loop on aerosol kind
    … …  
    380380
    381381      DO lg = 1,nlayer
    382         DO ig = 1,ngrid
     382        DO ig = 1, ngrid
    383383!         2.1 Effective radius index and kx calculation
    384384          var_tmp=reffrad(ig,lg,iaer)/refftabmin
    … …  
    725725          ENDIF                  ! --------------------------------
    726726        ENDDO !nlayermx
    727       ENDDO !ngridmx
     727      ENDDO !ngrid
    728728
    729729!==================================================================

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

    r305 r787  
    2222  real Play(1:nlayermx),Plev(1:nlayermx+1)
    2323  real Qsurf,Q(1:nlayermx)
    24   real muvar(ngridmx,nlayermx+1)
     24  real muvar(1,nlayermx+1)
    2525
    2626  ! internal

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

    r728 r787  
    55          fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,        &
    66          OLR_nu,OSR_nu,                                       &
    7           reffrad,tau_col,cloudfrac,totcloudfrac,              &
     7          reffrad,nueffrad,tau_col,cloudfrac,totcloudfrac,     &
    88          clearsky,firstcall,lastcall)
    99
    … …  
    1616      use radii_mod, only : su_aer_radii,co2_reffrad,h2o_reffrad,dust_reffrad,h2so4_reffrad
    1717      use aerosol_mod
     18
     19      USE tracer_h
    1820
    1921      implicit none
    … …  
    3739#include "comcstfi.h"
    3840#include "callkeys.h"
    39 #include "tracer.h"
    4041
    4142!-----------------------------------------------------------------------
    … …  
    6061!     Globally varying aerosol optical properties on GCM grid
    6162!     Not needed everywhere so not in radcommon_h
    62       REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
    63       REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
    64       REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind)
    65 
    66       REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
    67       REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
    68       REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind)
    69 
    70       REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
    71       REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
    72 
    73 !      REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
    74 !      REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these...
    75 
    76       REAL reffrad(ngrid,nlayer,naerkind)
    77       REAL, SAVE :: nueffrad(ngridmx,nlayermx,naerkind)
    78 
    79 !     OUTPUT
    80       REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW
    81       REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW
    82       REAL fluxsurf_lw(ngridmx)   ! incident LW flux to surf (W/m2)
    83       REAL fluxtop_lw(ngridmx)    ! outgoing LW flux to space (W/m2)
    84       REAL fluxsurf_sw(ngridmx)   ! incident SW flux to surf (W/m2)
    85       REAL fluxabs_sw(ngridmx)    ! SW flux absorbed by planet (W/m2)
    86       REAL fluxtop_dn(ngridmx)    ! incident top of atmosphere SW flux (W/m2)
     63      REAL :: QVISsQREF3d(ngrid,nlayermx,L_NSPECTV,naerkind)
     64      REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
     65      REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
     66
     67      REAL :: QIRsQREF3d(ngrid,nlayermx,L_NSPECTI,naerkind)
     68      REAL :: omegaIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
     69      REAL :: gIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
     70
     71!      REAL :: omegaREFvis3d(ngrid,nlayermx,naerkind)
     72!      REAL :: omegaREFir3d(ngrid,nlayermx,naerkind) ! not sure of the point of these...
     73
     74      !!! this is a local instance of a variable saved in physiq.F and being an argument
     75      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
     76      REAL, INTENT(INOUT) :: nueffrad(ngrid,nlayermx,naerkind)
     77
     78!!     OUTPUT
     79      REAL dtsw(ngrid,nlayermx) ! heating rate (K/s) due to SW
     80      REAL dtlw(ngrid,nlayermx) ! heating rate (K/s) due to LW
     81      REAL fluxsurf_lw(ngrid)   ! incident LW flux to surf (W/m2)
     82      REAL fluxtop_lw(ngrid)    ! outgoing LW flux to space (W/m2)
     83      REAL fluxsurf_sw(ngrid)   ! incident SW flux to surf (W/m2)
     84      REAL fluxabs_sw(ngrid)    ! SW flux absorbed by planet (W/m2)
     85      REAL fluxtop_dn(ngrid)    ! incident top of atmosphere SW flux (W/m2)
    8786      REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
    8887      REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
     88
    8989!-----------------------------------------------------------------------
    9090!     Declaration of the variables required by correlated-k subroutines
    … …  
    126126
    127127      real*8 qvar(L_LEVELS)          ! mixing ratio of variable component (mol/mol)
    128       REAL pq(ngridmx,nlayer,nq)
    129       REAL qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
     128      REAL pq(ngrid,nlayer,nq)
     129      REAL qsurf(ngrid,nq)       ! tracer on surface (e.g. kg.m-2)
    130130      integer nq
    131131
    … …  
    140140      save qxvaer, qsvaer, gvaer
    141141      save qxiaer, qsiaer, giaer
    142       save QREFvis3d, QREFir3d
     142
     143      !REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
     144      !REAL :: QREFir3d(ngrid,nlayermx,naerkind)
     145      !save QREFvis3d, QREFir3d
     146      real, dimension(:,:,:), save, allocatable :: QREFvis3d
     147      real, dimension(:,:,:), save, allocatable :: QREFir3d
    143148
    144149      REAL tau_col(ngrid) ! diagnostic from aeropacity
    … …  
    165170
    166171!     for H2O cloud fraction in aeropacity
    167       real cloudfrac(ngridmx,nlayermx)
    168       real totcloudfrac(ngridmx)
     172      real cloudfrac(ngrid,nlayermx)
     173      real totcloudfrac(ngrid)
    169174      logical clearsky
    170175
    171176      ! for weird cloud test
    172       real pqtest(ngridmx,nlayer,nq)
     177      real pqtest(ngrid,nlayer,nq)
    173178
    174179      real maxrad, minrad
    … …  
    178183
    179184!     included by RW for runaway greenhouse 1D study
    180       real muvar(ngridmx,nlayermx+1)
     185      real muvar(ngrid,nlayermx+1)
    181186      real vtmp(nlayermx)
    182187      REAL*8 muvarrad(L_LEVELS)
    … …  
    195200
    196201      if(firstcall) then
     202
     203         ALLOCATE(QREFvis3d(ngrid,nlayermx,naerkind))
     204         ALLOCATE(QREFir3d(ngrid,nlayermx,naerkind))
    197205
    198206         call system('rm -f surf_vals_long.out')
    … …  
    204212            call abort
    205213         endif
    206          call su_aer_radii(reffrad,nueffrad)
     214         call su_aer_radii(ngrid,reffrad,nueffrad)
    207215         
    208216         
    … …  
    233241         OSR_nu(:,:) = 0.
    234242
    235          if (ngridmx.eq.1) then
     243         if (ngrid.eq.1) then
    236244           PRINT*, 'Simulate global averaged conditions ?'
    237245           global1d = .false. ! default value
    … …  
    266274
    267275         if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! treat condensed co2 particles.
    268             call co2_reffrad(pq,reffrad)
    269             print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
    270             print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
     276            call co2_reffrad(ngrid,nq,pq,reffrad)
     277            print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
     278            print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
    271279         end if
    272280         if ((iaer.eq.iaero_h2o).and.water) then ! treat condensed water particles. to be generalized for other aerosols
    273             call h2o_reffrad(pq,pt,reffrad,nueffrad)
    274             print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
    275             print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngridmx,1:nlayermx,iaer))/1.e-6,' um'
     281            call h2o_reffrad(ngrid,nq,pq,pt,reffrad,nueffrad)
     282            print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
     283            print*,'Min. H2O cloud particle size = ',minval(reffrad(1:ngrid,1:nlayermx,iaer))/1.e-6,' um'
    276284         endif
    277285         if(iaer.eq.iaero_dust)then
    278             call dust_reffrad(reffrad)
     286            call dust_reffrad(ngrid,reffrad)
    279287            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
    280288         endif
    281289         if(iaer.eq.iaero_h2so4)then
    282             call h2so4_reffrad(reffrad)
     290            call h2so4_reffrad(ngrid,reffrad)
    283291            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
    284292         endif
    … …  
    304312!-----------------------------------------------------------------------
    305313!     Starting Big Loop over every GCM column
    306       do ig=1,ngridmx
     314      do ig=1,ngrid
    307315
    308316!=======================================================================
    … …  
    452460      endif
    453461
    454       if ((ngridmx.eq.1).and.(global1d)) then       ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
     462      if ((ngrid.eq.1).and.(global1d)) then       ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight
    455463         acosz = cos(pi*szangle/180.0)
    456464         print*,'acosz=',acosz,', szangle=',szangle
    … …  
    550558         !i_var=igcm_h2o_vap
    551559         i_var=1
    552          if(nqmx.gt.1)then
     560         if(nq.gt.1)then
    553561            print*,'Need 1 tracer only to run kcm1d.e'
    554562            stop
    … …  
    653661            fluxtoplanet=0.
    654662
    655             if((ngridmx.eq.1).and.(global1d))then
     663            if((ngrid.eq.1).and.(global1d))then
    656664               do nw=1,L_NSPECTV
    657665                  stel_fract(nw)= stel(nw) * 0.25 / acosz

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

    r486 r787  
    22     $                pplev,zlev, pt,
    33     &                pq, pdqfi, pdqsed,pdqs_sed,nq,rfall)
     4
     5      USE tracer_h
     6
    47      IMPLICIT NONE
    58
    … …  
    2528#include "dimphys.h"
    2629#include "comcstfi.h"
    27 #include "tracer.h"
    2830#include "callkeys.h"
    29 
    30 #include "fisice.h"
    3131
    3232c   arguments:
    … …  
    5353      INTEGER l,ig, iq
    5454
    55       real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers
    56       real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2)
    57       real epaisseur (ngridmx,nlayermx) ! Layer thickness (m)
    58       real wq(ngridmx,nlayermx+1) ! displaced tracer mass (kg.m-2)
    59 c      real dens(ngridmx,nlayermx) ! Mean density of the ice part. accounting for dust core
    60       real rfall(ngridmx,nlayermx)
     55      real zqi(ngrid,nlayermx,nq) ! to locally store tracers
     56      real masse (ngrid,nlayermx) ! Layer mass (kg.m-2)
     57      real epaisseur (ngrid,nlayermx) ! Layer thickness (m)
     58      real wq(ngrid,nlayermx+1) ! displaced tracer mass (kg.m-2)
     59c      real dens(ngrid,nlayermx) ! Mean density of the ice part. accounting for dust core
     60      real rfall(ngrid,nlayermx)
    6161
    6262
    … …  
    6969
    7070      IF (firstcall) THEN
    71          IF(ngrid.NE.ngridmx) THEN
    72             PRINT*,'STOP dans callsedim'
    73             PRINT*,'probleme de dimensions :'
    74             PRINT*,'ngrid  =',ngrid
    75             PRINT*,'ngridmx  =',ngridmx
    76             STOP
    77          ENDIF
    78      
    7971        firstcall=.false.
    8072      ENDIF ! of IF (firstcall)
    … …  
    9688!  of general tracers.]
    9789 
    98       zqi(1:ngrid,1:nlay,1:nqmx) = 0.
     90      zqi(1:ngrid,1:nlay,1:nq) = 0.
    9991      do iq=1,nq
    10092       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
    … …  
    10496
    10597          do l=1,nlay
    106             do ig=1,ngrid
     98            do ig=1, ngrid
    10799              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
    108100            enddo
    … …  
    128120!======================================================================
    129121
    130           do ig=1,ngrid 
     122          do ig=1,ngrid
    131123!     Ehouarn: with new way of tracking tracers by name, this is simply
    132124            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep

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

    r728 r787  
    1010      use radii_mod, only : co2_reffrad
    1111      use aerosol_mod, only : iaero_co2
     12      USE surfdat_h
     13      USE comgeomfi_h
     14      USE tracer_h
     15
    1216
    1317      implicit none
    … …  
    5559#include "dimphys.h"
    5660#include "comcstfi.h"
    57 #include "surfdat.h"
    58 #include "comgeomfi.h"
    5961#include "comvert.h"
    6062#include "callkeys.h"
    61 #include "tracer.h"
    6263
    6364!-----------------------------------------------------------------------
    … …  
    8081      REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer)
    8182      REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer)
    82       REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
     83      REAL pq(ngrid,nlayer,nq),pdq(ngrid,nlayer,nq)
    8384      REAL pdqc(ngrid,nlayer,nq), pdqsc(ngrid,nq)
    8485
    … …  
    9293      INTEGER l,ig,icap,ilay,it,iq
    9394
    94       REAL*8 zt(ngridmx,nlayermx)
    95       REAL zq(ngridmx,nlayermx,nqmx)
     95      REAL*8 zt(ngrid,nlayermx)
     96      REAL zq(ngrid,nlayermx,nq)
    9697      REAL zcpi
    97       REAL ztcond (ngridmx,nlayermx)
    98       REAL ztnuc (ngridmx,nlayermx)
    99       REAL ztcondsol(ngridmx)
    100       REAL zdiceco2(ngridmx)
    101       REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx)
    102       REAL zfallice(ngridmx), Mfallice(ngridmx)
     98      REAL ztcond (ngrid,nlayermx)
     99      REAL ztnuc (ngrid,nlayermx)
     100      REAL ztcondsol(ngrid)
     101      REAL zdiceco2(ngrid)
     102      REAL zcondicea(ngrid,nlayermx), zcondices(ngrid)
     103      REAL zfallice(ngrid), Mfallice(ngrid)
    103104      REAL zmflux(nlayermx+1)
    104105      REAL zu(nlayermx),zv(nlayermx)
    105       REAL ztsrf(ngridmx)
     106      REAL ztsrf(ngrid)
    106107      REAL ztc(nlayermx), ztm(nlayermx+1)
    107108      REAL zum(nlayermx+1) , zvm(nlayermx+1)
    108       LOGICAL condsub(ngridmx)
     109      LOGICAL condsub(ngrid)
    109110      REAL subptimestep
    110111      Integer Ntime
    111       real masse (ngridmx,nlayermx), w(ngridmx,nlayermx,nqmx)
    112       real wq(ngridmx,nlayermx+1)
     112      real masse (ngrid,nlayermx), w(ngrid,nlayermx,nq)
     113      real wq(ngrid,nlayermx+1)
    113114      real vstokes,reff
    114115
    115116!     Special diagnostic variables
    116       real tconda1(ngridmx,nlayermx)
    117       real tconda2(ngridmx,nlayermx)
    118       real zdtsig (ngridmx,nlayermx)
    119       real zdt (ngridmx,nlayermx)
     117      real tconda1(ngrid,nlayermx)
     118      real tconda2(ngrid,nlayermx)
     119      real zdtsig (ngrid,nlayermx)
     120      real zdt (ngrid,nlayermx)
    120121
    121122!-----------------------------------------------------------------------
    122123!     Saved local variables
    123124
    124       REAL emisref(ngridmx)
    125125      REAL latcond
    126126      REAL ccond
    127127      REAL cpice
    128       SAVE emisref, cpice
     128      REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref
     129      SAVE cpice
    129130      SAVE latcond,ccond
    130131
    … …  
    155156      call zerophys(ngrid*nlayer*nq,pdqc)
    156157      call zerophys(ngrid*nlayer*nq,pdtc)
    157       call zerophys(ngridmx*nlayermx*nqmx,zq)
    158       call zerophys(ngridmx*nlayermx,zt)
     158      call zerophys(ngrid*nlayermx*nq,zq)
     159      call zerophys(ngrid*nlayermx,zt)
    159160
    160161!     Initialisation
    161162      IF (firstcall) THEN
    162163
     164         ALLOCATE(emisref(ngrid)) !! this should be deallocated in lastcall...
     165
    163166         ! find CO2 ice tracer
    164          do iq=1,nqmx
     167         do iq=1,nq
    165168            tracername=noms(iq)
    166169            if (tracername.eq."co2_ice") then
    … …  
    287290!     Calculate the mass of each atmospheric layer (kg.m-2)
    288291      do  ilay=1,nlayer
    289          do ig=1, ngrid
     292         DO ig=1,ngrid
    290293            masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g
    291294         end do
    … …  
    312315           
    313316!     sedimentation computed from radius computed from q in module radii_mod
    314          call co2_reffrad(zq,reffrad)
     317         call co2_reffrad(ngrid,nq,zq,reffrad)
    315318         
    316319         do  ilay=1,nlayer
    317             do ig=1, ngrid
     320            DO ig=1,ngrid
    318321
    319322               reff = reffrad(ig,ilay,iaero_co2)
    … …  
    337340!     Progressively accumulating the flux to the ground
    338341!     Mfallice is the total amount of ice fallen to the ground
    339          do ig=1,ngrid
     342         DO ig=1,ngrid
    340343            Mfallice(ig) =  Mfallice(ig) + wq(ig,i_co2ice)
    341344         end do
    … …  
    463466!     Surface albedo and emissivity of the ground below the snow (emisref)
    464467!     --------------------------------------------------------------------
    465       do ig =1,ngrid
     468      DO ig=1,ngrid
    466469         IF(ig.GT.ngrid/2+1) THEN
    467470            icap=2

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

    r253 r787  
    55     &                   pduadj,pdvadj,pdhadj,
    66     &                   pdqadj)
     7
     8      USE tracer_h
    79
    810      implicit none
    … …  
    3133#include "comcstfi.h"
    3234#include "callkeys.h"
    33 #include "tracer.h"
    3435
    3536
    … …  
    5455
    5556      INTEGER ig,i,l,l1,l2,jj
    56       INTEGER jcnt, jadrs(ngridmx)
     57      INTEGER jcnt, jadrs(ngrid)
    5758
    5859      REAL sig(nlayermx+1),sdsig(nlayermx),dsig(nlayermx)
    59       REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
    60       REAL zh(ngridmx,nlayermx)
    61       REAL zu2(ngridmx,nlayermx),zv2(ngridmx,nlayermx)
    62       REAL zh2(ngridmx,nlayermx), zhc(ngridmx,nlayermx)
     60      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
     61      REAL zh(ngrid,nlayermx)
     62      REAL zu2(ngrid,nlayermx),zv2(ngrid,nlayermx)
     63      REAL zh2(ngrid,nlayermx), zhc(ngrid,nlayermx)
    6364      REAL zhm,zsm,zdsm,zum,zvm,zalpha,zhmc
    6465
    … …  
    6667      INTEGER iq,ico2
    6768      save ico2
    68       REAL zq(ngridmx,nlayermx,nqmx), zq2(ngridmx,nlayermx,nqmx)
    69       REAL zqm(nqmx),zqco2m
     69      REAL zq(ngrid,nlayermx,nq), zq2(ngrid,nlayermx,nq)
     70      REAL zqm(nq),zqco2m
    7071      real m_co2, m_noco2, A , B
    7172      save A, B
    … …  
    7374      real mtot1, mtot2 , mm1, mm2
    7475       integer l1ref, l2ref
    75       LOGICAL vtest(ngridmx),down,firstcall
     76      LOGICAL vtest(ngrid),down,firstcall
    7677      save firstcall
    7778      data firstcall/.true./
    … …  
    8788
    8889      IF (firstcall) THEN
    89         IF(ngrid.NE.ngridmx) THEN
    90            PRINT*
    91            PRINT*,'STOP in convadj'
    92            PRINT*,'ngrid    =',ngrid
    93            PRINT*,'ngridmx  =',ngridmx
    94         ENDIF
    9590        ico2=0
    9691        if (tracer) then
    9792!     Prepare Special treatment if one of the tracers is CO2 gas
    98            do iq=1,nqmx
     93           do iq=1,nq
    9994             if (noms(iq).eq."co2") then
    10095                print*,'dont go there'
    … …  
    401396        do iq=1, nq
    402397          do  l=1,nlay
    403             do ig=1, ngrid
     398            DO ig=1,ngrid
    404399              pdqadj(ig,l,iq)=(zq2(ig,l,iq)-zq(ig,l,iq))/ptimestep
    405400            end do

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

    r588 r787  
    8585      INTEGER    i,j,k
    8686
    87       INTEGER klatdat,ngridmxgdat
    88       PARAMETER (klatdat=180,ngridmxgdat=360)
     87      INTEGER klatdat,ngridmixgdat
     88      PARAMETER (klatdat=180,ngridmixgdat=360)
    8989
    9090c    on passe une grille en rlonu rlatv et im+1 jm a interp_horiz)
    … …  
    245245      do j=1,jmdp1
    246246          do i=1,imd
    247               zdataS(i+imdp1*(j-1)) = zdata(i+ngridmxgdat*(j-1))
     247              zdataS(i+imdp1*(j-1)) = zdata(i+ngridmixgdat*(j-1))
    248248          enddo
    249           zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmxgdat*(j-1))
     249          zdataS(imdp1+imdp1*(j-1)) = zdata(1+ngridmixgdat*(j-1))
    250250      enddo
    251251

  • trunk/LMDZ.GENERIC/libf/phystd/dimphys.h

    r786 r787  
    1212      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 cm, depth = 104.8 m => earth case
    1313!-----------------------------------------------------------------------
     14
     15      common/dimphys/cursor
     16      integer cursor  !! this is the position of point 1 on the globe
     17                      !! it is 1 obviously for sequential runs
     18                      !! but might differ for paraller runs
     19                      !! it is used in phyetat0 mostly
     20

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

    r650 r787  
    1       subroutine evap(dtime,pt, pq, pdq, pdt, 
     1      subroutine evap(ngrid,nq,dtime,pt, pq, pdq, pdt, 
    22     $     dqevap,dtevap, qevap, tevap)
    33       
    44      use watercommon_h
     5      USE tracer_h
    56
    67      implicit none
    … …  
    89#include "dimensions.h"
    910#include "dimphys.h"
    10 #include "tracer.h"
    1111#include "comcstfi.h"
    1212
    … …  
    2424!==================================================================
    2525
     26      INTEGER ngrid,nq
     27
    2628! Arguments:
    27       REAL pt(ngridmx,nlayermx)
    28       REAL pq(ngridmx,nlayermx,nqmx)
    29       REAL pdt(ngridmx,nlayermx)
    30       REAL pdq(ngridmx,nlayermx,nqmx)
    31       REAL dqevap(ngridmx,nlayermx)
    32       REAL dtevap(ngridmx,nlayermx)
    33       REAL qevap(ngridmx,nlayermx,nqmx)
     29      REAL pt(ngrid,nlayermx)
     30      REAL pq(ngrid,nlayermx,nq)
     31      REAL pdt(ngrid,nlayermx)
     32      REAL pdq(ngrid,nlayermx,nq)
     33      REAL dqevap(ngrid,nlayermx)
     34      REAL dtevap(ngrid,nlayermx)
     35      REAL qevap(ngrid,nlayermx,nq)
    3436      REAL dtime
    3537 
    3638! Local:
    37       REAL tevap(ngridmx,nlayermx)
     39      REAL tevap(ngrid,nlayermx)
    3840      REAL zlvdcp
    3941      REAL zlsdcp
    … …  
    4648
    4749      DO l=1,nlayermx
    48         DO ig=1,ngridmx
     50        DO ig=1,ngrid
    4951         qevap(ig,l,igcm_h2o_vap)=pq(ig,l,igcm_h2o_vap)
    5052     s                            +pdq(ig,l,igcm_h2o_vap)*dtime
    … …  
    5658
    5759      DO l = 1, nlayermx 
    58         DO ig = 1, ngridmx
     60        DO ig = 1, ngrid
    5961         zlvdcp=RLVTT/RCPD!/(1.0+RVTMP2*qevap(ig,l,igcm_h2o_vap))
    6062         zlsdcp=RLSTT/RCPD!/(1.0+RVTMP2*qevap(ig,l,igcm_h2o_vap))

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

    r253 r787  
    1       subroutine forceWCfn(pplev,pt,dq,dqs)
     1      subroutine forceWCfn(ngrid,nq,pplev,pt,dq,dqs)
     2
     3      USE tracer_h
    24
    35      implicit none
    … …  
    1921#include "dimphys.h"
    2022#include "comcstfi.h"
    21 #include "tracer.h"
     23
     24      INTEGER ngrid,nq
    2225
    2326      real masse, Wtot, Wdiff
    2427
    25       real pplev(ngridmx,nlayermx+1)
    26       real pt(ngridmx)
     28      real pplev(ngrid,nlayermx+1)
     29      real pt(ngrid)
    2730
    28       real dqs(ngridmx,nqmx)
    29       real dq(ngridmx,nlayermx,nqmx)
     31      real dqs(ngrid,nq)
     32      real dq(ngrid,nlayermx,nq)
    3033
    3134      integer iq, ig, ilay
    3235
    33       do iq=1,nqmx
    34         do ig=1,ngridmx
     36      do iq=1,nq
     37        do ig=1,ngrid
    3538           Wtot = 0.0
    3639           do ilay=1,nlayermx

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

    r773 r787  
    1 subroutine hydrol(ptimestep,rnat,tsurf,          &
     1subroutine hydrol(ngrid,nq,ptimestep,rnat,tsurf,          &
    22     qsurf,dqsurf,dqs_hyd,pcapcal,                &
    33     albedo0,albedo,mu0,pdtsurf,pdtsurf_hyd,hice)
    44
    55  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
     6  USE surfdat_h
     7  use comdiurn_h
     8  USE comgeomfi_h
     9  USE tracer_h
    610
    711  implicit none
    … …  
    3640#include "comcstfi.h"
    3741#include "callkeys.h"
    38 #include "tracer.h"
    39 #include "fisice.h"
    40 #include "comgeomfi.h"
    41 #include "comdiurn.h"
    42 #include "surfdat.h"
     42
     43      integer ngrid,nq
    4344
    4445!     Inputs
    … …  
    6061!     Arguments
    6162!     ---------
    62       integer rnat(ngridmx) ! I changed this to integer (RW)
    63       real runoff(ngridmx)
    64       save runoff
     63      integer rnat(ngrid) ! I changed this to integer (RW)
     64      real,dimension(:),allocatable,save :: runoff
    6565      real totalrunoff, tsea, oceanarea
    6666      save oceanarea
    6767
    6868      real ptimestep
    69       real mu0(ngridmx)
    70       real qsurf(ngridmx,nqmx), tsurf(ngridmx)
    71       real dqsurf(ngridmx,nqmx), pdtsurf(ngridmx)
    72       real hice(ngridmx)
    73       real albedo0(ngridmx), albedo(ngridmx)
     69      real mu0(ngrid)
     70      real qsurf(ngrid,nq), tsurf(ngrid)
     71      real dqsurf(ngrid,nq), pdtsurf(ngrid)
     72      real hice(ngrid)
     73      real albedo0(ngrid), albedo(ngrid)
    7474
    7575      real oceanarea2
    … …  
    7777!     Output
    7878!     ------
    79       real dqs_hyd(ngridmx,nqmx)
    80       real pdtsurf_hyd(ngridmx)
     79      real dqs_hyd(ngrid,nq)
     80      real pdtsurf_hyd(ngrid)
    8181
    8282!     Local
    … …  
    8585      integer ig,iq, icap ! wld like to remove icap
    8686      real fsnoi, subli, fauxo
    87       real twater(ngridmx)
    88       real pcapcal(ngridmx)
    89       real hicebis(ngridmx)
    90       real zqsurf(ngridmx,nqmx)
    91       real ztsurf(ngridmx)
     87      real twater(ngrid)
     88      real pcapcal(ngrid)
     89      real hicebis(ngrid)
     90      real zqsurf(ngrid,nq)
     91      real ztsurf(ngrid)
    9292
    9393      integer ivap, iliq, iice
    … …  
    104104
    105105      if(firstcall)then
     106
     107         ALLOCATE(runoff(ngrid))
    106108
    107109         ivap=igcm_h2o_vap
    … …  
    126128!     Total ocean surface area
    127129         oceanarea=0.
    128          do ig=1,ngridmx
     130         do ig=1,ngrid
    129131            if(rnat(ig).eq.0)then
    130132               oceanarea=oceanarea+area(ig)
    … …  
    150152!     ------------------------------------------
    151153
    152       do ig=1,ngridmx
     154      do ig=1,ngrid
    153155         ztsurf(ig) = tsurf(ig) + ptimestep*pdtsurf(ig)
    154156         pdtsurf_hyd(ig)=0.0
    155          do iq=1,nqmx
     157         do iq=1,nq
    156158            zqsurf(ig,iq) = qsurf(ig,iq) + ptimestep*dqsurf(ig,iq)
    157159         enddo   
    158160      enddo
    159161 
    160       do ig=1,ngridmx
    161          do iq=1,nqmx
     162      do ig=1,ngrid
     163         do iq=1,nq
    162164            dqs_hyd(ig,iq) = 0.0
    163165         enddo
    164166      enddo
    165167
    166       do ig = 1, ngridmx
     168      do ig = 1, ngrid
    167169
    168170!     Ocean
    … …  
    257259
    258260               runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
    259                if(ngridmx.gt.1)then ! runoff only exists in 3D
     261               if(ngrid.gt.1)then ! runoff only exists in 3D
    260262                  if(runoff(ig).ne.0.0)then
    261263                     zqsurf(ig,iliq) = mx_eau_sol
    … …  
    283285         endif
    284286
    285       end do ! ig=1,ngridmx
     287      end do ! ig=1,ngrid
    286288
    287289!     perform crude bulk averaging of temperature in ocean
    … …  
    290292
    291293         oceanarea2=0.       
    292          DO ig=1,ngridmx
     294         DO ig=1,ngrid
    293295            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then
    294296               oceanarea2=oceanarea2+area(ig)*pcapcal(ig)
    … …  
    297299       
    298300         tsea=0.
    299          DO ig=1,ngridmx
     301         DO ig=1,ngrid
    300302            if((rnat(ig).eq.0).and.(hice(ig).eq.0.))then       
    301303               tsea=tsea+ztsurf(ig)*area(ig)*pcapcal(ig)/oceanarea2
    … …  
    303305         END DO
    304306
    305          DO ig=1,ngridmx
     307         DO ig=1,ngrid
    306308            if((rnat(ig).eq.0).and.(hice(ig).eq.0))then
    307309               pdtsurf_hyd(ig) = pdtsurf_hyd(ig) + (tsea-ztsurf(ig))/oceantime
    … …  
    318320
    319321         totalrunoff=0.
    320          do ig=1,ngridmx
     322         do ig=1,ngrid
    321323            if (rnat(ig).eq.1) then
    322324               totalrunoff = totalrunoff + area(ig)*runoff(ig)
    … …  
    324326         enddo
    325327         
    326          do ig=1,ngridmx
     328         do ig=1,ngrid
    327329            if (rnat(ig).eq.0) then
    328330               zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
    … …  
    338340
    339341         icap=1
    340          do ig=1,ngridmx
     342         do ig=1,ngrid
    341343            if (qsurf(ig,igcm_co2_ice).gt.0) then
    342344               albedo(ig) = albedice(icap)
    … …  
    347349
    348350
    349       do ig=1,ngridmx
     351      do ig=1,ngrid
    350352         dqs_hyd(ig,iliq)=(zqsurf(ig,iliq) - qsurf(ig,iliq))/ptimestep
    351353         dqs_hyd(ig,iice)=(zqsurf(ig,iice) - qsurf(ig,iice))/ptimestep
    352354      enddo
    353355
    354       call writediagfi(ngridmx,'runoff','Runoff amount',' ',2,runoff)
     356      call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)
    355357
    356358      return

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

    r728 r787  
    66      use radinc_h, only : naerkind
    77      use datafile_mod, only: datadir
     8      use comdiurn_h
     9
     10      !! to be conservative, but why this include here?
     11      use surfdat_h
     12      use comsaison_h
     13
     14      USE comgeomfi_h
    815
    916!=======================================================================
    … …  
    5057#include "planete.h"
    5158#include "comcstfi.h"
    52 #include "comsaison.h"
    53 #include "comdiurn.h"
    54 #include "comgeomfi.h"
    5559#include "callkeys.h"
    56 #include "surfdat.h"
    5760
    5861
    … …  
    6164 
    6265      INTEGER ngrid,nlayer
    63       REAL plat(ngrid),plon(ngrid),parea(ngridmx)
     66      REAL plat(ngrid),plon(ngrid),parea(ngrid)
    6467      integer day_ini
    6568      INTEGER ig,ierr
    … …  
    8790      avocado = 6.02214179e23   ! added by RW
    8891
     92      cursor = 1 ! added by AS in dimphys. 1 for sequential runs.
     93
    8994! --------------------------------------------------------
    9095!     The usual Tests
    … …  
    95100         PRINT*,'nlayer     = ',nlayer
    96101         PRINT*,'nlayermx   = ',nlayermx
    97          STOP
    98       ENDIF
    99 
    100       IF (ngrid.NE.ngridmx) THEN
    101          PRINT*,'STOP in inifis'
    102          PRINT*,'Probleme de dimensions :'
    103          PRINT*,'ngrid     = ',ngrid
    104          PRINT*,'ngridmx   = ',ngridmx
    105102         STOP
    106103      ENDIF
    … …  
    299296! Test of incompatibility:
    300297! if kastprof used, we must be in 1D
    301          if (kastprof.and.(ngridmx.gt.1)) then
     298         if (kastprof.and.(ngrid.gt.1)) then
    302299           print*,'kastprof can only be used in 1D!'
    303300           call abort
    … …  
    350347! Test of incompatibility:
    351348! if testradtimes used, we must be in 1D
    352          if (testradtimes.and.(ngridmx.gt.1)) then
     349         if (testradtimes.and.(ngrid.gt.1)) then
    353350           print*,'testradtimes can only be used in 1D!'
    354351           call abort
    … …  
    595592!     ------------------------
    596593
     594      ! ALLOCATE ARRAYS IN comgeomfi_h
     595      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid))
     596      IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid))
     597      IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid))
     598
    597599      CALL SCOPY(ngrid,plon,1,long,1)
    598600      CALL SCOPY(ngrid,plat,1,lati,1)
    599601      CALL SCOPY(ngrid,parea,1,area,1)
    600       totarea=SSUM(ngridmx,area,1)
     602      totarea=SSUM(ngrid,area,1)
     603
     604      !! those are defined in comdiurn_h.F90
     605      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
     606      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
     607      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
     608      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
    601609
    602610      DO ig=1,ngrid

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

    r253 r787  
    1       SUBROUTINE initracer()
    2 
    3 
     1      SUBROUTINE initracer(ngrid,nq,nametrac)
     2
     3      use surfdat_h
     4      USE comgeomfi_h
     5      USE tracer_h
    46      IMPLICIT NONE
    57c=======================================================================
    … …  
    1315c   Test of dimension :
    1416c   Initialize COMMON tracer in tracer.h, using tracer names provided
    15 c   by the dynamics in "advtrac.h"
     17c   by the argument nametrac
    1618c
    1719c   author: F.Forget
    … …  
    2527#include "comcstfi.h"
    2628#include "callkeys.h"
    27 #include "tracer.h"
    28 #include "advtrac.h"
    29 #include "comgeomfi.h"
    30 #include "watercap.h"
    31 #include "fisice.h"
    32 
    33 !      real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g.  kg.m-2)
    34 !      real co2ice(ngridmx)           ! co2 ice mass on surface (e.g.  kg.m-2)
     29
     30      integer :: ngrid,nq
     31
     32!      real qsurf(ngrid,nq)       ! tracer on surface (e.g.  kg.m-2)
     33!      real co2ice(ngrid)           ! co2 ice mass on surface (e.g.  kg.m-2)
    3534      character(len=20) :: txt ! to store some text
    3635      integer iq,ig,count
    … …  
    3837!      logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
    3938
     39      character*20 nametrac(nq)   ! name of the tracer from dynamics
     40
    4041
    4142c-----------------------------------------------------------------------
    42 c  radius(nqmx)      ! aerosol particle radius (m)
    43 c  rho_q(nqmx)       ! tracer densities (kg.m-3)
    44 c  qext(nqmx)        ! Single Scat. Extinction coeff at 0.67 um
    45 c  alpha_lift(nqmx)  ! saltation vertical flux/horiz flux ratio (m-1)
    46 c  alpha_devil(nqmx) ! lifting coeeficient by dust devil
     43c  radius(nq)      ! aerosol particle radius (m)
     44c  rho_q(nq)       ! tracer densities (kg.m-3)
     45c  qext(nq)        ! Single Scat. Extinction coeff at 0.67 um
     46c  alpha_lift(nq)  ! saltation vertical flux/horiz flux ratio (m-1)
     47c  alpha_devil(nq) ! lifting coeeficient by dust devil
    4748c  rho_dust          ! Mars dust density
    4849c  rho_ice           ! Water ice density
    … …  
    5152c-----------------------------------------------------------------------
    5253
     54       !! we allocate once for all arrays in common in tracer_h.F90
     55       !! (supposedly those are not used before call to initracer)
     56       ALLOCATE(noms(nq))
     57       ALLOCATE(mmol(nq))
     58       ALLOCATE(radius(nq))
     59       ALLOCATE(rho_q(nq))
     60       ALLOCATE(qext(nq))
     61       ALLOCATE(alpha_lift(nq))
     62       ALLOCATE(alpha_devil(nq))
     63       ALLOCATE(qextrhor(nq))
     64       ALLOCATE(igcm_dustbin(nq))
     65       !! initialization
     66       alpha_lift(:)=0.
     67       alpha_devil(:)=0.
     68
    5369! Initialization: get tracer names from the dynamics and check if we are
    5470!                 using 'old' tracer convention ('q01',q02',...)
    … …  
    5773
    5874! copy tracer names from dynamics
    59         do iq=1,nqmx
    60           noms(iq)=tnom(iq)
     75        do iq=1,nq
     76          noms(iq)=nametrac(iq)
    6177        enddo
    6278
    … …  
    6581      ! 0. initialize tracer indexes to zero:
    6682      ! NB: igcm_* indexes are commons in 'tracer.h'
    67       do iq=1,nqmx
     83      do iq=1,nq
    6884        igcm_dustbin(iq)=0
    6985      enddo
    … …  
    97113      count=0
    98114!      if (dustbin.gt.0) then
    99 !        do iq=1,nqmx
     115!        do iq=1,nq
    100116!          txt=" "
    101117!          write(txt,'(a4,i2.2)')'dust',count+1   
    … …  
    105121!            mmol(iq)=100.
    106122!          endif
    107 !        enddo !do iq=1,nqmx
     123!        enddo !do iq=1,nq
    108124!      endif ! of if (dustbin.gt.0)
    109125
    110126
    111127!      if (doubleq) then
    112 !        do iq=1,nqmx
     128!        do iq=1,nq
    113129!          if (noms(iq).eq."dust_mass") then
    114130!            igcm_dust_mass=iq
    … …  
    122138!      endif ! of if (doubleq)
    123139      ! 2. find chemistry and water tracers
    124       do iq=1,nqmx
     140      do iq=1,nq
    125141        if (noms(iq).eq."co2") then
    126142          igcm_co2=iq
    … …  
    147163!          write(*,*) 'h2o_ice: count=',count
    148164        endif
    149       enddo ! of do iq=1,nqmx
     165      enddo ! of do iq=1,nq
    150166     
    151167      ! check that we identified all tracers:
    152       if (count.ne.nqmx) then
     168      if (count.ne.nq) then
    153169        write(*,*) "initracer: found only ",count," tracers"
    154         write(*,*) "               expected ",nqmx
     170        write(*,*) "               expected ",nq
    155171        do iq=1,count
    156172          write(*,*)'      ',iq,' ',trim(noms(iq))
    … …  
    159175      else
    160176        write(*,*) "initracer: found all expected tracers, namely:"
    161         do iq=1,nqmx
     177        do iq=1,nq
    162178          write(*,*)'      ',iq,' ',trim(noms(iq))
    163179        enddo
    … …  
    168184c     Initialisation tracers ....
    169185c------------------------------------------------------------
    170       call zerophys(nqmx,rho_q)
     186      call zerophys(nq,rho_q)
    171187
    172188      rho_dust=2500.  ! Mars dust density (kg.m-3)
    … …  
    182198c$$$c       iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
    183199c$$$
    184 c$$$        if( (nqmx.lt.2).or.(water.and.(nqmx.lt.3)) ) then
    185 c$$$          write(*,*)'initracer: nqmx is too low : nqmx=', nqmx
     200c$$$        if( (nq.lt.2).or.(water.and.(nq.lt.3)) ) then
     201c$$$          write(*,*)'initracer: nq is too low : nq=', nq
    186202c$$$          write(*,*)'water= ',water,' doubleq= ',doubleq   
    187203c$$$        end if
    … …  
    259275
    260276
    261 !         if(ngridmx.eq.1)
     277!         if(ngrid.eq.1)
    262278
    263279
    264280!     to be modified for BC+ version?
    265          do ig=1,ngridmx
    266            if (ngridmx.ne.1) watercaptag(ig)=.false.
     281
     282         !! this is defined in surfdat_h.F90
     283         IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
     284         IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
     285
     286         do ig=1,ngrid
     287           if (ngrid.ne.1) watercaptag(ig)=.false.
    267288           dryness(ig) = 1.
    268289         enddo
    … …  
    274295c Perennial H20 north cap defined by watercaptag=true (allows surface to be
    275296c hollowed by sublimation in vdifc).
    276 !         do ig=1,ngridmx
     297!         do ig=1,ngrid
    277298!           if (lati(ig)*180./pi.gt.84) then
    278 !             if (ngridmx.ne.1) watercaptag(ig)=.true.
     299!             if (ngrid.ne.1) watercaptag(ig)=.true.
    279300!             dryness(ig) = 1.
    280301c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
    281302c             if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
    282 c               if (ngridmx.ne.1) watercaptag(ig)=.true.
     303c               if (ngrid.ne.1) watercaptag(ig)=.true.
    283304c               dryness(ig) = 1.
    284305c             endif
    285306c             if (lati(ig)*180./pi.ge.85) then
    286 c               if (ngridmx.ne.1) watercaptag(ig)=.true.
     307c               if (ngrid.ne.1) watercaptag(ig)=.true.
    287308c               dryness(ig) = 1.
    288309c             endif
    289310!           endif  ! (lati>80 deg)
    290 !         end do ! (ngridmx)
     311!         end do ! (ngrid)
    291312!        ENDIF ! (caps)
    292313
    293 !         if(iceparty.and.(nqmx.ge.2)) then
     314!         if(iceparty.and.(nq.ge.2)) then
    294315
    295316           radius(igcm_h2o_ice)=3.e-6
    … …  
    303324
    304325
    305 !         elseif(iceparty.and.(nqmx.lt.2)) then
    306 !            write(*,*) 'nqmx is too low : nqmx=', nqmx
     326!         elseif(iceparty.and.(nq.lt.2)) then
     327!            write(*,*) 'nq is too low : nq=', nq
    307328!            write(*,*) 'water= ',water,' iceparty= ',iceparty
    308329!         endif

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

    r135 r787  
    11      SUBROUTINE iniwrite(nid,idayref,phis)
     2
     3      USE comsoil_h
     4
    25      IMPLICIT NONE
    36
    … …  
    3134#include "serre.h"
    3235#include"dimphys.h"
    33 #include"comsoil.h"
    3436
    3537c   Arguments:

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

    r533 r787  
    33      use radinc_h, only: L_NSPECTI
    44      use radcommon_h, only: WNOI,DWNI
     5      use comsoil_h
    56
    67      implicit none
    … …  
    3536#include "serre.h"
    3637#include"dimphys.h"
    37 #include"comsoil.h"
    3838
    3939c   Arguments:

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

    r533 r787  
    33      use radinc_h, only: L_NSPECTV
    44      use radcommon_h, only: WNOV,DWNV
     5      use comsoil_h
    56
    67      implicit none
    … …  
    3536#include "serre.h"
    3637#include"dimphys.h"
    37 #include"comsoil.h"
    3838
    3939c   Arguments:

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

    r135 r787  
    1 subroutine iniwritesoil(nid)
     1subroutine iniwritesoil(nid,ngrid)
     2
     3use comsoil_h
    24
    35! initialization routine for 'writediagoil'. Here we create/define
    … …  
    1214#include"comcstfi.h"
    1315#include"comgeom.h"
    14 #include"comsoil.h"
    1516#include"netcdf.inc"
    1617
    1718! Arguments:
     19integer,intent(in) :: ngrid
    1820integer,intent(in) :: nid ! NetCDF output file ID
    1921
    … …  
    246248  do i=1,iip1
    247249    data3(i,1,l)=inertiedat(1,l)
    248     data3(i,jjp1,l)=inertiedat(ngridmx,l)
     250    data3(i,jjp1,l)=inertiedat(ngrid,l)
    249251  enddo
     252  !!! THIS WILL NOT WORK IN PARALLEL !!!!
    250253  ! rest of the grid
    251254  do j=2,jjm

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

    r716 r787  
    55  use watercommon_h, only: mH2O
    66  use ioipsl_getincom
     7  use comsaison_h
    78
    89  implicit none
    … …  
    3132#include "comcstfi.h"
    3233#include "planete.h"
    33 #include "comsaison.h"
    3434#include "control.h"
    35 #include "advtrac.h"
    3635
    3736  ! --------------------------------------------------------------
    … …  
    7170  real dTstrat
    7271  real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
    73   real OLR_nu(ngridmx,L_NSPECTI)
    74   real OSR_nu(ngridmx,L_NSPECTV)
     72  real OLR_nu(1,L_NSPECTI)
     73  real OSR_nu(1,L_NSPECTV)
    7574  real Eatmtot
    7675
    7776  integer ierr
    7877  logical firstcall,lastcall,global1d
     78
     79  character*20 nametrac(nqmx)   ! name of the tracer (no need for adv trac common)
    7980
    8081  ! --------------
    … …  
    9293  nlayer=nlayermx
    9394  nlevel=nlayer+1
     95
     96  !! this is defined in comsaison_h
     97  ALLOCATE(mu0(1))
     98  ALLOCATE(fract(1))
     99
    94100
    95101
    … …  
    224230        do iq=1,nq
    225231           ! minimal version, just read in the tracer names, 1 per line
    226            read(90,*,iostat=ierr) tnom(iq)
     232           read(90,*,iostat=ierr) nametrac(iq)
    227233           if (ierr.ne.0) then
    228234              write(*,*) 'kcm1d: error reading tracer names...'
    … …  
    232238     endif
    233239
    234      call initracer()
     240     call initracer(1,nq,nametrac)
    235241
    236242  endif
    237243
    238244
    239   do iq=1,nqmx
     245  do iq=1,nq
    240246     do ilay=1,nlayer
    241247        q(ilay,iq) = 0.
    … …  
    243249  enddo
    244250
    245   do iq=1,nqmx
     251  do iq=1,nq
    246252     qsurf(iq) = 0.
    247253  enddo
    … …  
    265271
    266272     !    Run radiative transfer
    267      call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
     273     call callcorrk(1,nlayer,q,nq,qsurf,      &
    268274          albedo,emis,mu0,plev,play,temp,                    &
    269275          tsurf,fract,dist_star,aerosol,muvar,         &
    270276          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    271           fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,tau_col,  &
     277          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,nueffrad,tau_col,  &
    272278          cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
    273279
    … …  
    298304  firstcall=.false.
    299305  lastcall=.true.
    300   call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
     306  call callcorrk(1,nlayer,q,nq,qsurf,      &
    301307       albedo,emis,mu0,plev,play,temp,                    &
    302308       tsurf,fract,dist_star,aerosol,muvar,         &
    303309       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    304310       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,             &
    305        reffrad,tau_col,  &
     311       reffrad,nueffrad,tau_col,  &
    306312       cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
    307313

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

    r786 r787  
    1       subroutine largescale(ptimestep, pplev, pplay, pt, pq,   &
     1      subroutine largescale(ngrid,nq,ptimestep, pplev, pplay, pt, pq,   &
    22                        pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)
    33
    44      use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
    55          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
     6      USE tracer_h
    67      IMPLICIT none
    78
    … …  
    2324#include "comcstfi.h"
    2425
    25 #include "fisice.h"
    2626#include "callkeys.h"
    27 #include "tracer.h"
    2827
    29 
     28      INTEGER ngrid,nq
    3029
    3130!     Arguments
    3231      REAL ptimestep                 ! intervalle du temps (s)
    33       REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche
    34       REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche
    35       REAL pt(ngridmx,nlayermx)      ! temperature (K)
    36       real pq(ngridmx,nlayermx,nqmx) ! tracer mixing ratio (kg/kg)
    37       REAL pdt(ngridmx,nlayermx)     ! physical temperature tenedency (K/s)
    38       REAL pdq(ngridmx,nlayermx,nqmx)! physical tracer tenedency (K/s)
    39       REAL pdtlsc(ngridmx,nlayermx)  ! incrementation de la temperature (K)
    40       REAL pdqvaplsc(ngridmx,nlayermx) ! incrementation de la vapeur d'eau
    41       REAL pdqliqlsc(ngridmx,nlayermx) ! incrementation de l'eau liquide
    42       REAL rneb(ngridmx,nlayermx)    ! fraction nuageuse
     32      REAL pplev(ngrid,nlayermx+1) ! pression a inter-couche
     33      REAL pplay(ngrid,nlayermx)   ! pression au milieu de couche
     34      REAL pt(ngrid,nlayermx)      ! temperature (K)
     35      real pq(ngrid,nlayermx,nq) ! tracer mixing ratio (kg/kg)
     36      REAL pdt(ngrid,nlayermx)     ! physical temperature tenedency (K/s)
     37      REAL pdq(ngrid,nlayermx,nq)! physical tracer tenedency (K/s)
     38      REAL pdtlsc(ngrid,nlayermx)  ! incrementation de la temperature (K)
     39      REAL pdqvaplsc(ngrid,nlayermx) ! incrementation de la vapeur d'eau
     40      REAL pdqliqlsc(ngrid,nlayermx) ! incrementation de l'eau liquide
     41      REAL rneb(ngrid,nlayermx)    ! fraction nuageuse
    4342
    4443
    … …  
    5352      INTEGER,PARAMETER :: nitermax=1000
    5453      REAL,PARAMETER :: alpha=.5,qthreshold=1.e-6
    55       REAL zt(ngridmx), zq(ngridmx)
    56       REAL zcond(ngridmx),zcond_iter
    57       REAL zdelq(ngridmx)
    58       REAL zqs(ngridmx), zdqs(ngridmx)
     54      REAL zt(ngrid), zq(ngrid)
     55      REAL zcond(ngrid),zcond_iter
     56      REAL zdelq(ngrid)
     57      REAL zqs(ngrid), zdqs(ngrid)
    5958      REAL psat_tmp
    6059     
    6160! evaporation calculations
    62       REAL dqevap(ngridmx,nlayermx),dtevap(ngridmx,nlayermx)     
    63       REAL qevap(ngridmx,nlayermx,nqmx)
    64       REAL tevap(ngridmx,nlayermx)
     61      REAL dqevap(ngrid,nlayermx),dtevap(ngrid,nlayermx)     
     62      REAL qevap(ngrid,nlayermx,nq)
     63      REAL tevap(ngrid,nlayermx)
    6564
    66       REAL zcor(ngridmx), zdelta(ngridmx), zcvm5(ngridmx)
    67       REAL zx_q(ngridmx)
     65      REAL zcor(ngrid), zdelta(ngrid), zcvm5(ngrid)
     66      REAL zx_q(ngrid)
    6867      REAL Nmix_local,zfice
    6968
    7069!     GCM -----> subroutine variables, initialisation of outputs
    7170
    72       pdtlsc(1:ngridmx,1:nlayermx)  = 0.0
    73       pdqvaplsc(1:ngridmx,1:nlayermx)  = 0.0
    74       pdqliqlsc(1:ngridmx,1:nlayermx) = 0.0
    75       rneb(1:ngridmx,1:nlayermx) = 0.0
     71      pdtlsc(1:ngrid,1:nlayermx)  = 0.0
     72      pdqvaplsc(1:ngrid,1:nlayermx)  = 0.0
     73      pdqliqlsc(1:ngrid,1:nlayermx) = 0.0
     74      rneb(1:ngrid,1:nlayermx) = 0.0
    7675
    7776
    7877      ! Evaporate cloud water/ice
    79       call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
     78      call evap(ngrid,nq,ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
    8079      ! note: we use qevap but not tevap in largescale/moistadj
    8180            ! otherwise is a big mess
    … …  
    8584   DO k = nlayermx, 1, -1
    8685
    87       zt(1:ngridmx)=pt(1:ngridmx,k)+(pdt(1:ngridmx,k)+dtevap(1:ngridmx,k))*ptimestep
    88       zq(1:ngridmx)=qevap(1:ngridmx,k,igcm_h2o_vap) !liquid water is included in qevap
     86      zt(1:ngrid)=pt(1:ngrid,k)+(pdt(1:ngrid,k)+dtevap(1:ngrid,k))*ptimestep
     87      zq(1:ngrid)=qevap(1:ngrid,k,igcm_h2o_vap) !liquid water is included in qevap
    8988
    9089!     Calculer la vapeur d'eau saturante et
    9190!     determiner la condensation partielle
    92       DO i = 1, ngridmx
     91      DO i = 1, ngrid
    9392
    9493         if(zt(i).le.15.) then
    … …  
    150149
    151150!     Tendances de t et q
    152          pdqvaplsc(1:ngridmx,k)  = dqevap(1:ngridmx,k) - zcond(1:ngridmx)
    153          pdqliqlsc(1:ngridmx,k) = - pdqvaplsc(1:ngridmx,k)
    154          pdtlsc(1:ngridmx,k)  = pdqliqlsc(1:ngridmx,k)*RLVTT/RCPD
     151         pdqvaplsc(1:ngrid,k)  = dqevap(1:ngrid,k) - zcond(1:ngrid)
     152         pdqliqlsc(1:ngrid,k) = - pdqvaplsc(1:ngrid,k)
     153         pdtlsc(1:ngrid,k)  = pdqliqlsc(1:ngrid,k)*RLVTT/RCPD
    155154
    156155   Enddo ! k= nlayermx, 1, -1

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

    r786 r787  
    55                                                   
    66       USE watercommon_h, Only: Tsat_water,RLVTT
     7       use surfdat_h
     8       USE comgeomfi_h
     9       USE tracer_h
    710
    811       IMPLICIT NONE
    … …  
    4952#include "dimphys.h"
    5053#include "comcstfi.h"
    51 #include "surfdat.h"
    52 #include "comgeomfi.h"
    5354#include "comvert.h"
    5455#include "paramet.h"
    5556#include "callkeys.h"
    56 #include "tracer.h"
    5757
    5858!-----------------------------------------------------------------------
    … …  
    6161      INTEGER ngrid, nlayer, nq   
    6262      REAL ptimestep
    63       REAL pcapcal(ngridmx)
    64       INTEGER rnat(ngridmx)     
     63      REAL pcapcal(ngrid)
     64      INTEGER rnat(ngrid)     
    6565      REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
    6666      REAL pt(ngrid,nlayer),pdt(ngrid,nlayer)
    … …  
    7171      REAL pdmassmr(ngrid,nlayer)
    7272      REAL pdumr(ngrid,nlayer) , pdvmr(ngrid,nlayer)
    73       REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq)
    74       REAL pqs(ngridmx,nq)
     73      REAL pq(ngrid,nlayer,nq),pdq(ngrid,nlayer,nq)
     74      REAL pqs(ngrid,nq)
    7575      REAL pdqmr(ngrid,nlayer,nq),pdqsmr(ngrid,nq)
    7676      REAL pdpsrfmr(ngrid),pdtsrfmr(ngrid)
    … …  
    8080
    8181!    Boiling/sublimation
    82       REAL Tsat(ngridmx),zmassboil(ngridmx)
     82      REAL Tsat(ngrid),zmassboil(ngrid)
    8383
    8484!    vertical reorganization of sigma levels
    … …  
    8787!    Dummy variables     
    8888      INTEGER n,l,ig,iq
    89       REAL zdtsig(ngridmx,nlayermx)
    90       REAL zmass(ngridmx,nlayermx),zzmass(nlayermx),w(nlayermx+1)
    91       REAL zdmass_sum(ngridmx,nlayermx+1)
     89      REAL zdtsig(ngrid,nlayermx)
     90      REAL zmass(ngrid,nlayermx),zzmass(nlayermx),w(nlayermx+1)
     91      REAL zdmass_sum(ngrid,nlayermx+1)
    9292      REAL zmflux(nlayermx+1)
    9393      REAL zq1(nlayermx)
    94       REAL ztsrf(ngridmx)
     94      REAL ztsrf(ngrid)
    9595      REAL ztm(nlayermx+1)
    9696      REAL zum(nlayermx+1) , zvm(nlayermx+1)
    97       REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1)
     97      REAL zqm(nlayermx+1,nq),zqm1(nlayermx+1)
    9898
    9999!   local saved variables
    … …  
    120120!     Surface tracer Tendencies set to 0
    121121!     -------------------------------------
    122       pdqsmr(1:ngridmx,1:nqmx)=0.
    123 
    124       ztsrf(1:ngridmx) = ptsrf(1:ngridmx) + pdtsrf(1:ngridmx)*ptimestep
     122      pdqsmr(1:ngrid,1:nq)=0.
     123
     124      ztsrf(1:ngrid) = ptsrf(1:ngrid) + pdtsrf(1:ngrid)*ptimestep
    125125
    126126
    … …  
    135135
    136136      if (water) then
    137          do ig=1,ngridmx
     137         do ig=1,ngrid
    138138            call Tsat_water(pplev(ig,1)+zdmass_sum(ig,1)*g*ptimestep,Tsat(ig))
    139139         enddo
    140140               call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat)
    141141         
    142          do ig=1,ngridmx
     142         do ig=1,ngrid
    143143            if (ztsrf(ig).gt.Tsat(ig)) then
    144144               zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT/ptimestep
    … …  
    163163!    """"""""""""""""""""""""""""""""""""
    164164         
    165       pdpsrfmr(1:ngridmx) = (zdmass_sum(1:ngridmx,1)+zmassboil(1:ngridmx))*g
    166 
    167       do ig = 1, ngridmx
     165      pdpsrfmr(1:ngrid) = (zdmass_sum(1:ngrid,1)+zmassboil(1:ngrid))*g
     166
     167      do ig = 1, ngrid
    168168        IF(ABS(pdpsrfmr(ig)*ptimestep).GT.pplev(ig,1)) THEN
    169169         PRINT*,'STOP in condens'
    … …  
    187187         zzu(1:nlayermx)  = pu(ig,1:nlayermx) + pdu(ig,1:nlayermx) * ptimestep
    188188         zzv(1:nlayermx)  = pv(ig,1:nlayermx) + pdv(ig,1:nlayermx) * ptimestep
    189          zzq(1:nlayermx,1:nqmx)=pq(ig,1:nlayermx,1:nqmx)+pdq(ig,1:nlayermx,1:nqmx)*ptimestep ! must add the water that has fallen???
     189         zzq(1:nlayermx,1:nq)=pq(ig,1:nlayermx,1:nq)+pdq(ig,1:nlayermx,1:nq)*ptimestep ! must add the water that has fallen???
    190190
    191191!  Mass fluxes of air through the sigma levels (kg.m-2.s-1)  (>0 when up)
    … …  
    214214         zum(1) = 0. 
    215215         zvm(1) = 0. 
    216          zqm(1,1:nqmx)=0. ! most tracer do not condense !
     216         zqm(1,1:nq)=0. ! most tracer do not condense !
    217217         if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface
    218218         
    … …  
    222222         call vl1d(zzu ,2.,zzmass,w,zum)
    223223         call vl1d(zzv ,2.,zzmass,w,zvm)
    224          do iq=1,nqmx
     224         do iq=1,nq
    225225           zq1(1:nlayermx)=zzq(1:nlayermx,iq)
    226226           zqm1(1)=zqm(1,iq)
    … …  
    247247         zum(nlayer+1)= zzu(nlayer)  ! should not be used, but...
    248248         zvm(nlayer+1)= zzv(nlayer)  ! should not be used, but...
    249          zqm(nlayer+1,1:nqmx)= zzq(nlayer,1:nqmx)
     249         zqm(nlayer+1,1:nq)= zzq(nlayer,1:nq)
    250250 
    251251!        Tendencies on T, U, V, Q
    … …  
    267267
    268268!        Tendencies on Q
    269          do iq=1,nqmx
     269         do iq=1,nq
    270270            DO l=1,nlayer
    271271               pdqmr(ig,l,iq)= (1/zzmass(l)) *   &

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

    r773 r787  
    1 subroutine moistadj(pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
     1subroutine moistadj(ngrid, nq, pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
    22
    33  use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
     4  USE tracer_h
    45
    56  implicit none
    … …  
    2122#include "dimensions.h"
    2223#include "dimphys.h"
    23 #include "tracer.h"
    2424#include "comcstfi.h"
    2525
    26       REAL pt(ngridmx,nlayermx)            ! temperature (K)
    27       REAL pq(ngridmx,nlayermx,nqmx)       ! tracer (kg/kg)
    28       REAL pdq(ngridmx,nlayermx,nqmx)
    29 
    30       REAL pdqmana(ngridmx,nlayermx,nqmx)  ! tendency of tracers (kg/kg.s-1)
    31       REAL pdtmana(ngridmx,nlayermx)       ! temperature increment
     26      INTEGER ngrid, nq
     27
     28      REAL pt(ngrid,nlayermx)            ! temperature (K)
     29      REAL pq(ngrid,nlayermx,nq)       ! tracer (kg/kg)
     30      REAL pdq(ngrid,nlayermx,nq)
     31
     32      REAL pdqmana(ngrid,nlayermx,nq)  ! tendency of tracers (kg/kg.s-1)
     33      REAL pdtmana(ngrid,nlayermx)       ! temperature increment
    3234
    3335!     local variables
    34       REAL zt(ngridmx,nlayermx)      ! temperature (K)
    35       REAL zq(ngridmx,nlayermx)      ! humidite specifique (kg/kg)
    36       REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa)
    37       REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche (Pa)
    38 
    39       REAL d_t(ngridmx,nlayermx)     ! temperature increment
    40       REAL d_q(ngridmx,nlayermx)     ! incrementation pour vapeur d'eau
    41       REAL d_ql(ngridmx,nlayermx)    ! incrementation pour l'eau liquide
    42       REAL rneb(ngridmx,nlayermx) ! cloud fraction
     36      REAL zt(ngrid,nlayermx)      ! temperature (K)
     37      REAL zq(ngrid,nlayermx)      ! humidite specifique (kg/kg)
     38      REAL pplev(ngrid,nlayermx+1) ! pression a inter-couche (Pa)
     39      REAL pplay(ngrid,nlayermx)   ! pression au milieu de couche (Pa)
     40
     41      REAL d_t(ngrid,nlayermx)     ! temperature increment
     42      REAL d_q(ngrid,nlayermx)     ! incrementation pour vapeur d'eau
     43      REAL d_ql(ngrid,nlayermx)    ! incrementation pour l'eau liquide
     44      REAL rneb(ngrid,nlayermx) ! cloud fraction
    4345      REAL ptimestep
    4446
    … …  
    5254      INTEGER, PARAMETER :: niter = 1
    5355      INTEGER k1, k1p, k2, k2p
    54       LOGICAL itest(ngridmx)
    55       REAL delta_q(ngridmx, nlayermx)
     56      LOGICAL itest(ngrid)
     57      REAL delta_q(ngrid, nlayermx)
    5658      REAL cp_new_t(nlayermx)
    5759      REAL cp_delta_t(nlayermx)
    5860      REAL v_cptj(nlayermx), v_cptjk1, v_ssig
    59       REAL v_cptt(ngridmx,nlayermx), v_p, v_t, zpsat
    60       REAL zqs(ngridmx,nlayermx), zdqs(ngridmx,nlayermx)
    61       REAL zq1(ngridmx), zq2(ngridmx)
    62       REAL gamcpdz(ngridmx,2:nlayermx)
     61      REAL v_cptt(ngrid,nlayermx), v_p, v_t, zpsat
     62      REAL zqs(ngrid,nlayermx), zdqs(ngrid,nlayermx)
     63      REAL zq1(ngrid), zq2(ngrid)
     64      REAL gamcpdz(ngrid,2:nlayermx)
    6365      REAL zdp, zdpm
    6466
    … …  
    6668      REAL zflo ! flotabilite
    6769
    68       REAL local_q(ngridmx,nlayermx),local_t(ngridmx,nlayermx)
     70      REAL local_q(ngrid,nlayermx),local_t(ngrid,nlayermx)
    6971
    7072      REAL zdelta, zcor, zcvm5
    … …  
    9496
    9597!     GCM -----> subroutine variables
    96       zq(1:ngridmx,1:nlayermx)    = pq(1:ngridmx,1:nlayermx,i_h2o)+ pdq(1:ngridmx,1:nlayermx,i_h2o)*ptimestep
    97       zt(1:ngridmx,1:nlayermx)    = pt(1:ngridmx,1:nlayermx)
    98       pdqmana(1:ngridmx,1:nlayermx,1:nqmx)=0.0
    99 
    100       DO k = 1, nlayermx
    101        DO i = 1, ngridmx
     98      zq(1:ngrid,1:nlayermx)    = pq(1:ngrid,1:nlayermx,i_h2o)+ pdq(1:ngrid,1:nlayermx,i_h2o)*ptimestep
     99      zt(1:ngrid,1:nlayermx)    = pt(1:ngrid,1:nlayermx)
     100      pdqmana(1:ngrid,1:nlayermx,1:nq)=0.0
     101
     102      DO k = 1, nlayermx
     103       DO i = 1, ngrid
    102104         if(zq(i,k).lt.0.)then
    103105            zq(i,k)=0.0
    … …  
    106108      ENDDO
    107109     
    108       local_q(1:ngridmx,1:nlayermx) = zq(1:ngridmx,1:nlayermx)
    109       local_t(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
    110       rneb(1:ngridmx,1:nlayermx) = 0.0
    111       d_ql(1:ngridmx,1:nlayermx) = 0.0
    112       d_t(1:ngridmx,1:nlayermx)  = 0.0
    113       d_q(1:ngridmx,1:nlayermx)  = 0.0
     110      local_q(1:ngrid,1:nlayermx) = zq(1:ngrid,1:nlayermx)
     111      local_t(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx)
     112      rneb(1:ngrid,1:nlayermx) = 0.0
     113      d_ql(1:ngrid,1:nlayermx) = 0.0
     114      d_t(1:ngrid,1:nlayermx)  = 0.0
     115      d_q(1:ngrid,1:nlayermx)  = 0.0
    114116
    115117!     Calculate v_cptt
    116118      DO k = 1, nlayermx
    117          DO i = 1, ngridmx
     119         DO i = 1, ngrid
    118120            v_cptt(i,k) = RCPD * local_t(i,k)
    119121            v_t = MAX(local_t(i,k),15.)
    … …  
    127129!     Calculate Gamma * Cp * dz: (gamma is the critical gradient)
    128130      DO k = 2, nlayermx
    129          DO i = 1, ngridmx
     131         DO i = 1, ngrid
    130132            zdp = pplev(i,k)-pplev(i,k+1)
    131133            zdpm = pplev(i,k-1)-pplev(i,k)
    … …  
    138140
    139141!------------------------------------ modification of unstable profile
    140       DO 9999 i = 1, ngridmx
     142      DO 9999 i = 1, ngrid
     143
    141144      itest(i) = .FALSE.
    142145
    … …  
    266269
    267270      DO k = 1, nlayermx
    268       DO i = 1, ngridmx
     271      DO i = 1, ngrid
    269272         IF (itest(i)) THEN
    270273         delta_q(i,k) = local_q(i,k) - zq(i,k)
    … …  
    278281! diminue et d'une maniere proportionnelle a cet diminution):
    279282
    280       DO i = 1, ngridmx
     283      DO i = 1, ngrid
    281284         IF (itest(i)) THEN
    282285         zq1(i) = 0.0
    … …  
    285288      ENDDO
    286289      DO k = 1, nlayermx
    287       DO i = 1, ngridmx
     290      DO i = 1, ngrid
    288291         IF (itest(i)) THEN
    289292         zdp = pplev(i,k)-pplev(i,k+1)
    … …  
    294297      ENDDO
    295298      DO k = 1, nlayermx
    296       DO i = 1, ngridmx
     299      DO i = 1, ngrid
    297300         IF (itest(i)) THEN
    298301           IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
    … …  
    304307
    305308      DO k = 1, nlayermx
    306       DO i = 1, ngridmx
     309      DO i = 1, ngrid
    307310          local_q(i, k) = MAX(local_q(i, k), seuil_vap)
    308311      ENDDO
    … …  
    310313
    311314      DO k = 1, nlayermx
    312       DO i = 1, ngridmx
     315      DO i = 1, ngrid
    313316         d_t(i,k) = local_t(i,k) - zt(i,k)
    314317         d_q(i,k) = local_q(i,k) - zq(i,k)
    … …  
    318321!     now subroutine -----> GCM variables
    319322      DO k = 1, nlayermx
    320          DO i = 1, ngridmx
     323         DO i = 1, ngrid
    321324           
    322325            pdtmana(i,k)       = d_t(i,k)/ptimestep

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

    r486 r787  
    3535c    Traceurs :
    3636      real pqi(ngrid,nlay)  ! traceur   (e.g. ?/kg)
    37       real wq(ngridmx,nlay+1)  ! flux de traceur durant timestep (?/m-2)
     37      real wq(ngrid,nlay+1)  ! flux de traceur durant timestep (?/m-2)
    3838     
    3939c   local:
    … …  
    4848c    Traceurs :
    4949c    ~~~~~~~~
    50       real traversee (ngridmx,nlayermx)
    51       real vstokes(ngridmx,nlayermx)
    52       real w(ngridmx,nlayermx)
     50      real traversee (ngrid,nlayermx)
     51      real vstokes(ngrid,nlayermx)
     52      real w(ngrid,nlayermx)
    5353      real ptop, dztop, Ep, Stra
    5454
    … …  
    7474
    7575      IF (firstcall) THEN
    76          IF(ngrid.NE.ngridmx) THEN
    77             PRINT*,'STOP dans newsedim'
    78             PRINT*,'probleme de dimensions :'
    79             PRINT*,'ngrid  =',ngrid
    80             PRINT*,'ngridmx  =',ngridmx
    81             STOP
    82          ENDIF
    8376         firstcall=.false.
    8477

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

    r253 r787  
    1 subroutine newtrelax(mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall)
     1subroutine newtrelax(ngrid,mu0,sinlat,popsk,temp,pplay,pplev,dtrad,firstcall)
    22       
    33  implicit none
    … …  
    2020!     
    2121!==================================================================
    22  
     22 
     23  integer ngrid
     24 
    2325  ! Input
    24   real mu0(ngridmx)                    ! cosine of sun incident angle
    25   real sinlat(ngridmx)                 ! sine of latitude
    26   real temp(ngridmx,nlayermx)          ! temperature at each layer (K)
    27   real pplay(ngridmx,nlayermx)         ! pressure at each layer (Pa)
    28   real pplev(ngridmx,nlayermx+1)       ! pressure at each level (Pa)
    29   real popsk(ngridmx,nlayermx)         ! pot. T to T converter
     26  real mu0(ngrid)                    ! cosine of sun incident angle
     27  real sinlat(ngrid)                 ! sine of latitude
     28  real temp(ngrid,nlayermx)          ! temperature at each layer (K)
     29  real pplay(ngrid,nlayermx)         ! pressure at each layer (Pa)
     30  real pplev(ngrid,nlayermx+1)       ! pressure at each level (Pa)
     31  real popsk(ngrid,nlayermx)         ! pot. T to T converter
    3032
    3133  ! Output
    32   real dtrad(ngridmx,nlayermx)
     34  real dtrad(ngrid,nlayermx)
    3335
    3436  ! Internal
    35   real Trelax(ngridmx,nlayermx), Trelax_V, Trelax_H
    36   save Trelax
     37  real Trelax_V, Trelax_H
     38  real,allocatable,dimension(:,:),save :: Trelax
    3739
    3840  real T_trop ! relaxation temperature at tropopause (K)
    … …  
    5153  if(firstcall) then
    5254
     55     ALLOCATE(Trelax(ngrid,nlayermx))
     56
    5357     print*,'-----------------------------------------------------'
    5458     print*,'| ATTENTION: You are using a Newtonian cooling scheme'
    … …  
    5862
    5963     if(tidallocked)then
    60         do ig=1,ngridmx
     64        do ig=1,ngrid
    6165
    6266           T_surf = 126. + 239.*mu0(ig)
    … …  
    8387
    8488        do l=1,nlayermx
    85            do ig=1,ngridmx
     89           do ig=1,ngrid
    8690
    8791              ! vertically varying component
    … …  
    111115  ! Calculate radiative forcing
    112116  do l=1,nlayermx
    113      do ig=1,ngridmx
     117     do ig=1,ngrid
    114118        dtrad(ig,l) = -(temp(ig,l) - Trelax(ig,l)) / tau_relax
    115119        if(temp(ig,l).gt.500.)then ! Trelax(ig,l))then
    … …  
    122126  enddo
    123127
    124   call writediagfi(ngridmx,'Tref','rad forc temp','K',3,Trelax)
    125   !call writediagfi(ngridmx,'ThetaZ','stellar zenith angle','deg',2,mu0)
     128  call writediagfi(ngrid,'Tref','rad forc temp','K',3,Trelax)
     129  !call writediagfi(ngrid,'ThetaZ','stellar zenith angle','deg',2,mu0)
    126130
    127131  return

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

    r728 r787  
    1       SUBROUTINE phyetat0 (fichnom,tab0,Lmodif,nsoil,nq,
     1      SUBROUTINE phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq,
    22     .           day_ini,time,
    33     .           tsurf,tsoil,emis,q2,qsurf,cloudfrac,totcloudfrac,hice)
     4
     5      USE surfdat_h
     6      USE comgeomfi_h
     7      USE tracer_h
     8
    49      implicit none
     10
    511c======================================================================
    612c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
    … …  
    1117#include "dimensions.h"
    1218#include "dimphys.h"
    13 #include "comgeomfi.h"
    14 #include "surfdat.h"
    1519#include "planete.h"
    1620#include "comcstfi.h"
    17 #include"advtrac.h"
    18 
     21
     22      INTEGER ngrid
    1923c======================================================================
    2024      INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
    … …  
    3337
    3438!  outputs:
    35       real tsurf(ngridmx,nbsrf) ! surface temperature
    36       real tsoil(ngridmx,nsoil,nbsrf) ! soil temperature
    37       real emis(ngridmx) ! surface emissivity
    38       real q2(ngridmx, llm+1) !
    39       real qsurf(ngridmx,nq) ! tracers on surface
    40 !      real co2ice(ngridmx) ! co2 ice cover
    41       real cloudfrac(ngridmx,nlayermx)
    42       real hice(ngridmx), totcloudfrac(ngridmx)
    43 
     39      real tsurf(ngrid,nbsrf) ! surface temperature
     40      real tsoil(ngrid,nsoil,nbsrf) ! soil temperature
     41      real emis(ngrid) ! surface emissivity
     42      real q2(ngrid, llm+1) !
     43      real qsurf(ngrid,nq) ! tracers on surface
     44!      real co2ice(ngrid) ! co2 ice cover
     45      real cloudfrac(ngrid,nlayermx)
     46      real hice(ngrid), totcloudfrac(ngrid)
    4447
    4548
    … …  
    7174      character(len=30) :: txt ! to store some text
    7275
    73  
     76!! added variables by AS to allow to read only slices of startfi
     77      real :: toto(ngrid)
     78      integer :: sta   !! subscript in starti where we start to read
     79      integer, dimension(2) :: sta2d
     80      integer, dimension(2) :: siz2d
     81
     82c AS: get the cursor that is stored in dimphys.h
     83c ... this allows to read only a part of startfi horiz grid
     84      sta = cursor
     85
     86c
     87c ALLOCATE ARRAYS IN surfdat_h
     88c
     89      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
     90      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
     91      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
     92      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
     93      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
     94      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
     95      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
     96
    7497c
    7598c Ouvrir le fichier contenant l etat initial:
    … …  
    85108!                    qsurf02, ...)
    86109      count=0
    87       do iq=1,nqmx
     110      do iq=1,nq
    88111        txt= " "
    89112        write(txt,'(a5,i2.2)')'qsurf',iq
    … …  
    97120        endif
    98121      enddo
    99       if (count.eq.nqmx) then
     122      if (count.eq.nq) then
    100123        write(*,*) "phyetat0:tracers seem to follow old naming ",
    101124     &             "convention (qsurf01,qsurf02,...)"
    … …  
    106129
    107130c modifications possibles des variables de tab_cntrl
    108       PRINT*
    109131      write(*,*) 'TABFI de phyeta0',Lmodif,tab0
    110       call tabfi (nid,Lmodif,tab0,day_ini,lmax,p_rad,
     132      call tabfi (ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
    111133     .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
     134
    112135c
    113136c Lecture des latitudes (coordonnees):
    … …  
    119142      ENDIF
    120143#ifdef NC_DOUBLE
    121       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, lati)
    122 #else
    123       ierr = NF_GET_VAR_REAL(nid, nvarid, lati)
     144      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
     145#else
     146      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
    124147#endif
    125148      IF (ierr.NE.NF_NOERR) THEN
    … …  
    136159      ENDIF
    137160#ifdef NC_DOUBLE
    138       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, long)
    139 #else
    140       ierr = NF_GET_VAR_REAL(nid, nvarid, long)
     161      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
     162#else
     163      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
    141164#endif
    142165      IF (ierr.NE.NF_NOERR) THEN
    … …  
    153176      ENDIF
    154177#ifdef NC_DOUBLE
    155       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, area)
    156 #else
    157       ierr = NF_GET_VAR_REAL(nid, nvarid, area)
     178      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
     179#else
     180      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
    158181#endif
    159182      IF (ierr.NE.NF_NOERR) THEN
    … …  
    175198      ENDIF
    176199#ifdef NC_DOUBLE
    177       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, phisfi)
    178 #else
    179       ierr = NF_GET_VAR_REAL(nid, nvarid, phisfi)
     200      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,phisfi)
     201#else
     202      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,phisfi)
    180203#endif
    181204      IF (ierr.NE.NF_NOERR) THEN
    … …  
    197220      ENDIF
    198221#ifdef NC_DOUBLE
    199       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, albedodat)
    200 #else
    201       ierr = NF_GET_VAR_REAL(nid, nvarid, albedodat)
     222      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,albedodat)
     223#else
     224      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,albedodat)
    202225#endif
    203226      IF (ierr.NE.NF_NOERR) THEN
    … …  
    211234      PRINT*,'Albedo du sol nu <albedodat>:', xmin, xmax
    212235c
    213 c Lecture de l''inertie thermique du sol:
    214 c
    215 !      ierr = NF_INQ_VARID (nid, "inertiedat", nvarid)
    216 !      IF (ierr.NE.NF_NOERR) THEN
    217 !         PRINT*, 'phyetat0: Le champ <inertiedat> est absent'
    218 !         CALL abort
    219 !      ENDIF
    220 !#ifdef NC_DOUBLE
    221 !      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, inertiedat)
    222 !#else
    223 !      ierr = NF_GET_VAR_REAL(nid, nvarid, inertiedat)
    224 !#endif
    225 !      IF (ierr.NE.NF_NOERR) THEN
    226 !         PRINT*, 'phyetat0: Lecture echouee pour <inertiedat>'
    227 !         CALL abort
    228 !      ENDIF
    229 !      xmin = 1.0E+20
    230 !      xmax = -1.0E+20
    231 !      xmin = MINVAL(inertiedat)
    232 !      xmax = MAXVAL(inertiedat)
    233 !      PRINT*,'Inertie thermique du sol <inertiedat>:', xmin, xmax
    234 c
    235236c ZMEA
    236237c
    … …  
    241242      ENDIF
    242243#ifdef NC_DOUBLE
    243       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zmea)
    244 #else
    245       ierr = NF_GET_VAR_REAL(nid, nvarid, zmea)
     244      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zmea)
     245#else
     246      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zmea)
    246247#endif
    247248      IF (ierr.NE.NF_NOERR) THEN
    … …  
    251252      xmin = 1.0E+20
    252253      xmax = -1.0E+20
    253       DO i = 1, ngridmx
     254      DO i = 1, ngrid
    254255         xmin = MIN(zmea(i),xmin)
    255256         xmax = MAX(zmea(i),xmax)
    … …  
    265266      ENDIF
    266267#ifdef NC_DOUBLE
    267       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zstd)
    268 #else
    269       ierr = NF_GET_VAR_REAL(nid, nvarid, zstd)
     268      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zstd)
     269#else
     270      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zstd)
    270271#endif
    271272      IF (ierr.NE.NF_NOERR) THEN
    … …  
    275276      xmin = 1.0E+20
    276277      xmax = -1.0E+20
    277       DO i = 1, ngridmx
     278      DO i = 1, ngrid
    278279         xmin = MIN(zstd(i),xmin)
    279280         xmax = MAX(zstd(i),xmax)
    … …  
    289290      ENDIF
    290291#ifdef NC_DOUBLE
    291       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zsig)
    292 #else
    293       ierr = NF_GET_VAR_REAL(nid, nvarid, zsig)
     292      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zsig)
     293#else
     294      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zsig)
    294295#endif
    295296      IF (ierr.NE.NF_NOERR) THEN
    … …  
    299300      xmin = 1.0E+20
    300301      xmax = -1.0E+20
    301       DO i = 1, ngridmx
     302      DO i = 1, ngrid
    302303         xmin = MIN(zsig(i),xmin)
    303304         xmax = MAX(zsig(i),xmax)
    … …  
    313314      ENDIF
    314315#ifdef NC_DOUBLE
    315       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zgam)
    316 #else
    317       ierr = NF_GET_VAR_REAL(nid, nvarid, zgam)
     316      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zgam)
     317#else
     318      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zgam)
    318319#endif
    319320      IF (ierr.NE.NF_NOERR) THEN
    … …  
    323324      xmin = 1.0E+20
    324325      xmax = -1.0E+20
    325       DO i = 1, ngridmx
     326      DO i = 1, ngrid
    326327         xmin = MIN(zgam(i),xmin)
    327328         xmax = MAX(zgam(i),xmax)
    … …  
    337338      ENDIF
    338339#ifdef NC_DOUBLE
    339       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, zthe)
    340 #else
    341       ierr = NF_GET_VAR_REAL(nid, nvarid, zthe)
     340      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zthe)
     341#else
     342      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zthe)
    342343#endif
    343344      IF (ierr.NE.NF_NOERR) THEN
    … …  
    347348      xmin = 1.0E+20
    348349      xmax = -1.0E+20
    349       DO i = 1, ngridmx
     350      DO i = 1, ngrid
    350351         xmin = MIN(zthe(i),xmin)
    351352         xmax = MAX(zthe(i),xmax)
    … …  
    413414         PRINT*, '          J ignore donc les autres temperatures TS**'
    414415#ifdef NC_DOUBLE
    415          ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsurf(1,1))
    416 #else
    417          ierr = NF_GET_VAR_REAL(nid, nvarid, tsurf(1,1))
     416         ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, tsurf)
     417#else
     418         ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, tsurf)
    418419#endif
    419420         IF (ierr.NE.NF_NOERR) THEN
    … …  
    428429         IF (nbsrf >= 2) THEN
    429430            DO nsrf = 2, nbsrf
    430                DO i = 1, ngridmx
     431               DO i = 1, ngrid
    431432                  tsurf(i,nsrf) = tsurf(i,1)
    432433               ENDDO
    … …  
    448449!               PRINT*, "phyetat0: Le champ <tsoil> est absent"
    449450!               PRINT*, "          Il prend donc la valeur de surface"
    450 !               DO i=1, ngridmx
     451!               DO i=1, ngrid
    451452!                  tsoil(i,isoil,nsrf)=tsurf(i,nsrf)
    452453!               ENDDO
    … …  
    478479      ENDIF
    479480#ifdef NC_DOUBLE
    480       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, emis)
    481 #else
    482       ierr = NF_GET_VAR_REAL(nid, nvarid, emis)
     481      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta,ngrid, emis)
     482#else
     483      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, emis)
    483484#endif
    484485      IF (ierr.NE.NF_NOERR) THEN
    … …  
    501502!         CALL abort
    502503      ENDIF
    503 #ifdef NC_DOUBLE
    504       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, cloudfrac)
    505 #else
    506       ierr = NF_GET_VAR_REAL(nid, nvarid, cloudfrac)
     504      sta2d = (/sta,1/)
     505      siz2d = (/ngrid, nlayermx/)
     506#ifdef NC_DOUBLE
     507      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,cloudfrac)
     508#else
     509      ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,cloudfrac)
    507510#endif
    508511      IF (ierr.NE.NF_NOERR) THEN
    … …  
    528531      ENDIF
    529532#ifdef NC_DOUBLE
    530       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, totcloudfrac)
    531 #else
    532       ierr = NF_GET_VAR_REAL(nid, nvarid, totcloudfrac)
     533      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, totcloudfrac)
     534#else
     535      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, totcloudfrac)
    533536#endif
    534537      IF (ierr.NE.NF_NOERR) THEN
    … …  
    553556      ENDIF
    554557#ifdef NC_DOUBLE
    555       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, hice)
    556 #else
    557       ierr = NF_GET_VAR_REAL(nid, nvarid, hice)
     558      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, hice)
     559#else
     560      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, hice)
    558561#endif
    559562      IF (ierr.NE.NF_NOERR) THEN
    … …  
    576579         CALL abort
    577580      ENDIF
    578 #ifdef NC_DOUBLE
    579       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, q2)
    580 #else
    581       ierr = NF_GET_VAR_REAL(nid, nvarid, q2)
     581      sta2d = (/sta,1/)
     582      siz2d = (/ngrid, llm+1/)
     583#ifdef NC_DOUBLE
     584      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,q2)
     585#else
     586      ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,q2)
    582587#endif
    583588      IF (ierr.NE.NF_NOERR) THEN
    … …  
    604609             write(txt,'(a5,i2.2)')'qsurf',iq
    605610           ELSE
    606              txt=tnom(iq)
     611             txt=noms(iq)
    607612!             if (txt.eq."h2o_vap") then
    608613               ! There is no surface tracer for h2o_vap;
    … …  
    618623     &                  ' not found in file'
    619624             write(*,*) trim(txt), ' set to 0'
    620              do ig=1,ngridmx
     625             do ig=1,ngrid
    621626               qsurf(ig,iq)=0.
    622627             end do
    623628             nqold=min(iq-1,nqold)
    624629           ELSE
    625 #ifdef NC_DOUBLE
    626              ierr = NF_GET_VAR_DOUBLE(nid, nvarid,qsurf(1,iq))
    627 #else
    628              ierr = NF_GET_VAR_REAL(nid, nvarid,qsurf(1,iq))
     630           toto(:) = 0.
     631#ifdef NC_DOUBLE
     632           ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,toto)
     633#else
     634           ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,toto)
    629635#endif
    630636             IF (ierr.NE.NF_NOERR) THEN
    … …  
    632638               CALL abort
    633639             ENDIF
     640             qsurf(1:ngrid,iq) = toto(1:ngrid)
    634641           ENDIF
    635642           xmin = 1.0E+20
    636643           xmax = -1.0E+20
    637            xmin = MINVAL(qsurf(1:ngridmx,iq))
    638            xmax = MAXVAL(qsurf(1:ngridmx,iq))
     644           xmin = MINVAL(qsurf(1:ngrid,iq))
     645           xmax = MAXVAL(qsurf(1:ngrid,iq))
    639646           PRINT*,'tracer on surface <',trim(txt),'>:',xmin,xmax
    640647         ENDDO
    … …  
    642649c        case when new tracer are added in addition to old ones
    643650             write(*,*)'qsurf 1 to ', nqold,'were already present'
    644              write(*,*)'qsurf ', nqold+1,' to ', nqmx,'are new'
     651             write(*,*)'qsurf ', nqold+1,' to ', nq,'are new'
    645652             write(*,*)' and initialized to zero'
    646              qsurf(:,nqold+1:nqmx)=0.0
     653             qsurf(:,nqold+1:nq)=0.0
    647654!            yes=' '
    648655!            do while ((yes.ne.'y').and.(yes.ne.'n'))
    649656!             write(*,*) 'Would you like to reindex qsurf # 1 ->',nqold
    650 !             write(*,*) 'to #',nqmx-nqold+1,'->', nqmx,'   (y or n) ?'
     657!             write(*,*) 'to #',nq-nqold+1,'->', nq,'   (y or n) ?'
    651658!             read(*,fmt='(a)') yes
    652659!            end do
    653660!            if (yes.eq.'y') then
    654661!              write(*,*) 'OK, let s reindex qsurf'
    655 !                 do ig=1,ngridmx
    656 !                    do iq=nqmx,nqmx-nqold+1,-1
    657 !                       qsurf(ig,iq)=qsurf(ig,iq-nqmx+nqold)
     662!                 do ig=1,ngrid
     663!                    do iq=nq,nq-nqold+1,-1
     664!                       qsurf(ig,iq)=qsurf(ig,iq-nq+nqold)
    658665!                    end do
    659 !                    do iq=nqmx-nqold,1,-1
     666!                    do iq=nq-nqold,1,-1
    660667!                       qsurf(ig,iq)= 0.
    661668!                    end do
    … …  
    668675! as well as thermal inertia and volumetric heat capacity
    669676
    670       call soil_settings(nid,ngridmx,nsoil,tsurf,tsoil)
     677      PRINT*,'SOIL INIT'
     678      call soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
    671679c
    672680c Fermer le fichier:
    … …  
    674682      ierr = NF_CLOSE(nid)
    675683c
     684
    676685      RETURN
    677686      END

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

    r726 r787  
    1       subroutine physdem1(filename,lonfi,latfi,nsoil,nq,
     1      subroutine physdem1(ngrid,filename,lonfi,latfi,nsoil,nq,
    22     .                   phystep,day_ini,
    33     .                   time,tsurf,tsoil,emis,q2,qsurf,
    44     .                   airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe,
    5      .                   cloudfrac,totcloudfrac,hice)
    6 
     5     .                   cloudfrac,totcloudfrac,hice,nametrac)
     6
     7      use surfdat_h
     8      use comsoil_h
     9      USE comgeomfi_h
    710 
    811      implicit none
    … …  
    3639#include "netcdf.inc"
    3740#include "dimphys.h"
    38 #include"advtrac.h"
    3941#include"callkeys.h"
     42
     43      INTEGER :: ngrid
    4044
    4145      INTEGER nid,iq
    … …  
    5054
    5155      REAL phystep,time
    52       REAL latfi(ngridmx), lonfi(ngridmx)
    53 !      REAL champhys(ngridmx)
    54       REAL tsurf(ngridmx)
     56      REAL latfi(ngrid), lonfi(ngrid)
     57!      REAL champhys(ngrid)
     58      REAL tsurf(ngrid)
    5559      INTEGER length
    5660      PARAMETER (length=100)
    … …  
    5862
    5963!     added by BC
    60       REAL hice(ngridmx),cloudfrac(ngridmx,nlayermx)
    61       REAL totcloudfrac(ngridmx)
     64      REAL hice(ngrid),cloudfrac(ngrid,nlayermx)
     65      REAL totcloudfrac(ngrid)
     66
     67!     AS: name of tracers added as an argument to avoid using nqmx in commons (i.e. adv trac)
     68!       nota: physdem1 is used both in physiq and newstart.
     69      character*20 nametrac(nq)   ! name of the tracer
    6270
    6371
    … …  
    6876#include "clesph0.h"
    6977#include "fxyprim.h"
    70 #include "comgeomfi.h"
    71 #include "surfdat.h"
    72 #include "comsoil.h"
    7378#include "planete.h"
    7479#include "comcstfi.h"
    7580
    76       real tsoil(ngridmx,nsoil),emis(ngridmx)
    77       real q2(ngridmx, llm+1),qsurf(ngridmx,nq)
    78       real airefi(ngridmx)
    79       real alb(ngridmx),ith(ngridmx,nsoil)
    80       real pzmea(ngridmx),pzstd(ngridmx)
    81       real pzsig(ngridmx),pzgam(ngridmx),pzthe(ngridmx)
     81      real tsoil(ngrid,nsoil),emis(ngrid)
     82      real q2(ngrid, llm+1),qsurf(ngrid,nq)
     83      real airefi(ngrid)
     84      real alb(ngrid),ith(ngrid,nsoil)
     85      real pzmea(ngrid),pzstd(ngrid)
     86      real pzsig(ngrid),pzgam(ngrid),pzthe(ngrid)
    8287      integer ig
    8388
    … …  
    9297
    9398      ! copy airefi(:) to area(:)
    94       CALL SCOPY(ngridmx,airefi,1,area,1)
     99      CALL SCOPY(ngrid,airefi,1,area,1)
    95100      ! note: area() is defined in comgeomfi.h
    96101
    97       DO ig=1,ngridmx
     102      DO ig=1,ngrid
    98103         albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat.h
    99104         zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat.h
    … …  
    127132      endif
    128133c
    129       ierr = NF_DEF_DIM (nid,"physical_points",ngridmx,idim2)
     134      ierr = NF_DEF_DIM (nid,"physical_points",ngrid,idim2)
    130135      if (ierr.ne.NF_NOERR) then
    131136        WRITE(6,*)'physdem1: Problem defining physical_points dimension'
    … …  
    169174      ENDDO
    170175
    171       write(*,*) "physdem1: ngridmx: ",ngridmx
     176      write(*,*) "physdem1: ngrid: ",ngrid
    172177
    173178ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    … …  
    175180ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    176181c Informations on the physics grid
    177       tab_cntrl(1) = float(ngridmx)  ! number of nodes on physics grid
     182      tab_cntrl(1) = float(ngrid)  ! number of nodes on physics grid
    178183      tab_cntrl(2) = float(nlayermx) ! number of atmospheric layers
    179184      tab_cntrl(3) = day_ini + int(time)         ! initial day
    … …  
    618623!                    qsurf02, ...)
    619624      count=0
    620       do iq=1,nqmx
     625      do iq=1,nq
    621626        txt= " "
    622627        write(txt,'(a1,i2.2)')'q',iq
    623         if (txt.ne.tnom(iq)) then ! use tracer names stored in dynamics
     628        if (txt.ne.nametrac(iq)) then ! use tracer names stored in dynamics
    624629          ! did not find old tracer name
    625630          exit ! might as well stop here
    … …  
    629634        endif
    630635      enddo
    631       if (count.eq.nqmx) then
     636      if (count.eq.nq) then
    632637        write(*,*) "physdem1:tracers seem to follow old naming ",
    633638     &             "convention (qsurf01,qsurf02,...)"
    … …  
    638643        !oldtracernames=.true.
    639644        ! Moreover, if computing water cycle with ice, move surface ice
    640         ! back to qsurf(nqmx)
     645        ! back to qsurf(nq)
    641646        !IF (iceparty) THEN
    642         !  write(*,*)'physdem1: moving surface water ice to index ',nqmx
    643         !  qsurf(1:ngridmx,nqmx)=qsurf(1:ngridmx,nqmx-1)
    644         !  qsurf(1:ngridmx,nqmx-1)=0
     647        !  write(*,*)'physdem1: moving surface water ice to index ',nq
     648        !  qsurf(1:ngrid,nq)=qsurf(1:ngrid,nq-1)
     649        !  qsurf(1:ngrid,nq-1)=0
    645650        !ENDIF
    646       endif ! of if (count.eq.nqmx)
     651      endif ! of if (count.eq.nq)
    647652
    648653      IF(nq.GE.1) THEN
    … …  
    650655        if (.not.oldtracernames) then
    651656         do iq=1,nq
    652            if (tnom(iq).eq."h2o_vap") then
     657           if (nametrac(iq).eq."h2o_vap") then
    653658             i_h2o_vap=iq
    654659           endif
    655            if (tnom(iq).eq."h2o_ice") then
     660           if (nametrac(iq).eq."h2o_ice") then
    656661             i_h2o_ice=iq
    657662           endif
    … …  
    671676             write(txt,'(a5,i2.2)')'qsurf',iq
    672677           ELSE
    673              txt=tnom(iq)
     678             txt=nametrac(iq)
    674679
    675680

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

    r775 r787  
    11      subroutine physiq(ngrid,nlayer,nq,   &
     2                  nametrac,                &
    23                  firstcall,lastcall,      &
    34                  pday,ptime,ptimestep,    &
    … …  
    1314      use radii_mod, only: h2o_reffrad, co2_reffrad
    1415      use aerosol_mod
     16      use surfdat_h
     17      use comdiurn_h
     18      use comsaison_h
     19      use comsoil_h
     20      USE comgeomfi_h
     21      USE tracer_h
     22
    1523      implicit none
    1624
    … …  
    103111!           New turbulent diffusion scheme: J. Leconte (2012)
    104112!           Loops converted to F90 matrix format: J. Leconte (2012)
     113!           No more ngridmx/nqmx, F90 commons and adaptation to parallel: A. Spiga (2012)
    105114!     
    106115!==================================================================
    … …  
    112121#include "dimensions.h"
    113122#include "dimphys.h"
    114 #include "comgeomfi.h"
    115 #include "surfdat.h"
    116 #include "comsoil.h"
    117 #include "comdiurn.h"
    118123#include "callkeys.h"
    119124#include "comcstfi.h"
    120125#include "planete.h"
    121 #include "comsaison.h"
    122126#include "control.h"
    123 #include "tracer.h"
    124 #include "watercap.h"
    125127#include "netcdf.inc"
    126128
    … …  
    132134      integer ngrid,nlayer,nq
    133135      real ptimestep
    134       real pplev(ngridmx,nlayer+1),pplay(ngridmx,nlayer)
    135       real pphi(ngridmx,nlayer)
    136       real pu(ngridmx,nlayer),pv(ngridmx,nlayer)
    137       real pt(ngridmx,nlayer),pq(ngridmx,nlayer,nq)
    138       real pw(ngridmx,nlayer)        ! pvervel transmitted by dyn3d
    139       real zh(ngridmx,nlayermx)      ! potential temperature (K)
     136      real pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
     137      real pphi(ngrid,nlayer)
     138      real pu(ngrid,nlayer),pv(ngrid,nlayer)
     139      real pt(ngrid,nlayer),pq(ngrid,nlayer,nq)
     140      real pw(ngrid,nlayer)        ! pvervel transmitted by dyn3d
     141      real zh(ngrid,nlayermx)      ! potential temperature (K)
    140142      logical firstcall,lastcall
    141143
    … …  
    144146      logical tracerdyn
    145147
     148      character*20 nametrac(nq)   ! name of the tracer from dynamics
     149
    146150!   outputs:
    147151!   --------
    148152!     physical tendencies
    149       real pdu(ngridmx,nlayer),pdv(ngridmx,nlayer)
    150       real pdt(ngridmx,nlayer),pdq(ngridmx,nlayer,nq)
    151       real pdpsrf(ngridmx) ! surface pressure tendency
     153      real pdu(ngrid,nlayer),pdv(ngrid,nlayer)
     154      real pdt(ngrid,nlayer),pdq(ngrid,nlayer,nq)
     155      real pdpsrf(ngrid) ! surface pressure tendency
    152156
    153157
    … …  
    157161!     "longrefvis" set in dimradmars.h , for one of the "naerkind"  kind of
    158162!      aerosol optical properties:
    159 !      real aerosol(ngridmx,nlayermx,naerkind)
     163!      real aerosol(ngrid,nlayermx,naerkind)
    160164!     this is now internal to callcorrk and hence no longer needed here
    161165
    162166      integer day_ini                ! Initial date of the run (sol since Ls=0)
    163167      integer icount                 ! counter of calls to physiq during the run.
    164       real tsurf(ngridmx)            ! Surface temperature (K)
    165       real tsoil(ngridmx,nsoilmx)    ! sub-surface temperatures (K)
    166       real albedo(ngridmx)           ! Surface albedo
    167 
    168       real albedo0(ngridmx)          ! Surface albedo
    169       integer rnat(ngridmx)          ! added by BC
    170       save rnat
    171 
    172       real emis(ngridmx)             ! Thermal IR surface emissivity
    173       real dtrad(ngridmx,nlayermx)   ! Net atm. radiative heating rate (K.s-1)
    174       real fluxrad_sky(ngridmx)      ! rad. flux from sky absorbed by surface (W.m-2)
    175       real fluxrad(ngridmx)          ! Net radiative surface flux (W.m-2)
    176       real capcal(ngridmx)           ! surface heat capacity (J m-2 K-1)
    177       real fluxgrd(ngridmx)          ! surface conduction flux (W.m-2)
    178       real qsurf(ngridmx,nqmx)       ! tracer on surface (e.g. kg.m-2)
    179       real q2(ngridmx,nlayermx+1)    ! Turbulent Kinetic Energy
     168      real, dimension(:),allocatable,save ::  tsurf  ! Surface temperature (K)
     169      real, dimension(:,:),allocatable,save ::  tsoil  ! sub-surface temperatures (K)
     170      real, dimension(:),allocatable,save :: albedo ! Surface albedo
     171
     172      real,dimension(:),allocatable,save :: albedo0 ! Surface albedo
     173      integer,dimension(:),allocatable,save :: rnat ! added by BC
     174
     175      real,dimension(:),allocatable,save :: emis ! Thermal IR surface emissivity
     176      real,dimension(:,:),allocatable,save :: dtrad ! Net atm. radiative heating rate (K.s-1)
     177      real,dimension(:),allocatable,save :: fluxrad_sky ! rad. flux from sky absorbed by surface (W.m-2)
     178      real,dimension(:),allocatable,save :: fluxrad ! Net radiative surface flux (W.m-2)
     179      real,dimension(:),allocatable,save :: capcal ! surface heat capacity (J m-2 K-1)
     180      real,dimension(:),allocatable,save :: fluxgrd ! surface conduction flux (W.m-2)
     181      real,dimension(:,:),allocatable,save :: qsurf ! tracer on surface (e.g. kg.m-2)
     182      real,dimension(:,:),allocatable,save :: q2 ! Turbulent Kinetic Energy
    180183
    181184      save day_ini, icount
    182       save tsurf,tsoil
    183       save albedo0,albedo,emis,q2
    184       save capcal,fluxgrd,dtrad,fluxrad,fluxrad_sky,qsurf
    185185
    186186! Local variables :
    … …  
    189189!     aerosol (dust or ice) extinction optical depth  at reference wavelength
    190190!     for the "naerkind" optically active aerosols:
    191       real aerosol(ngridmx,nlayermx,naerkind)
     191      real aerosol(ngrid,nlayermx,naerkind)
    192192
    193193      character*80 fichier
    194194      integer l,ig,ierr,iq,iaer
    195195
    196       real fluxsurf_lw(ngridmx)      ! incident LW (IR) surface flux (W.m-2)
    197       real fluxsurf_sw(ngridmx)      ! incident SW (stellar) surface flux (W.m-2)
    198       real fluxtop_lw(ngridmx)       ! Outgoing LW (IR) flux to space (W.m-2)
    199       real fluxabs_sw(ngridmx)       ! Absorbed SW (stellar) flux (W.m-2)
    200 
    201       real fluxtop_dn(ngridmx)
    202       real fluxdyn(ngridmx)          ! horizontal heat transport by dynamics     
    203       real OLR_nu(ngridmx,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
    204       real OSR_nu(ngridmx,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
    205  
    206       real,save :: sensibFlux(ngridmx)   ! turbulent flux given by the atm to the surface
    207  
    208       save fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn,fluxdyn,OLR_nu,OSR_nu
    209 
     196      !!! this is saved for diagnostic
     197      real,dimension(:),allocatable,save :: fluxsurf_lw ! incident LW (IR) surface flux (W.m-2)
     198      real,dimension(:),allocatable,save :: fluxsurf_sw ! incident SW (stellar) surface flux (W.m-2)
     199      real,dimension(:),allocatable,save :: fluxtop_lw ! Outgoing LW (IR) flux to space (W.m-2)
     200      real,dimension(:),allocatable,save :: fluxabs_sw ! Absorbed SW (stellar) flux (W.m-2)
     201      real,dimension(:),allocatable,save :: fluxtop_dn
     202      real,dimension(:),allocatable,save :: fluxdyn ! horizontal heat transport by dynamics 
     203      real,dimension(:,:),allocatable,save :: OLR_nu ! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1)
     204      real,dimension(:,:),allocatable,save :: OSR_nu ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1)
     205      real,dimension(:,:),allocatable,save :: zdtlw ! (K/s)
     206      real,dimension(:,:),allocatable,save :: zdtsw ! (K/s)
     207      real,dimension(:),allocatable,save :: sensibFlux ! turbulent flux given by the atm to the surface
    210208
    211209      real zls                       ! solar longitude (rad)
    212210      real zday                      ! date (time since Ls=0, in martian days)
    213       real zzlay(ngridmx,nlayermx)   ! altitude at the middle of the layers
    214       real zzlev(ngridmx,nlayermx+1) ! altitude at layer boundaries
     211      real zzlay(ngrid,nlayermx)   ! altitude at the middle of the layers
     212      real zzlev(ngrid,nlayermx+1) ! altitude at layer boundaries
    215213      real latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
    216214
    217215!     Tendencies due to various processes:
    218       real dqsurf(ngridmx,nqmx)
    219       real,save :: zdtlw(ngridmx,nlayermx)                    ! (K/s)
    220       real,save :: zdtsw(ngridmx,nlayermx)                    ! (K/s)
    221 
    222       real cldtlw(ngridmx,nlayermx)                           ! (K/s) LW heating rate for clear areas
    223       real cldtsw(ngridmx,nlayermx)                           ! (K/s) SW heating rate for clear areas
    224       real zdtsurf(ngridmx)                                   ! (K/s)
    225       real dtlscale(ngridmx,nlayermx)
    226       real zdvdif(ngridmx,nlayermx),zdudif(ngridmx,nlayermx)  ! (m.s-2)
    227       real zdhdif(ngridmx,nlayermx), zdtsdif(ngridmx)         ! (K/s)
    228       real zdtdif(ngridmx,nlayermx)                           ! (K/s)
    229       real zdvadj(ngridmx,nlayermx),zduadj(ngridmx,nlayermx)  ! (m.s-2)
    230       real zdhadj(ngridmx,nlayermx)                           ! (K/s)
    231       real zdtgw(ngridmx,nlayermx)                            ! (K/s)
    232       real zdtmr(ngridmx,nlayermx)
    233       real zdugw(ngridmx,nlayermx),zdvgw(ngridmx,nlayermx)    ! (m.s-2)
    234       real zdtc(ngridmx,nlayermx),zdtsurfc(ngridmx)
    235       real zdvc(ngridmx,nlayermx),zduc(ngridmx,nlayermx)
    236       real zdumr(ngridmx,nlayermx),zdvmr(ngridmx,nlayermx)
    237       real zdtsurfmr(ngridmx)
     216      real dqsurf(ngrid,nq)
     217      real cldtlw(ngrid,nlayermx)                           ! (K/s) LW heating rate for clear areas
     218      real cldtsw(ngrid,nlayermx)                           ! (K/s) SW heating rate for clear areas
     219      real zdtsurf(ngrid)                                   ! (K/s)
     220      real dtlscale(ngrid,nlayermx)
     221      real zdvdif(ngrid,nlayermx),zdudif(ngrid,nlayermx)  ! (m.s-2)
     222      real zdhdif(ngrid,nlayermx), zdtsdif(ngrid)         ! (K/s)
     223      real zdtdif(ngrid,nlayermx)                             ! (K/s)
     224      real zdvadj(ngrid,nlayermx),zduadj(ngrid,nlayermx)  ! (m.s-2)
     225      real zdhadj(ngrid,nlayermx)                           ! (K/s)
     226      real zdtgw(ngrid,nlayermx)                            ! (K/s)
     227      real zdtmr(ngrid,nlayermx)
     228      real zdugw(ngrid,nlayermx),zdvgw(ngrid,nlayermx)    ! (m.s-2)
     229      real zdtc(ngrid,nlayermx),zdtsurfc(ngrid)
     230      real zdvc(ngrid,nlayermx),zduc(ngrid,nlayermx)
     231      real zdumr(ngrid,nlayermx),zdvmr(ngrid,nlayermx)
     232      real zdtsurfmr(ngrid)
    238233     
    239       real zdmassmr(ngridmx,nlayermx),zdpsrfmr(ngridmx)
    240 
    241       real zdqdif(ngridmx,nlayermx,nqmx), zdqsdif(ngridmx,nqmx)
    242       real zdqevap(ngridmx,nlayermx)
    243       real zdqsed(ngridmx,nlayermx,nqmx), zdqssed(ngridmx,nqmx)
    244       real zdqdev(ngridmx,nlayermx,nqmx), zdqsdev(ngridmx,nqmx)
    245       real zdqadj(ngridmx,nlayermx,nqmx)
    246       real zdqc(ngridmx,nlayermx,nqmx)
    247       real zdqmr(ngridmx,nlayermx,nqmx),zdqsurfmr(ngridmx,nqmx)
    248       real zdqlscale(ngridmx,nlayermx,nqmx)
    249       real zdqslscale(ngridmx,nqmx)
    250       real zdqchim(ngridmx,nlayermx,nqmx)
    251       real zdqschim(ngridmx,nqmx)
    252 
    253       real zdteuv(ngridmx,nlayermx)    ! (K/s)
    254       real zdtconduc(ngridmx,nlayermx) ! (K/s)
    255       real zdumolvis(ngridmx,nlayermx)
    256       real zdvmolvis(ngridmx,nlayermx)
    257       real zdqmoldiff(ngridmx,nlayermx,nqmx)
     234      real zdmassmr(ngrid,nlayermx),zdpsrfmr(ngrid)
     235
     236      real zdqdif(ngrid,nlayermx,nq), zdqsdif(ngrid,nq)
     237      real zdqevap(ngrid,nlayermx)
     238      real zdqsed(ngrid,nlayermx,nq), zdqssed(ngrid,nq)
     239      real zdqdev(ngrid,nlayermx,nq), zdqsdev(ngrid,nq)
     240      real zdqadj(ngrid,nlayermx,nq)
     241      real zdqc(ngrid,nlayermx,nq)
     242      real zdqmr(ngrid,nlayermx,nq),zdqsurfmr(ngrid,nq)
     243      real zdqlscale(ngrid,nlayermx,nq)
     244      real zdqslscale(ngrid,nq)
     245      real zdqchim(ngrid,nlayermx,nq)
     246      real zdqschim(ngrid,nq)
     247
     248      real zdteuv(ngrid,nlayermx)    ! (K/s)
     249      real zdtconduc(ngrid,nlayermx) ! (K/s)
     250      real zdumolvis(ngrid,nlayermx)
     251      real zdvmolvis(ngrid,nlayermx)
     252      real zdqmoldiff(ngrid,nlayermx,nq)
    258253
    259254!     Local variables for local calculations:
    260       real zflubid(ngridmx)
    261       real zplanck(ngridmx),zpopsk(ngridmx,nlayermx)
    262       real zdum1(ngridmx,nlayermx)
    263       real zdum2(ngridmx,nlayermx)
     255      real zflubid(ngrid)
     256      real zplanck(ngrid),zpopsk(ngrid,nlayermx)
     257      real zdum1(ngrid,nlayermx)
     258      real zdum2(ngrid,nlayermx)
    264259      real ztim1,ztim2,ztim3, z1,z2
    265260      real ztime_fin
    266       real zdh(ngridmx,nlayermx)
     261      real zdh(ngrid,nlayermx)
    267262      integer length
    268263      parameter (length=100)
    … …  
    270265! local variables only used for diagnostics (output in file "diagfi" or "stats")
    271266! ------------------------------------------------------------------------------
    272       real ps(ngridmx), zt(ngridmx,nlayermx)
    273       real zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
    274       real zq(ngridmx,nlayermx,nqmx)
     267      real ps(ngrid), zt(ngrid,nlayermx)
     268      real zu(ngrid,nlayermx),zv(ngrid,nlayermx)
     269      real zq(ngrid,nlayermx,nq)
    275270      character*2 str2
    276271      character*5 str5
    277       real zdtadj(ngridmx,nlayermx)
    278       real zdtdyn(ngridmx,nlayermx),ztprevious(ngridmx,nlayermx)
    279       save ztprevious
    280       real reff(ngridmx,nlayermx) ! effective dust radius (used if doubleq=T)
     272      real zdtadj(ngrid,nlayermx)
     273      real zdtdyn(ngrid,nlayermx)
     274      real,allocatable,dimension(:,:),save :: ztprevious
     275      real reff(ngrid,nlayermx) ! effective dust radius (used if doubleq=T)
    281276      real qtot1,qtot2            ! total aerosol mass
    282277      integer igmin, lmin
    … …  
    285280      real zplev(ngrid,nlayermx+1),zplay(ngrid,nlayermx)
    286281      real zstress(ngrid), cd
    287       real hco2(nqmx), tmean, zlocal(nlayermx)
    288       real vmr(ngridmx,nlayermx) ! volume mixing ratio
     282      real hco2(nq), tmean, zlocal(nlayermx)
     283      real vmr(ngrid,nlayermx) ! volume mixing ratio
    289284
    290285      real time_phys
    291286
    292287!     reinstated by RW for diagnostic
    293       real tau_col(ngridmx)
    294       save tau_col
     288      real,allocatable,dimension(:),save :: tau_col
    295289     
    296290!     included by RW to reduce insanity of code
    … …  
    298292
    299293!     included by RW to compute tracer column densities
    300       real qcol(ngridmx,nqmx)
     294      real qcol(ngrid,nq)
    301295
    302296!     included by RW for H2O precipitation
    303       real zdtrain(ngridmx,nlayermx)
    304       real zdqrain(ngridmx,nlayermx,nqmx)
    305       real zdqsrain(ngridmx)
    306       real zdqssnow(ngridmx)
    307       real rainout(ngridmx)
     297      real zdtrain(ngrid,nlayermx)
     298      real zdqrain(ngrid,nlayermx,nq)
     299      real zdqsrain(ngrid)
     300      real zdqssnow(ngrid)
     301      real rainout(ngrid)
    308302
    309303!     included by RW for H2O Manabe scheme
    310       real dtmoist(ngridmx,nlayermx)
    311       real dqmoist(ngridmx,nlayermx,nqmx)
    312 
    313       real qvap(ngridmx,nlayermx)
    314       real dqvaplscale(ngridmx,nlayermx)
    315       real dqcldlscale(ngridmx,nlayermx)
    316       real rneb_man(ngridmx,nlayermx)
    317       real rneb_lsc(ngridmx,nlayermx)
     304      real dtmoist(ngrid,nlayermx)
     305      real dqmoist(ngrid,nlayermx,nq)
     306
     307      real qvap(ngrid,nlayermx)
     308      real dqvaplscale(ngrid,nlayermx)
     309      real dqcldlscale(ngrid,nlayermx)
     310      real rneb_man(ngrid,nlayermx)
     311      real rneb_lsc(ngrid,nlayermx)
    318312
    319313!     included by RW to account for surface cooling by evaporation
    320       real dtsurfh2olat(ngridmx)
     314      real dtsurfh2olat(ngrid)
    321315
    322316
    323317!     to test energy conservation (RW+JL)
    324       real mass(ngridmx,nlayermx),massarea(ngridmx,nlayermx)
     318      real mass(ngrid,nlayermx),massarea(ngrid,nlayermx)
    325319      real dEtot, dEtots, AtmToSurf_TurbFlux
    326320      real dEtotSW, dEtotsSW, dEtotLW, dEtotsLW
    327       real dEzRadsw(ngridmx,nlayermx),dEzRadlw(ngridmx,nlayermx),dEzdiff(ngridmx,nlayermx)
    328       real dEdiffs(ngridmx),dEdiff(ngridmx)
    329       real madjdE(ngridmx), lscaledE(ngridmx)
     321      real dEzRadsw(ngrid,nlayermx),dEzRadlw(ngrid,nlayermx),dEzdiff(ngrid,nlayermx)
     322      real dEdiffs(ngrid),dEdiff(ngrid)
     323      real madjdE(ngrid), lscaledE(ngrid)
    330324!JL12 conservation test for mean flow kinetic energy has been disabled temporarily
    331325     
    … …  
    333327
    334328!     included by BC for evaporation     
    335       real qevap(ngridmx,nlayermx,nqmx)
    336       real tevap(ngridmx,nlayermx)
    337       real dqevap1(ngridmx,nlayermx)
    338       real dtevap1(ngridmx,nlayermx)
     329      real qevap(ngrid,nlayermx,nq)
     330      real tevap(ngrid,nlayermx)
     331      real dqevap1(ngrid,nlayermx)
     332      real dtevap1(ngrid,nlayermx)
    339333
    340334!     included by BC for hydrology
    341       real hice(ngridmx)
     335      real hice(ngrid)
    342336
    343337!     included by RW to test water conservation (by routine)
    … …  
    349343
    350344!     included by RW for RH diagnostic
    351       real qsat(ngridmx,nlayermx), RH(ngridmx,nlayermx), H2Omaxcol(ngridmx),psat_tmp
     345      real qsat(ngrid,nlayermx), RH(ngrid,nlayermx), H2Omaxcol(ngrid),psat_tmp
    352346
    353347!     included by RW for hydrology
    354       real dqs_hyd(ngridmx,nqmx)
    355       real zdtsurf_hyd(ngridmx)
     348      real dqs_hyd(ngrid,nq)
     349      real zdtsurf_hyd(ngrid)
    356350
    357351!     included by RW for water cycle conservation tests
    … …  
    360354!     included by BC for double radiative transfer call
    361355      logical clearsky
    362       real zdtsw1(ngridmx,nlayermx)
    363       real zdtlw1(ngridmx,nlayermx)
    364       real fluxsurf_lw1(ngridmx)     
    365       real fluxsurf_sw1(ngridmx) 
    366       real fluxtop_lw1(ngridmx)   
    367       real fluxabs_sw1(ngridmx) 
     356      real zdtsw1(ngrid,nlayermx)
     357      real zdtlw1(ngrid,nlayermx)
     358      real fluxsurf_lw1(ngrid)     
     359      real fluxsurf_sw1(ngrid) 
     360      real fluxtop_lw1(ngrid)   
     361      real fluxabs_sw1(ngrid) 
    368362      real tau_col1(ngrid)
    369363      real OLR_nu1(ngrid,L_NSPECTI)
    … …  
    372366
    373367!     included by BC for cloud fraction computations
    374       real,save :: cloudfrac(ngridmx,nlayermx)
    375       real,save :: totcloudfrac(ngridmx)
     368      real,allocatable,dimension(:,:),save :: cloudfrac
     369      real,allocatable,dimension(:),save :: totcloudfrac
    376370
    377371!     included by RW for vdifc water conservation test
    378372      real nconsMAX
    379       real vdifcncons(ngridmx)
    380       real cadjncons(ngridmx)
    381 
    382 !      double precision qsurf_hist(ngridmx,nqmx)
    383       real qsurf_hist(ngridmx,nqmx)
    384       save qsurf_hist
     373      real vdifcncons(ngrid)
     374      real cadjncons(ngrid)
     375
     376!      double precision qsurf_hist(ngrid,nq)
     377      real,allocatable,dimension(:,:),save :: qsurf_hist
    385378
    386379!     included by RW for temp convadj conservation test
    … …  
    392385
    393386!     included by RW for runway greenhouse 1D study
    394       real muvar(ngridmx,nlayermx+1)
     387      real muvar(ngrid,nlayermx+1)
    395388
    396389!     included by RW for variable H2O particle sizes
    397       real,save :: reffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius (m)
    398       real,save :: nueffrad(ngridmx,nlayermx,naerkind) ! aerosol effective radius variance
    399       real :: reffh2oliq(ngridmx,nlayermx) ! liquid water particles effective radius (m)
    400       real :: reffh2oice(ngridmx,nlayermx) ! water ice particles effective radius (m)
    401       real reffH2O(ngridmx,nlayermx)
    402       real reffcol(ngridmx,naerkind)
     390      real,allocatable,dimension(:,:,:),save :: reffrad ! aerosol effective radius (m)
     391      real,allocatable,dimension(:,:,:),save :: nueffrad ! aerosol effective radius variance
     392      real :: nueffrad_dummy(ngrid,nlayermx,naerkind) !! AS. This is temporary. Check below why.
     393      real :: reffh2oliq(ngrid,nlayermx) ! liquid water particles effective radius (m)
     394      real :: reffh2oice(ngrid,nlayermx) ! water ice particles effective radius (m)
     395      real reffH2O(ngrid,nlayermx)
     396      real reffcol(ngrid,naerkind)
    403397
    404398!     included by RW for sourceevol
    405       real, save :: ice_initial(ngridmx)!, isoil
     399      real,allocatable,dimension(:),save :: ice_initial
    406400      real delta_ice,ice_tot
    407       real, save :: ice_min(ngridmx)
     401      real,allocatable,dimension(:),save :: ice_min
    408402     
    409403      integer num_run
    … …  
    411405     
    412406!=======================================================================
    413      
     407
     408      !!! ALLOCATE SAVED ARRAYS
     409      IF (.not. ALLOCATED(tsurf)) ALLOCATE(tsurf(ngrid))
     410      IF (.not. ALLOCATED(tsoil)) ALLOCATE(tsoil(ngrid,nsoilmx))   
     411      IF (.not. ALLOCATED(albedo)) ALLOCATE(albedo(ngrid))         
     412      IF (.not. ALLOCATED(albedo0)) ALLOCATE(albedo0(ngrid))         
     413      IF (.not. ALLOCATED(rnat)) ALLOCATE(rnat(ngrid))         
     414      IF (.not. ALLOCATED(emis)) ALLOCATE(emis(ngrid))   
     415      IF (.not. ALLOCATED(dtrad)) ALLOCATE(dtrad(ngrid,nlayermx))
     416      IF (.not. ALLOCATED(fluxrad_sky)) ALLOCATE(fluxrad_sky(ngrid))
     417      IF (.not. ALLOCATED(fluxrad)) ALLOCATE(fluxrad(ngrid))   
     418      IF (.not. ALLOCATED(capcal)) ALLOCATE(capcal(ngrid))   
     419      IF (.not. ALLOCATED(fluxgrd)) ALLOCATE(fluxgrd(ngrid)) 
     420      IF (.not. ALLOCATED(qsurf)) ALLOCATE(qsurf(ngrid,nq)) 
     421      IF (.not. ALLOCATED(q2)) ALLOCATE(q2(ngrid,nlayermx+1))
     422      IF (.not. ALLOCATED(ztprevious)) ALLOCATE(ztprevious(ngrid,nlayermx))
     423      IF (.not. ALLOCATED(cloudfrac)) ALLOCATE(cloudfrac(ngrid,nlayermx))
     424      IF (.not. ALLOCATED(totcloudfrac)) ALLOCATE(totcloudfrac(ngrid))
     425      IF (.not. ALLOCATED(qsurf_hist)) ALLOCATE(qsurf_hist(ngrid,nq))
     426      IF (.not. ALLOCATED(reffrad)) ALLOCATE(reffrad(ngrid,nlayermx,naerkind))
     427      IF (.not. ALLOCATED(nueffrad)) ALLOCATE(nueffrad(ngrid,nlayermx,naerkind))
     428      IF (.not. ALLOCATED(ice_initial)) ALLOCATE(ice_initial(ngrid))
     429      IF (.not. ALLOCATED(ice_min)) ALLOCATE(ice_min(ngrid))
     430
     431      IF (.not. ALLOCATED(fluxsurf_lw)) ALLOCATE(fluxsurf_lw(ngrid))
     432      IF (.not. ALLOCATED(fluxsurf_sw)) ALLOCATE(fluxsurf_sw(ngrid))
     433      IF (.not. ALLOCATED(fluxtop_lw)) ALLOCATE(fluxtop_lw(ngrid))
     434      IF (.not. ALLOCATED(fluxabs_sw)) ALLOCATE(fluxabs_sw(ngrid))
     435      IF (.not. ALLOCATED(fluxtop_dn)) ALLOCATE(fluxtop_dn(ngrid))
     436      IF (.not. ALLOCATED(fluxdyn)) ALLOCATE(fluxdyn(ngrid))
     437      IF (.not. ALLOCATED(OLR_nu)) ALLOCATE(OLR_nu(ngrid,L_NSPECTI))
     438      IF (.not. ALLOCATED(OSR_nu)) ALLOCATE(OSR_nu(ngrid,L_NSPECTV))
     439      IF (.not. ALLOCATED(sensibFlux)) ALLOCATE(sensibFlux(ngrid))
     440      IF (.not. ALLOCATED(zdtlw)) ALLOCATE(zdtlw(ngrid,nlayermx))
     441      IF (.not. ALLOCATED(zdtsw)) ALLOCATE(zdtsw(ngrid,nlayermx))
     442      IF (.not. ALLOCATED(tau_col)) ALLOCATE(tau_col(ngrid))
     443
     444!=======================================================================
    414445
    415446! 1. Initialisation
    … …  
    420451      if (firstcall) then
    421452
     453         print*,'FIRSTCALL INITIALIZATION'
     454
     455         !! this is defined in comsaison_h
     456         ALLOCATE(mu0(ngrid))
     457         ALLOCATE(fract(ngrid))
    422458
    423459!        variables set to 0
    … …  
    438474         tracerdyn=tracer
    439475         if (tracer) then
    440             call initracer()
     476            call initracer(ngrid,nq,nametrac)
    441477         endif ! end tracer
    442478
    … …  
    445481!        read startfi (initial parameters)
    446482!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    447          call phyetat0("startfi.nc",0,0,nsoilmx,nq,           &
     483         call phyetat0(ngrid,"startfi.nc",0,0,nsoilmx,nq,   &
    448484               day_ini,time_phys,tsurf,tsoil,emis,q2,qsurf,   &
    449485               cloudfrac,totcloudfrac,hice)
    … …  
    461497!        Initialize albedo and orbital calculation
    462498!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    463          call surfini(ngrid,qsurf,albedo0)
     499         call surfini(ngrid,nq,qsurf,albedo0)
    464500         call iniorbit(apoastr,periastr,year_day,peri_day,obliquit)
    465501
    … …  
    474510!        ~~~~~~~~~~~~~~~
    475511         if (callsoil) then
    476             call soil(ngrid,nsoilmx,firstcall,inertiedat,    &
     512            call soil(ngrid,nsoilmx,firstcall,lastcall,inertiedat, &
    477513                ptimestep,tsurf,tsoil,capcal,fluxgrd)
    478514         else
    … …  
    506542!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    507543         rnat(:)=1
    508          do ig=1,ngridmx
     544         do ig=1,ngrid
    509545!            if(iceball.or.oceanball.or.(inertiedat(ig,1).gt.1.E4))then
    510546            if(inertiedat(ig,1).gt.1.E4)then
    … …  
    562598                            ! need epsi for the wvp definitions in callcorrk.F
    563599
     600         print*,'end FIRSTCALL INITIALIZATION'
     601
    564602      endif        !  (end of "if firstcall")
    565603
    … …  
    568606! ---------------------------------------------------
    569607
    570       if (ngrid.NE.ngridmx) then
    571          print*,'STOP in PHYSIQ'
    572          print*,'Probleme de dimensions :'
    573          print*,'ngrid     = ',ngrid
    574          print*,'ngridmx   = ',ngridmx
    575          stop
    576       endif
     608      print*,'ALL INITIALIZATION'
    577609
    578610!     Initialize various variables
    579611!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    580612
    581       pdu(1:ngridmx,1:nlayermx) = 0.0
    582       pdv(1:ngridmx,1:nlayermx) = 0.0
     613      pdu(1:ngrid,1:nlayermx) = 0.0
     614      pdv(1:ngrid,1:nlayermx) = 0.0
    583615      if ( .not.nearco2cond ) then
    584          pdt(1:ngridmx,1:nlayermx) = 0.0
     616         pdt(1:ngrid,1:nlayermx) = 0.0
    585617      endif
    586       pdq(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
    587       pdpsrf(1:ngridmx)       = 0.0
    588       zflubid(1:ngridmx)      = 0.0
    589       zdtsurf(1:ngridmx)      = 0.0
    590       dqsurf(1:ngridmx,1:nqmx)= 0.0
     618      pdq(1:ngrid,1:nlayermx,1:nq) = 0.0
     619      pdpsrf(1:ngrid)       = 0.0
     620      zflubid(1:ngrid)      = 0.0
     621      zdtsurf(1:ngrid)      = 0.0
     622      dqsurf(1:ngrid,1:nq)= 0.0
    591623
    592624      zday=pday+ptime ! compute time, in sols (and fraction thereof)
    … …  
    603635!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    604636
    605       zzlay(1:ngridmx,1:nlayermx)=pphi(1:ngridmx,1:nlayermx)/g
    606 
    607       zzlev(1:ngridmx,1)=0.
    608       zzlev(1:ngridmx,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
     637      zzlay(1:ngrid,1:nlayermx)=pphi(1:ngrid,1:nlayermx)/g
     638
     639      zzlev(1:ngrid,1)=0.
     640      zzlev(1:ngrid,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
    609641
    610642      do l=2,nlayer
    … …  
    621653
    622654      do l=1,nlayer         
    623          do ig=1,ngrid 
     655         do ig=1,ngrid
    624656            zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
    625657            zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
    … …  
    628660         enddo
    629661      enddo
     662
     663      print*,'end ALL INITIALIZATION'
    630664
    631665!-----------------------------------------------------------------------
    … …  
    684718                      zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    685719                      fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,               &
    686                       reffrad,tau_col,cloudfrac,totcloudfrac,            &
     720                      reffrad,nueffrad,tau_col,cloudfrac,totcloudfrac,   &
    687721                      clearsky,firstcall,lastcall)     
    688722
    … …  
    698732                      zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, &
    699733                      fluxabs_sw1,fluxtop_dn,OLR_nu1,OSR_nu1,              &
    700                       reffrad,tau_col1,cloudfrac,totcloudfrac,             &
     734                      reffrad,nueffrad,tau_col1,cloudfrac,totcloudfrac,    &
    701735                      clearsky,firstcall,lastcall)
    702736                 clearsky = .false.  ! just in case.     
    … …  
    725759!             Radiative flux from the sky absorbed by the surface (W.m-2)
    726760              GSR=0.0
    727               fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)+fluxsurf_sw(1:ngridmx)*(1.-albedo(1:ngridmx))
     761              fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)+fluxsurf_sw(1:ngrid)*(1.-albedo(1:ngrid))
    728762
    729763              if(noradsurf)then ! no lower surface; SW flux just disappears
    730                  GSR = SUM(fluxsurf_sw(1:ngridmx)*area(1:ngridmx))/totarea
    731                  fluxrad_sky(1:ngridmx)=emis(1:ngridmx)*fluxsurf_lw(1:ngridmx)
     764                 GSR = SUM(fluxsurf_sw(1:ngrid)*area(1:ngrid))/totarea
     765                 fluxrad_sky(1:ngrid)=emis(1:ngrid)*fluxsurf_lw(1:ngrid)
    732766                 print*,'SW lost in deep atmosphere = ',GSR,' W m^-2'
    733767              endif
    734768
    735769!             Net atmospheric radiative heating rate (K.s-1)
    736               dtrad(1:ngridmx,1:nlayermx)=zdtsw(1:ngridmx,1:nlayermx)+zdtlw(1:ngridmx,1:nlayermx)
     770              dtrad(1:ngrid,1:nlayermx)=zdtsw(1:ngrid,1:nlayermx)+zdtlw(1:ngrid,1:nlayermx)
    737771
    738772           elseif(newtonian)then
    … …  
    740774!          b) Call Newtonian cooling scheme
    741775!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    742               call newtrelax(mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
    743 
    744               zdtsurf(1:ngridmx) = +(pt(1:ngridmx,1)-tsurf(1:ngridmx))/ptimestep
     776              call newtrelax(ngrid,mu0,sinlat,zpopsk,pt,pplay,pplev,dtrad,firstcall)
     777
     778              zdtsurf(1:ngrid) = +(pt(1:ngrid,1)-tsurf(1:ngrid))/ptimestep
    745779              ! e.g. surface becomes proxy for 1st atmospheric layer ?
    746780
    … …  
    749783!          c) Atmosphere has no radiative effect
    750784!          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    751               fluxtop_dn(1:ngridmx)  = fract(1:ngridmx)*mu0(1:ngridmx)*Fat1AU/dist_star**2
     785              fluxtop_dn(1:ngrid)  = fract(1:ngrid)*mu0(1:ngrid)*Fat1AU/dist_star**2
    752786              if(ngrid.eq.1)then ! / by 4 globally in 1D case...
    753787                 fluxtop_dn(1)  = fract(1)*Fat1AU/dist_star**2/2.0
    754788              endif
    755               fluxsurf_sw(1:ngridmx) = fluxtop_dn(1:ngridmx)
    756               fluxrad_sky(1:ngridmx) = fluxtop_dn(1:ngridmx)*(1.-albedo(1:ngridmx))
    757               fluxtop_lw(1:ngridmx)  = emis(1:ngridmx)*sigma*tsurf(1:ngridmx)**4
     789              fluxsurf_sw(1:ngrid) = fluxtop_dn(1:ngrid)
     790              fluxrad_sky(1:ngrid) = fluxtop_dn(1:ngrid)*(1.-albedo(1:ngrid))
     791              fluxtop_lw(1:ngrid)  = emis(1:ngrid)*sigma*tsurf(1:ngrid)**4
    758792              ! radiation skips the atmosphere entirely
    759793
    760794
    761               dtrad(1:ngridmx,1:nlayermx)=0.0
     795              dtrad(1:ngrid,1:nlayermx)=0.0
    762796              ! hence no atmospheric radiative heating
    763797
    … …  
    765799
    766800        endif ! of if(mod(icount-1,iradia).eq.0)
    767 
     801       
    768802
    769803!    Transformation of the radiative tendencies
    770804!    ------------------------------------------
    771805
    772         zplanck(1:ngridmx)=tsurf(1:ngridmx)*tsurf(1:ngridmx)
    773         zplanck(1:ngridmx)=emis(1:ngridmx)*sigma*zplanck(1:ngridmx)*zplanck(1:ngridmx)
    774         fluxrad(1:ngridmx)=fluxrad_sky(1:ngridmx)-zplanck(1:ngridmx)
    775         pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+dtrad(1:ngridmx,1:nlayermx)
     806        zplanck(1:ngrid)=tsurf(1:ngrid)*tsurf(1:ngrid)
     807        zplanck(1:ngrid)=emis(1:ngrid)*sigma*zplanck(1:ngrid)*zplanck(1:ngrid)
     808        fluxrad(1:ngrid)=fluxrad_sky(1:ngrid)-zplanck(1:ngrid)
     809        pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+dtrad(1:ngrid,1:nlayermx)
    776810
    777811!-------------------------
    … …  
    794828!-------------------------
    795829
     830      print*,'end RADIATIVE TRANSFER'
    796831      endif ! of if (callrad)
    797832
    … …  
    802837      if (calldifv) then
    803838     
    804          zflubid(1:ngridmx)=fluxrad(1:ngridmx)+fluxgrd(1:ngridmx)
    805 
    806          zdum1(1:ngridmx,1:nlayermx)=0.0
    807          zdum2(1:ngridmx,1:nlayermx)=0.0
     839         zflubid(1:ngrid)=fluxrad(1:ngrid)+fluxgrd(1:ngrid)
     840
     841         zdum1(1:ngrid,1:nlayermx)=0.0
     842         zdum2(1:ngrid,1:nlayermx)=0.0
    808843
    809844
    … …  
    821856         else
    822857         
    823            zdh(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
     858           zdh(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)/zpopsk(1:ngrid,1:nlayermx)
    824859 
    825860           call vdifc(ngrid,nlayer,nq,rnat,zpopsk,   &
    … …  
    831866              sensibFlux,q2,zdqdif,zdqsdif,lastcall)
    832867
    833            zdtdif(1:ngridmx,1:nlayermx)=zdhdif(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
    834            zdqevap(1:ngridmx,1:nlayermx)=0.
     868           zdtdif(1:ngrid,1:nlayermx)=zdhdif(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
     869           zdqevap(1:ngrid,1:nlayermx)=0.
    835870
    836871         end if !turbdiff
    837872
    838          pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvdif(1:ngridmx,1:nlayermx)
    839          pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zdudif(1:ngridmx,1:nlayermx)
    840          pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtdif(1:ngridmx,1:nlayermx)
    841          zdtsurf(1:ngridmx)=zdtsurf(1:ngridmx)+zdtsdif(1:ngridmx)
     873         pdv(1:ngrid,1:nlayermx)=pdv(1:ngrid,1:nlayermx)+zdvdif(1:ngrid,1:nlayermx)
     874         pdu(1:ngrid,1:nlayermx)=pdu(1:ngrid,1:nlayermx)+zdudif(1:ngrid,1:nlayermx)
     875         pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtdif(1:ngrid,1:nlayermx)
     876         zdtsurf(1:ngrid)=zdtsurf(1:ngrid)+zdtsdif(1:ngrid)
    842877         if (tracer) then
    843            pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqdif(1:ngridmx,1:nlayermx,1:nqmx)
    844            dqsurf(1:ngridmx,1:nqmx)=dqsurf(1:ngridmx,1:nqmx) + zdqsdif(1:ngridmx,1:nqmx)
     878           pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqdif(1:ngrid,1:nlayermx,1:nq)
     879           dqsurf(1:ngrid,1:nq)=dqsurf(1:ngrid,1:nq) + zdqsdif(1:ngrid,1:nq)
    845880         end if ! of if (tracer)
    846881
    … …  
    899934         if(.not.newtonian)then
    900935
    901             zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + (fluxrad(1:ngridmx) + fluxgrd(1:ngridmx))/capcal(1:ngridmx)
    902 
    903          endif
    904 
     936            zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + (fluxrad(1:ngrid) + fluxgrd(1:ngrid))/capcal(1:ngrid)
     937
     938         endif
     939
     940      print*,'end TURBULENCE'
    905941      endif ! of if (calldifv)
    906942
    … …  
    912948      if(calladj) then
    913949
    914          zdh(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)/zpopsk(1:ngridmx,1:nlayermx)
    915          zduadj(1:ngridmx,1:nlayermx)=0.0
    916          zdvadj(1:ngridmx,1:nlayermx)=0.0
    917          zdhadj(1:ngridmx,1:nlayermx)=0.0
     950         zdh(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)/zpopsk(1:ngrid,1:nlayermx)
     951         zduadj(1:ngrid,1:nlayermx)=0.0
     952         zdvadj(1:ngrid,1:nlayermx)=0.0
     953         zdhadj(1:ngrid,1:nlayermx)=0.0
    918954
    919955
    … …  
    925961              zdqadj)
    926962
    927          pdu(1:ngridmx,1:nlayermx) = pdu(1:ngridmx,1:nlayermx) + zduadj(1:ngridmx,1:nlayermx)
    928          pdv(1:ngridmx,1:nlayermx) = pdv(1:ngridmx,1:nlayermx) + zdvadj(1:ngridmx,1:nlayermx)
    929          pdt(1:ngridmx,1:nlayermx)    = pdt(1:ngridmx,1:nlayermx) + zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx)
    930          zdtadj(1:ngridmx,1:nlayermx) = zdhadj(1:ngridmx,1:nlayermx)*zpopsk(1:ngridmx,1:nlayermx) ! for diagnostic only
     963         pdu(1:ngrid,1:nlayermx) = pdu(1:ngrid,1:nlayermx) + zduadj(1:ngrid,1:nlayermx)
     964         pdv(1:ngrid,1:nlayermx) = pdv(1:ngrid,1:nlayermx) + zdvadj(1:ngrid,1:nlayermx)
     965         pdt(1:ngrid,1:nlayermx)    = pdt(1:ngrid,1:nlayermx) + zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx)
     966         zdtadj(1:ngrid,1:nlayermx) = zdhadj(1:ngrid,1:nlayermx)*zpopsk(1:ngrid,1:nlayermx) ! for diagnostic only
    931967 
    932968         if(tracer) then
    933             pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqadj(1:ngridmx,1:nlayermx,1:nqmx)
     969            pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqadj(1:ngrid,1:nlayermx,1:nq)
    934970         end if
    935971
    … …  
    959995         endif
    960996         !-------------------------
    961                
     997         
     998      print*,'end CONVADJ'
    962999      endif ! of if(calladj)
    9631000
    … …  
    9801017
    9811018
    982          pdt(1:ngridmx,1:nlayermx)=pdt(1:ngridmx,1:nlayermx)+zdtc(1:ngridmx,1:nlayermx)
    983          pdv(1:ngridmx,1:nlayermx)=pdv(1:ngridmx,1:nlayermx)+zdvc(1:ngridmx,1:nlayermx)
    984          pdu(1:ngridmx,1:nlayermx)=pdu(1:ngridmx,1:nlayermx)+zduc(1:ngridmx,1:nlayermx)
    985          zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurfc(1:ngridmx)
    986 
    987          pdq(1:ngridmx,1:nlayermx,1:nqmx)=pdq(1:ngridmx,1:nlayermx,1:nqmx)+ zdqc(1:ngridmx,1:nlayermx,1:nqmx)
     1019         pdt(1:ngrid,1:nlayermx)=pdt(1:ngrid,1:nlayermx)+zdtc(1:ngrid,1:nlayermx)
     1020         pdv(1:ngrid,1:nlayermx)=pdv(1:ngrid,1:nlayermx)+zdvc(1:ngrid,1:nlayermx)
     1021         pdu(1:ngrid,1:nlayermx)=pdu(1:ngrid,1:nlayermx)+zduc(1:ngrid,1:nlayermx)
     1022         zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurfc(1:ngrid)
     1023
     1024         pdq(1:ngrid,1:nlayermx,1:nq)=pdq(1:ngrid,1:nlayermx,1:nq)+ zdqc(1:ngrid,1:nlayermx,1:nq)
    9881025         ! Note: we do not add surface co2ice tendency
    9891026         ! because qsurf(:,igcm_co2_ice) is updated in condens_co2cloud
    … …  
    10211058!             ----------------
    10221059
    1023                dqmoist(1:ngridmx,1:nlayermx,1:nqmx)=0.
    1024                dtmoist(1:ngridmx,1:nlayermx)=0.
    1025 
    1026                call moistadj(pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
    1027 
    1028                pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)   &
    1029                           +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)
    1030                pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) =pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)     &
    1031                           +dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_ice)
    1032                pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtmoist(1:ngridmx,1:nlayermx)
     1060               dqmoist(1:ngrid,1:nlayermx,1:nq)=0.
     1061               dtmoist(1:ngrid,1:nlayermx)=0.
     1062
     1063               call moistadj(ngrid,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
     1064
     1065               pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap)   &
     1066                          +dqmoist(1:ngrid,1:nlayermx,igcm_h2o_vap)
     1067               pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) =pdq(1:ngrid,1:nlayermx,igcm_h2o_ice)     &
     1068                          +dqmoist(1:ngrid,1:nlayermx,igcm_h2o_ice)
     1069               pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+dtmoist(1:ngrid,1:nlayermx)
    10331070
    10341071               !-------------------------
    … …  
    10361073               if(enertest)then
    10371074                  dEtot=cpp*SUM(massarea(:,:)*dtmoist(:,:))/totarea
    1038                   do ig = 1, ngrid
     1075                  do ig=1,ngrid
    10391076                     madjdE(ig) = cpp*SUM(mass(:,:)*dtmoist(:,:))
    10401077                  enddo
    … …  
    10551092!        --------------------------------
    10561093
    1057                call largescale(ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
    1058 
    1059                pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+dtlscale(1:ngridmx,1:nlayermx)
    1060                pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngridmx,1:nlayermx)
    1061                pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngridmx,1:nlayermx)
     1094               call largescale(ngrid,nq,ptimestep,pplev,pplay,pt,pq,pdt,pdq,dtlscale,dqvaplscale,dqcldlscale,rneb_lsc)
     1095
     1096               pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+dtlscale(1:ngrid,1:nlayermx)
     1097               pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap)+dqvaplscale(1:ngrid,1:nlayermx)
     1098               pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) = pdq(1:ngrid,1:nlayermx,igcm_h2o_ice)+dqcldlscale(1:ngrid,1:nlayermx)
    10621099
    10631100               !-------------------------
    10641101               ! test energy conservation
    10651102               if(enertest)then
    1066                   do ig = 1, ngrid
     1103                  do ig=1,ngrid
    10671104                     lscaledE(ig) = cpp*SUM(mass(:,:)*dtlscale(:,:))
    10681105                  enddo
    … …  
    10831120               ! compute cloud fraction
    10841121               do l = 1, nlayer
    1085                   do ig = 1,ngrid
     1122                  do ig=1,ngrid
    10861123                     cloudfrac(ig,l)=MAX(rneb_lsc(ig,l),rneb_man(ig,l))
    10871124                  enddo
    … …  
    10961133            if(waterrain)then
    10971134
    1098                zdqrain(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
    1099                zdqsrain(1:ngridmx)    = 0.0
    1100                zdqssnow(1:ngridmx)    = 0.0
    1101 
    1102                call rain(ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
     1135               zdqrain(1:ngrid,1:nlayermx,1:nq) = 0.0
     1136               zdqsrain(1:ngrid)    = 0.0
     1137               zdqssnow(1:ngrid)    = 0.0
     1138
     1139               call rain(ngrid,nq,ptimestep,pplev,pplay,pt,pdt,pq,pdq,            &
    11031140                   zdtrain,zdqrain,zdqsrain,zdqssnow,cloudfrac)
    11041141
    1105                pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_vap) &
    1106                      +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)
    1107                pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) = pdq(1:ngridmx,1:nlayermx,igcm_h2o_ice) &
    1108                      +zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_ice)
    1109                pdt(1:ngridmx,1:nlayermx) = pdt(1:ngridmx,1:nlayermx)+zdtrain(1:ngridmx,1:nlayermx)
    1110                dqsurf(1:ngridmx,igcm_h2o_vap) = dqsurf(1:ngridmx,igcm_h2o_vap)+zdqsrain(1:ngridmx) ! a bug was here
    1111                dqsurf(1:ngridmx,igcm_h2o_ice) = dqsurf(1:ngridmx,igcm_h2o_ice)+zdqssnow(1:ngridmx)
    1112                rainout(1:ngridmx)             = zdqsrain(1:ngridmx)+zdqssnow(1:ngridmx) ! diagnostic   
     1142               pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) = pdq(1:ngrid,1:nlayermx,igcm_h2o_vap) &
     1143                     +zdqrain(1:ngrid,1:nlayermx,igcm_h2o_vap)
     1144               pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) = pdq(1:ngrid,1:nlayermx,igcm_h2o_ice) &
     1145                     +zdqrain(1:ngrid,1:nlayermx,igcm_h2o_ice)
     1146               pdt(1:ngrid,1:nlayermx) = pdt(1:ngrid,1:nlayermx)+zdtrain(1:ngrid,1:nlayermx)
     1147               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid) ! a bug was here
     1148               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
     1149               rainout(1:ngrid)             = zdqsrain(1:ngrid)+zdqssnow(1:ngrid) ! diagnostic   
    11131150
    11141151
    … …  
    11471184!        -------------
    11481185        if (sedimentation) then
    1149            zdqsed(1:ngridmx,1:nlayermx,1:nqmx) = 0.0
    1150            zdqssed(1:ngridmx,1:nqmx)  = 0.0
     1186           zdqsed(1:ngrid,1:nlayermx,1:nq) = 0.0
     1187           zdqssed(1:ngrid,1:nq)  = 0.0
    11511188
    11521189
    … …  
    11791216              ! and as in rain.F90, whether it falls as rain or snow depends
    11801217              ! only on the surface temperature
    1181            pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqsed(1:ngridmx,1:nlayermx,1:nqmx)
    1182            dqsurf(1:ngridmx,1:nqmx) = dqsurf(1:ngridmx,1:nqmx) + zdqssed(1:ngridmx,1:nqmx)
     1218           pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqsed(1:ngrid,1:nlayermx,1:nq)
     1219           dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
    11831220
    11841221           !-------------------------
    … …  
    12051242         if(mass_redistrib) then
    12061243
    1207             zdmassmr(1:ngridmx,1:nlayermx) = mass(1:ngridmx,1:nlayermx) * &
    1208                (   zdqevap(1:ngridmx,1:nlayermx)                          &
    1209                  + zdqrain(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
    1210                  + dqmoist(1:ngridmx,1:nlayermx,igcm_h2o_vap)             &
    1211                  + dqvaplscale(1:ngridmx,1:nlayermx) )
     1244            zdmassmr(1:ngrid,1:nlayermx) = mass(1:ngrid,1:nlayermx) * &
     1245               (   zdqevap(1:ngrid,1:nlayermx)                          &
     1246                 + zdqrain(1:ngrid,1:nlayermx,igcm_h2o_vap)             &
     1247                 + dqmoist(1:ngrid,1:nlayermx,igcm_h2o_vap)             &
     1248                 + dqvaplscale(1:ngrid,1:nlayermx) )
    12121249                 
    12131250           
    1214             call writediagfi(ngridmx,"mass_evap","mass gain"," ",3,zdmassmr)
    1215             call writediagfi(ngridmx,"mass","mass"," ",3,mass)
     1251            call writediagfi(ngrid,"mass_evap","mass gain"," ",3,zdmassmr)
     1252            call writediagfi(ngrid,"mass","mass"," ",3,mass)
    12161253
    12171254            call mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
    … …  
    12211258         
    12221259
    1223             pdq(1:ngridmx,1:nlayermx,1:nqmx) = pdq(1:ngridmx,1:nlayermx,1:nqmx) + zdqmr(1:ngridmx,1:nlayermx,1:nqmx)
    1224             dqsurf(1:ngridmx,1:nqmx)         = dqsurf(1:ngridmx,1:nqmx) + zdqsurfmr(1:ngridmx,1:nqmx)
    1225             pdt(1:ngridmx,1:nlayermx)        = pdt(1:ngridmx,1:nlayermx) + zdtmr(1:ngridmx,1:nlayermx)
    1226             pdu(1:ngridmx,1:nlayermx)        = pdu(1:ngridmx,1:nlayermx) + zdumr(1:ngridmx,1:nlayermx)
    1227             pdv(1:ngridmx,1:nlayermx)        = pdv(1:ngridmx,1:nlayermx) + zdvmr(1:ngridmx,1:nlayermx)
    1228             pdpsrf(1:ngridmx)                = pdpsrf(1:ngridmx) + zdpsrfmr(1:ngridmx)
    1229             zdtsurf(1:ngridmx)               = zdtsurf(1:ngridmx) + zdtsurfmr(1:ngridmx)
     1260            pdq(1:ngrid,1:nlayermx,1:nq) = pdq(1:ngrid,1:nlayermx,1:nq) + zdqmr(1:ngrid,1:nlayermx,1:nq)
     1261            dqsurf(1:ngrid,1:nq)         = dqsurf(1:ngrid,1:nq) + zdqsurfmr(1:ngrid,1:nq)
     1262            pdt(1:ngrid,1:nlayermx)        = pdt(1:ngrid,1:nlayermx) + zdtmr(1:ngrid,1:nlayermx)
     1263            pdu(1:ngrid,1:nlayermx)        = pdu(1:ngrid,1:nlayermx) + zdumr(1:ngrid,1:nlayermx)
     1264            pdv(1:ngrid,1:nlayermx)        = pdv(1:ngrid,1:nlayermx) + zdvmr(1:ngrid,1:nlayermx)
     1265            pdpsrf(1:ngrid)                = pdpsrf(1:ngrid) + zdpsrfmr(1:ngrid)
     1266            zdtsurf(1:ngrid)               = zdtsurf(1:ngrid) + zdtsurfmr(1:ngrid)
    12301267           
    12311268!           print*,'after mass redistrib, q=',pq(211,1:nlayermx,igcm_h2o_vap)+ptimestep*pdq(211,1:nlayermx,igcm_h2o_vap)
    … …  
    12401277         if(hydrology)then
    12411278
    1242             call hydrol(ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
     1279            call hydrol(ngrid,nq,ptimestep,rnat,tsurf,qsurf,dqsurf,dqs_hyd,    &
    12431280               capcal,albedo0,albedo,mu0,zdtsurf,zdtsurf_hyd,hice)
    12441281            ! note: for now, also changes albedo in the subroutine
    12451282
    1246             zdtsurf(1:ngridmx) = zdtsurf(1:ngridmx) + zdtsurf_hyd(1:ngridmx)
    1247             qsurf(1:ngridmx,1:nqmx) = qsurf(1:ngridmx,1:nqmx) + ptimestep*dqs_hyd(1:ngridmx,1:nqmx)
     1283            zdtsurf(1:ngrid) = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
     1284            qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqs_hyd(1:ngrid,1:nq)
    12481285            ! when hydrology is used, other dqsurf tendencies are all added to dqs_hyd inside
    12491286
    … …  
    12691306         ELSE     ! of if (hydrology)
    12701307
    1271             qsurf(1:ngridmx,1:nqmx)=qsurf(1:ngridmx,1:nqmx)+ptimestep*dqsurf(1:ngridmx,1:nqmx)
     1308            qsurf(1:ngrid,1:nq)=qsurf(1:ngrid,1:nq)+ptimestep*dqsurf(1:ngrid,1:nq)
    12721309
    12731310         END IF   ! of if (hydrology)
    … …  
    12801317
    12811318         if(ice_update)then
    1282             ice_min(1:ngridmx)=min(ice_min(1:ngridmx),qsurf(1:ngridmx,igcm_h2o_ice))
     1319            ice_min(1:ngrid)=min(ice_min(1:ngrid),qsurf(1:ngrid,igcm_h2o_ice))
    12831320         endif
    12841321
    … …  
    12911328
    12921329!     Increment surface temperature
    1293       tsurf(1:ngridmx)=tsurf(1:ngridmx)+ptimestep*zdtsurf(1:ngridmx)
     1330      tsurf(1:ngrid)=tsurf(1:ngrid)+ptimestep*zdtsurf(1:ngrid)
    12941331
    12951332!     Compute soil temperatures and subsurface heat flux
    12961333      if (callsoil) then
    1297          call soil(ngrid,nsoilmx,.false.,inertiedat,        &
     1334         call soil(ngrid,nsoilmx,.false.,lastcall,inertiedat,   &
    12981335                ptimestep,tsurf,tsoil,capcal,fluxgrd)
    12991336      endif
    … …  
    13171354 
    13181355      ! temperature, zonal and meridional wind
    1319       zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx) + pdt(1:ngridmx,1:nlayermx)*ptimestep
    1320       zu(1:ngridmx,1:nlayermx) = pu(1:ngridmx,1:nlayermx) + pdu(1:ngridmx,1:nlayermx)*ptimestep
    1321       zv(1:ngridmx,1:nlayermx) = pv(1:ngridmx,1:nlayermx) + pdv(1:ngridmx,1:nlayermx)*ptimestep
     1356      zt(1:ngrid,1:nlayermx) = pt(1:ngrid,1:nlayermx) + pdt(1:ngrid,1:nlayermx)*ptimestep
     1357      zu(1:ngrid,1:nlayermx) = pu(1:ngrid,1:nlayermx) + pdu(1:ngrid,1:nlayermx)*ptimestep
     1358      zv(1:ngrid,1:nlayermx) = pv(1:ngrid,1:nlayermx) + pdv(1:ngrid,1:nlayermx)*ptimestep
    13221359
    13231360      ! diagnostic
    1324       zdtdyn(1:ngridmx,1:nlayermx)     = pt(1:ngridmx,1:nlayermx)-ztprevious(1:ngridmx,1:nlayermx)
    1325       ztprevious(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
     1361      zdtdyn(1:ngrid,1:nlayermx)     = pt(1:ngrid,1:nlayermx)-ztprevious(1:ngrid,1:nlayermx)
     1362      ztprevious(1:ngrid,1:nlayermx) = zt(1:ngrid,1:nlayermx)
    13261363
    13271364      if(firstcall)then
    1328          zdtdyn(1:ngridmx,1:nlayermx)=0.0
     1365         zdtdyn(1:ngrid,1:nlayermx)=0.0
    13291366      endif
    13301367
    … …  
    13351372
    13361373      ! tracers
    1337       zq(1:ngridmx,1:nlayermx,1:nqmx) = pq(1:ngridmx,1:nlayermx,1:nqmx) + pdq(1:ngridmx,1:nlayermx,1:nqmx)*ptimestep
     1374      zq(1:ngrid,1:nlayermx,1:nq) = pq(1:ngrid,1:nlayermx,1:nq) + pdq(1:ngrid,1:nlayermx,1:nq)*ptimestep
    13381375
    13391376      ! surface pressure
    1340       ps(1:ngridmx) = pplev(1:ngridmx,1) + pdpsrf(1:ngridmx)*ptimestep
     1377      ps(1:ngrid) = pplev(1:ngrid,1) + pdpsrf(1:ngrid)*ptimestep
    13411378
    13421379      ! pressure
    13431380      do l=1,nlayer
    1344           zplev(1:ngridmx,l) = pplev(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(:)
    1345           zplay(1:ngridmx,l) = pplay(1:ngridmx,l)/pplev(1:ngridmx,1)*ps(1:ngridmx)
     1381          zplev(1:ngrid,l) = pplev(1:ngrid,l)/pplev(1:ngrid,1)*ps(:)
     1382          zplay(1:ngrid,l) = pplay(1:ngrid,l)/pplev(1:ngrid,1)*ps(1:ngrid)
    13461383      enddo
    13471384
    … …  
    13731410         GND = SUM(area(:)*fluxgrd(:))/totarea
    13741411         DYN = SUM(area(:)*fluxdyn(:))/totarea
    1375          do ig=1,ngrid             
     1412         do ig=1,ngrid
    13761413            if(fluxtop_dn(ig).lt.0.0)then
    13771414               print*,'fluxtop_dn has gone crazy'
    … …  
    13851422         end do
    13861423                     
    1387          if(ngridmx.eq.1)then
     1424         if(ngrid.eq.1)then
    13881425            DYN=0.0
    13891426         endif
    … …  
    13991436
    14001437         if(meanOLR)then
    1401             if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
     1438            if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
    14021439               ! to record global radiative balance
    14031440               open(92,file="rad_bal.out",form='formatted',position='append')
    … …  
    14321469!     ---------------------------------------------------------
    14331470      if(tracer)then
    1434          qcol(1:ngridmx,1:nqmx)=0.0
     1471         qcol(1:ngrid,1:nq)=0.0
    14351472         do iq=1,nq
    1436            do ig=1,ngridmx
     1473           do ig=1,ngrid
    14371474              qcol(ig,iq) = SUM( zq(ig,1:nlayermx,iq) * mass(ig,1:nlayermx))
    14381475           enddo
    … …  
    14401477
    14411478         ! Generalised for arbitrary aerosols now. (LK)
    1442          reffcol(1:ngridmx,1:naerkind)=0.0
     1479         reffcol(1:ngrid,1:naerkind)=0.0
    14431480         if(co2cond.and.(iaero_co2.ne.0))then
    1444             call co2_reffrad(zq,reffrad)
    1445             do ig=1,ngridmx
     1481            call co2_reffrad(ngrid,nq,zq,reffrad)
     1482            do ig=1,ngrid
    14461483               reffcol(ig,iaero_co2) = SUM(zq(ig,1:nlayermx,igcm_co2_ice)*reffrad(ig,1:nlayermx,iaero_co2)*mass(ig,1:nlayermx))
    14471484            enddo
    14481485         endif
    14491486         if(water.and.(iaero_h2o.ne.0))then
    1450             call h2o_reffrad(zq,zt,reffrad,nueffrad)
    1451             do ig=1,ngridmx
     1487
     1488            !! AS: This modifies reffrad and nueffrad
     1489            !!     ... and those are saved in physiq
     1490            !!     To be conservative with previous versions,
     1491            !!     ... I put nueffrad_dummy so that nueffrad
     1492            !!     ... is not modified, but shouldn't it be modified actually?
     1493            call h2o_reffrad(ngrid,nq,zq,zt,reffrad,nueffrad_dummy)
     1494            do ig=1,ngrid
    14521495               reffcol(ig,iaero_h2o) = SUM(zq(ig,1:nlayermx,igcm_h2o_ice)*reffrad(ig,1:nlayermx,iaero_h2o)*mass(ig,1:nlayermx))
    14531496            enddo
    … …  
    14741517
    14751518         if(meanOLR)then
    1476             if((ngridmx.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
     1519            if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
    14771520               ! to record global water balance
    14781521               open(98,file="h2o_bal.out",form='formatted',position='append')
    … …  
    14901533      if(water)then
    14911534         do l = 1, nlayer
    1492             do ig = 1, ngrid
     1535            do ig=1,ngrid
    14931536!               call watersat(pt(ig,l),pplay(ig,l),qsat(ig,l))
    14941537               call Psat_water(zt(ig,l),pplay(ig,l),psat_tmp,qsat(ig,l))
    … …  
    15241567            if(ice_update)then
    15251568
    1526                do ig = 1, ngrid
     1569               do ig=1,ngrid
    15271570
    15281571                  delta_ice = (qsurf(ig,igcm_h2o_ice)-ice_initial(ig))
    … …  
    15501593
    15511594            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
    1552             call physdem1("restartfi.nc",long,lati,nsoilmx,nq,            &
     1595            call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
    15531596                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
    15541597                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
    1555                     cloudfrac,totcloudfrac,hice)
     1598                    cloudfrac,totcloudfrac,hice,noms)
    15561599         endif
    15571600
    … …  
    15941637             do iq=1,nq
    15951638                call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    1596                 call wstats(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     1639                call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    15971640                          'kg m^-2',2,qsurf(1,iq) )
    15981641
    1599                 call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
     1642                call wstats(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    16001643                          'kg m^-2',2,qcol(1,iq) )
    1601 !                call wstats(ngridmx,trim(noms(iq))//'_reff',                          &
     1644!                call wstats(ngrid,trim(noms(iq))//'_reff',                          &
    16021645!                          trim(noms(iq))//'_reff',                                   &
    16031646!                          'm',3,reffrad(1,1,iq))
    16041647              end do
    16051648            if (water) then
    1606                vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
     1649               vmr=zq(1:ngrid,1:nlayermx,igcm_h2o_vap)*mugaz/mmol(igcm_h2o_vap)
    16071650               call wstats(ngrid,"vmr_h2ovapor",                           &
    16081651                          "H2O vapour volume mixing ratio","mol/mol",       &
    … …  
    16471690        if(enertest) then
    16481691          if (calldifv) then
    1649              call writediagfi(ngridmx,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
    1650 !            call writediagfi(ngridmx,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
    1651 !            call writediagfi(ngridmx,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
    1652 !            call writediagfi(ngridmx,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
    1653              call writediagfi(ngridmx,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
     1692             call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
     1693!            call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
     1694!            call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
     1695!            call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
     1696             call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
    16541697          endif
    16551698          if (corrk) then
    1656 !            call writediagfi(ngridmx,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
    1657 !            call writediagfi(ngridmx,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
     1699!            call writediagfi(ngrid,"dEzradsw","radiative heating","w.m^-2",3,dEzradsw)
     1700!            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
    16581701          endif
    16591702          if(watercond) then
    16601703             call writediagfi(ngrid,"lscaledE","heat from largescale","W m-2",2,lscaledE)
    16611704             call writediagfi(ngrid,"madjdE","heat from moistadj","W m-2",2,madjdE)
    1662              call writediagfi(ngridmx,"RH","relative humidity"," ",3,RH)
    1663 !            call writediagfi(ngridmx,"qsatatm","atm qsat"," ",3,qsat)
    1664 !            call writediagfi(ngridmx,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
     1705             call writediagfi(ngrid,"RH","relative humidity"," ",3,RH)
     1706!            call writediagfi(ngrid,"qsatatm","atm qsat"," ",3,qsat)
     1707!            call writediagfi(ngrid,"h2o_max_col","maximum H2O column amount","kg.m^-2",2,H2Omaxcol)
    16651708          endif
    16661709        endif
    … …  
    16781721!     Output aerosols
    16791722        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
    1680           call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
     1723          call writediagfi(ngrid,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,iaero_co2))
    16811724        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
    1682           call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
     1725          call writediagfi(ngrid,'H2Oice_reff','H2Oice_reff','m',3,reffrad(:,:,iaero_h2o))
    16831726        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0) &
    1684           call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
     1727          call writediagfi(ngrid,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,iaero_co2))
    16851728        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0) &
    1686           call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
     1729          call writediagfi(ngrid,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,iaero_h2o))
    16871730
    16881731!     Output tracers
    … …  
    16901733        do iq=1,nq
    16911734          call writediagfi(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
    1692 !          call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     1735!          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    16931736!                          'kg m^-2',2,qsurf(1,iq) )
    1694           call writediagfi(ngridmx,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
     1737          call writediagfi(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  &
    16951738                          'kg m^-2',2,qsurf_hist(1,iq) )
    1696           call writediagfi(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
     1739          call writediagfi(ngrid,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    &
    16971740                          'kg m^-2',2,qcol(1,iq) )
    16981741
    … …  
    17051748
    17061749          if(waterrain)then
    1707              call writediagfi(ngridmx,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
    1708              call writediagfi(ngridmx,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
     1750             call writediagfi(ngrid,"rain","rainfall","kg m-2 s-1",2,zdqsrain)
     1751             call writediagfi(ngrid,"snow","snowfall","kg m-2 s-1",2,zdqssnow)
    17091752          endif
    17101753
    17111754          if(hydrology)then
    1712              call writediagfi(ngridmx,"hice","oceanic ice height","m",2,hice)
     1755             call writediagfi(ngrid,"hice","oceanic ice height","m",2,hice)
    17131756          endif
    17141757
    17151758          if(ice_update)then
    1716              call writediagfi(ngridmx,"ice_min","min annual ice","m",2,ice_min)
    1717              call writediagfi(ngridmx,"ice_ini","initial annual ice","m",2,ice_initial)
     1759             call writediagfi(ngrid,"ice_min","min annual ice","m",2,ice_min)
     1760             call writediagfi(ngrid,"ice_ini","initial annual ice","m",2,ice_initial)
    17181761          endif
    17191762
    … …  
    17251768          end do
    17261769
    1727           call writediagfi(ngridmx,"tau_col","Total aerosol optical depth","[]",2,tau_col)
     1770          call writediagfi(ngrid,"tau_col","Total aerosol optical depth","[]",2,tau_col)
    17281771
    17291772        enddo
    … …  
    17391782      icount=icount+1
    17401783
    1741       ! deallocate gas variables
    17421784      if (lastcall) then
     1785
     1786        ! deallocate gas variables
    17431787        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
    17441788        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
     1789
     1790        ! deallocate saved arrays
     1791        ! this is probably not that necessary
     1792        ! ... but probably a good idea to clean the place before we leave
     1793        IF ( ALLOCATED(tsurf)) DEALLOCATE(tsurf)
     1794        IF ( ALLOCATED(tsoil)) DEALLOCATE(tsoil)
     1795        IF ( ALLOCATED(albedo)) DEALLOCATE(albedo)
     1796        IF ( ALLOCATED(albedo0)) DEALLOCATE(albedo0)
     1797        IF ( ALLOCATED(rnat)) DEALLOCATE(rnat)
     1798        IF ( ALLOCATED(emis)) DEALLOCATE(emis)
     1799        IF ( ALLOCATED(dtrad)) DEALLOCATE(dtrad)
     1800        IF ( ALLOCATED(fluxrad_sky)) DEALLOCATE(fluxrad_sky)
     1801        IF ( ALLOCATED(fluxrad)) DEALLOCATE(fluxrad)
     1802        IF ( ALLOCATED(capcal)) DEALLOCATE(capcal)
     1803        IF ( ALLOCATED(fluxgrd)) DEALLOCATE(fluxgrd)
     1804        IF ( ALLOCATED(qsurf)) DEALLOCATE(qsurf)
     1805        IF ( ALLOCATED(q2)) DEALLOCATE(q2)
     1806        IF ( ALLOCATED(ztprevious)) DEALLOCATE(ztprevious)
     1807        IF ( ALLOCATED(cloudfrac)) DEALLOCATE(cloudfrac)
     1808        IF ( ALLOCATED(totcloudfrac)) DEALLOCATE(totcloudfrac)
     1809        IF ( ALLOCATED(qsurf_hist)) DEALLOCATE(qsurf_hist)
     1810        IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
     1811        IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)
     1812        IF ( ALLOCATED(ice_initial)) DEALLOCATE(ice_initial)
     1813        IF ( ALLOCATED(ice_min)) DEALLOCATE(ice_min)
     1814
     1815        IF ( ALLOCATED(fluxsurf_lw)) DEALLOCATE(fluxsurf_lw)
     1816        IF ( ALLOCATED(fluxsurf_sw)) DEALLOCATE(fluxsurf_sw)
     1817        IF ( ALLOCATED(fluxtop_lw)) DEALLOCATE(fluxtop_lw)
     1818        IF ( ALLOCATED(fluxabs_sw)) DEALLOCATE(fluxabs_sw)
     1819        IF ( ALLOCATED(fluxtop_dn)) DEALLOCATE(fluxtop_dn)
     1820        IF ( ALLOCATED(fluxdyn)) DEALLOCATE(fluxdyn)
     1821        IF ( ALLOCATED(OLR_nu)) DEALLOCATE(OLR_nu)
     1822        IF ( ALLOCATED(OSR_nu)) DEALLOCATE(OSR_nu)
     1823        IF ( ALLOCATED(sensibFlux)) DEALLOCATE(sensibFlux)
     1824        IF ( ALLOCATED(zdtlw)) DEALLOCATE(zdtlw)
     1825        IF ( ALLOCATED(zdtsw)) DEALLOCATE(zdtsw)
     1826        IF ( ALLOCATED(tau_col)) DEALLOCATE(tau_col)
     1827
     1828        !! this is defined in comsaison_h
     1829        IF ( ALLOCATED(mu0)) DEALLOCATE(mu0)
     1830        IF ( ALLOCATED(fract)) DEALLOCATE(fract)
     1831
     1832        !! this is defined in comsoil_h
     1833        IF ( ALLOCATED(layer)) DEALLOCATE(layer)
     1834        IF ( ALLOCATED(mlayer)) DEALLOCATE(mlayer)
     1835        IF ( ALLOCATED(inertiedat)) DEALLOCATE(inertiedat)
     1836
     1837        !! this is defined in surfdat_h
     1838        IF ( ALLOCATED(albedodat)) DEALLOCATE(albedodat)
     1839        IF ( ALLOCATED(phisfi)) DEALLOCATE(phisfi)
     1840        IF ( ALLOCATED(zmea)) DEALLOCATE(zmea)
     1841        IF ( ALLOCATED(zstd)) DEALLOCATE(zstd)
     1842        IF ( ALLOCATED(zsig)) DEALLOCATE(zsig)
     1843        IF ( ALLOCATED(zgam)) DEALLOCATE(zgam)
     1844        IF ( ALLOCATED(zthe)) DEALLOCATE(zthe)
     1845        IF ( ALLOCATED(dryness)) DEALLOCATE(dryness)
     1846        IF ( ALLOCATED(watercaptag)) DEALLOCATE(watercaptag)
     1847 
     1848        !! this is defined in tracer_h
     1849        IF ( ALLOCATED(noms)) DEALLOCATE(noms)
     1850        IF ( ALLOCATED(mmol)) DEALLOCATE(mmol)
     1851        IF ( ALLOCATED(radius)) DEALLOCATE(radius)
     1852        IF ( ALLOCATED(rho_q)) DEALLOCATE(rho_q)
     1853        IF ( ALLOCATED(qext)) DEALLOCATE(qext)
     1854        IF ( ALLOCATED(alpha_lift)) DEALLOCATE(alpha_lift)
     1855        IF ( ALLOCATED(alpha_devil)) DEALLOCATE(alpha_devil)
     1856        IF ( ALLOCATED(qextrhor)) DEALLOCATE(qextrhor)
     1857        IF ( ALLOCATED(igcm_dustbin)) DEALLOCATE(igcm_dustbin)
     1858
     1859        !! this is defined in comgeomfi_h
     1860        IF ( ALLOCATED(lati)) DEALLOCATE(lati)
     1861        IF ( ALLOCATED(long)) DEALLOCATE(long)
     1862        IF ( ALLOCATED(area)) DEALLOCATE(area)
     1863        IF ( ALLOCATED(sinlat)) DEALLOCATE(sinlat)
     1864        IF ( ALLOCATED(coslat)) DEALLOCATE(coslat)
     1865        IF ( ALLOCATED(sinlon)) DEALLOCATE(sinlon)
     1866        IF ( ALLOCATED(coslon)) DEALLOCATE(coslon)
     1867
    17451868      endif
    1746 
    17471869
    17481870      return

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

    r728 r787  
    2020
    2121!==================================================================
    22    subroutine su_aer_radii(reffrad,nueffrad)
     22   subroutine su_aer_radii(ngrid,reffrad,nueffrad)
    2323!==================================================================
    2424!     Purpose
    … …  
    3535      use radinc_h, only: naerkind
    3636      use aerosol_mod
    37       Implicit none
    38 
    39 #include "callkeys.h"
    40 #include "dimensions.h"
    41 #include "dimphys.h"
    42 #include "tracer.h"
    43 
    44       real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
    45       real, intent(inout) :: nueffrad(ngridmx,nlayermx,naerkind)     !variance     
     37      USE tracer_h
     38      Implicit none
     39
     40#include "callkeys.h"
     41#include "dimensions.h"
     42#include "dimphys.h"
     43
     44      integer :: ngrid
     45
     46      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
     47      real, intent(inout) :: nueffrad(ngrid,nlayermx,naerkind)     !variance     
    4648
    4749      logical, save :: firstcall=.true.
    … …  
    5658
    5759            if(iaer.eq.iaero_co2)then ! CO2 ice
    58                reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-4
    59                nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1
     60               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-4
     61               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
    6062            endif
    6163
    6264            if(iaer.eq.iaero_h2o)then ! H2O ice
    63                reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-5
    64                nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1
     65               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-5
     66               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
    6567            endif
    6668
    6769            if(iaer.eq.iaero_dust)then ! dust
    68                reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-5
    69                nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1
     70               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-5
     71               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
    7072            endif
    7173 
    7274            if(iaer.eq.iaero_h2so4)then ! H2O ice
    73                reffrad(1:ngridmx,1:nlayermx,iaer) = 1.e-6
    74                nueffrad(1:ngridmx,1:nlayermx,iaer) = 0.1
     75               reffrad(1:ngrid,1:nlayermx,iaer) = 1.e-6
     76               nueffrad(1:ngrid,1:nlayermx,iaer) = 0.1
    7577            endif
    7678
    … …  
    117119
    118120!==================================================================
    119    subroutine h2o_reffrad(pq,pt,reffrad,nueffrad)
     121   subroutine h2o_reffrad(ngrid,nq,pq,pt,reffrad,nueffrad)
    120122!==================================================================
    121123!     Purpose
    … …  
    131133      use watercommon_h, Only: T_h2O_ice_liq,T_h2O_ice_clouds,rhowater,rhowaterice
    132134      use aerosol_mod, only : iaero_h2o
    133       Implicit none
    134 
    135 #include "callkeys.h"
    136 #include "dimensions.h"
    137 #include "dimphys.h"
    138 #include "tracer.h"
     135      USE tracer_h
     136      Implicit none
     137
     138#include "callkeys.h"
     139#include "dimensions.h"
     140#include "dimphys.h"
    139141#include "comcstfi.h"
    140142
    141       real, intent(in) :: pq(ngridmx,nlayermx,nqmx) !tracer mixing ratios (kg/kg)
    142       real, intent(in) :: pt(ngridmx,nlayermx)      !temperature (K)
    143       real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
    144       real, intent(inout) :: nueffrad(ngridmx,nlayermx,naerkind) ! dispersion     
     143      integer :: ngrid,nq
     144
     145      real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
     146      real, intent(in) :: pt(ngrid,nlayermx)      !temperature (K)
     147      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
     148      real, intent(inout) :: nueffrad(ngrid,nlayermx,naerkind) ! dispersion     
    145149
    146150      integer :: ig,l
    … …  
    151155      if (radfixed) then
    152156         do l=1,nlayermx
    153             do ig=1,ngridmx
     157            do ig=1,ngrid
    154158               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
    155159               zfice = MIN(MAX(zfice,0.0),1.0)
    … …  
    160164      else
    161165         do l=1,nlayermx
    162             do ig=1,ngridmx
     166            do ig=1,ngrid
    163167               zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
    164168               zfice = MIN(MAX(zfice,0.0),1.0)
    … …  
    177181
    178182!==================================================================
    179    subroutine h2o_cloudrad(pql,reffliq,reffice)
     183   subroutine h2o_cloudrad(ngrid,pql,reffliq,reffice)
    180184!==================================================================
    181185!     Purpose
    … …  
    189193!==================================================================
    190194      use watercommon_h, Only: rhowater,rhowaterice
    191       Implicit none
    192 
    193 #include "callkeys.h"
    194 #include "dimensions.h"
    195 #include "dimphys.h"
    196 #include "tracer.h"
     195      USE tracer_h
     196      Implicit none
     197
     198#include "callkeys.h"
     199#include "dimensions.h"
     200#include "dimphys.h"
    197201#include "comcstfi.h"
    198202
    199       real, intent(in) :: pql(ngridmx,nlayermx) !condensed water mixing ratios (kg/kg)
    200       real, intent(out) :: reffliq(ngridmx,nlayermx),reffice(ngridmx,nlayermx)     !liquid and ice water particle radii (K)
     203      integer :: ngrid
     204
     205      real, intent(in) :: pql(ngrid,nlayermx) !condensed water mixing ratios (kg/kg)
     206      real, intent(out) :: reffliq(ngrid,nlayermx),reffice(ngrid,nlayermx)     !liquid and ice water particle radii (K)
    201207
    202208      real,external :: CBRT           
    … …  
    204210
    205211      if (radfixed) then
    206          reffliq(1:ngridmx,1:nlayermx)= rad_h2o
    207          reffice(1:ngridmx,1:nlayermx)= rad_h2o_ice
     212         reffliq(1:ngrid,1:nlayermx)= rad_h2o
     213         reffice(1:ngrid,1:nlayermx)= rad_h2o_ice
    208214      else
    209          reffliq(1:ngridmx,1:nlayermx)  = CBRT( 3*pql(1:ngridmx,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
    210          reffliq(1:ngridmx,1:nlayermx)  = min(max(reffliq(1:ngridmx,1:nlayermx),1.e-6),100.e-6)
    211          reffice(1:ngridmx,1:nlayermx)  = CBRT( 3*pql(1:ngridmx,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
    212          reffice(1:ngridmx,1:nlayermx)  = min(max(reffice(1:ngridmx,1:nlayermx),1.e-6),100.e-6)
     215         reffliq(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o*pi*rhowater) )
     216         reffliq(1:ngrid,1:nlayermx)  = min(max(reffliq(1:ngrid,1:nlayermx),1.e-6),100.e-6)
     217         reffice(1:ngrid,1:nlayermx)  = CBRT( 3*pql(1:ngrid,1:nlayermx)/(4*Nmix_h2o_ice*pi*rhowaterice) )
     218         reffice(1:ngrid,1:nlayermx)  = min(max(reffice(1:ngrid,1:nlayermx),1.e-6),100.e-6)
    213219      end if
    214220
    … …  
    219225
    220226!==================================================================
    221    subroutine co2_reffrad(pq,reffrad)
     227   subroutine co2_reffrad(ngrid,nq,pq,reffrad)
    222228!==================================================================
    223229!     Purpose
    … …  
    232238      use radinc_h, only: naerkind
    233239      use aerosol_mod, only : iaero_co2
    234       Implicit none
    235 
    236 #include "callkeys.h"
    237 #include "dimensions.h"
    238 #include "dimphys.h"
    239 #include "tracer.h"
     240      USE tracer_h
     241      Implicit none
     242
     243#include "callkeys.h"
     244#include "dimensions.h"
     245#include "dimphys.h"
    240246#include "comcstfi.h"
    241247
    242       real, intent(in) :: pq(ngridmx,nlayermx,nqmx) !tracer mixing ratios (kg/kg)
    243       real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
     248      integer :: ngrid,nq
     249
     250      real, intent(in) :: pq(ngrid,nlayermx,nq) !tracer mixing ratios (kg/kg)
     251      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
    244252
    245253      integer :: ig,l
    … …  
    250258
    251259      if (radfixed) then
    252          reffrad(1:ngridmx,1:nlayermx,iaero_co2) = 5.e-5 ! CO2 ice
     260         reffrad(1:ngrid,1:nlayermx,iaero_co2) = 5.e-5 ! CO2 ice
    253261      else
    254262         do l=1,nlayermx
    255             do ig=1,ngridmx
     263            do ig=1,ngrid
    256264               zrad = CBRT( 3*pq(ig,l,igcm_co2_ice)/(4*Nmix_co2*pi*rho_co2) )
    257265               reffrad(ig,l,iaero_co2) = min(max(zrad,1.e-6),100.e-6)
    … …  
    266274
    267275!==================================================================
    268    subroutine dust_reffrad(reffrad)
     276   subroutine dust_reffrad(ngrid,reffrad)
    269277!==================================================================
    270278!     Purpose
    … …  
    285293#include "dimphys.h"
    286294
    287       real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
     295      integer :: ngrid
     296
     297      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
    288298           
    289       reffrad(1:ngridmx,1:nlayermx,iaero_dust) = 2.e-6 ! dust
     299      reffrad(1:ngrid,1:nlayermx,iaero_dust) = 2.e-6 ! dust
    290300
    291301   end subroutine dust_reffrad
    … …  
    294304
    295305!==================================================================
    296    subroutine h2so4_reffrad(reffrad)
     306   subroutine h2so4_reffrad(ngrid,reffrad)
    297307!==================================================================
    298308!     Purpose
    … …  
    313323#include "dimphys.h"
    314324
    315       real, intent(inout) :: reffrad(ngridmx,nlayermx,naerkind)      !aerosols radii (K)
     325      integer :: ngrid
     326
     327      real, intent(inout) :: reffrad(ngrid,nlayermx,naerkind)      !aerosols radii (K)
    316328               
    317       reffrad(1:ngridmx,1:nlayermx,iaero_h2so4) = 1.e-6 ! h2so4
     329      reffrad(1:ngrid,1:nlayermx,iaero_h2so4) = 1.e-6 ! h2so4
    318330
    319331   end subroutine h2so4_reffrad

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

    r731 r787  
    1 subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
     1subroutine rain(ngrid,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
    22
    33
    … …  
    66  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
    77  use radii_mod, only: h2o_cloudrad
     8  USE tracer_h
    89  implicit none
    910
    … …  
    2425#include "dimensions.h"
    2526#include "dimphys.h"
    26 #include "tracer.h"
    2727#include "comcstfi.h"
    2828#include "callkeys.h"
    2929
     30      integer ngrid,nq
     31
    3032!     Pre-arguments (for universal model)
    31       real pq(ngridmx,nlayermx,nqmx)       ! tracer (kg/kg)
    32       real qsurf(ngridmx,nqmx)             ! tracer at the surface (kg.m-2)
    33       REAL pdt(ngridmx,nlayermx),pdq(ngridmx,nlayermx,nqmx)
    34 
    35       real dqrain(ngridmx,nlayermx,nqmx)   ! tendency of H2O precipitation (kg/kg.s-1)
    36       real dqsrain(ngridmx)                ! rain flux at the surface (kg.m-2.s-1)
    37       real dqssnow(ngridmx)                ! snow flux at the surface (kg.m-2.s-1)
    38       REAL d_t(ngridmx,nlayermx)           ! temperature increment
     33      real pq(ngrid,nlayermx,nq)       ! tracer (kg/kg)
     34      real qsurf(ngrid,nq)             ! tracer at the surface (kg.m-2)
     35      REAL pdt(ngrid,nlayermx),pdq(ngrid,nlayermx,nq)
     36
     37      real dqrain(ngrid,nlayermx,nq)   ! tendency of H2O precipitation (kg/kg.s-1)
     38      real dqsrain(ngrid)                ! rain flux at the surface (kg.m-2.s-1)
     39      real dqssnow(ngrid)                ! snow flux at the surface (kg.m-2.s-1)
     40      REAL d_t(ngrid,nlayermx)           ! temperature increment
    3941
    4042!     Arguments
    4143      REAL ptimestep                    ! time interval
    42       REAL pplev(ngridmx,nlayermx+1)    ! inter-layer pressure
    43       REAL pplay(ngridmx,nlayermx)      ! mid-layer pressure
    44       REAL t(ngridmx,nlayermx)          ! input temperature (K)
    45       REAL zt(ngridmx,nlayermx)         ! working temperature (K)
    46       REAL ql(ngridmx,nlayermx)         ! liquid water (Kg/Kg)
    47       REAL q(ngridmx,nlayermx)          ! specific humidity (Kg/Kg)
    48       REAL rneb(ngridmx,nlayermx)       ! cloud fraction
    49       REAL d_q(ngridmx,nlayermx)        ! water vapor increment
    50       REAL d_ql(ngridmx,nlayermx)       ! liquid water / ice increment
     44      REAL pplev(ngrid,nlayermx+1)    ! inter-layer pressure
     45      REAL pplay(ngrid,nlayermx)      ! mid-layer pressure
     46      REAL t(ngrid,nlayermx)          ! input temperature (K)
     47      REAL zt(ngrid,nlayermx)         ! working temperature (K)
     48      REAL ql(ngrid,nlayermx)         ! liquid water (Kg/Kg)
     49      REAL q(ngrid,nlayermx)          ! specific humidity (Kg/Kg)
     50      REAL rneb(ngrid,nlayermx)       ! cloud fraction
     51      REAL d_q(ngrid,nlayermx)        ! water vapor increment
     52      REAL d_ql(ngrid,nlayermx)       ! liquid water / ice increment
    5153
    5254!     Subroutine options
    … …  
    7880!     Local variables
    7981      INTEGER i, k, n
    80       REAL zqs(ngridmx,nlayermx),Tsat(ngridmx,nlayermx), zdelta, zcor
    81       REAL zrfl(ngridmx), zrfln(ngridmx), zqev, zqevt
    82 
    83       REAL zoliq(ngridmx)
    84       REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx)
    85       REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx)
    86 
    87       real reffh2oliq(ngridmx,nlayermx),reffh2oice(ngridmx,nlayermx)
     82      REAL zqs(ngrid,nlayermx),Tsat(ngrid,nlayermx), zdelta, zcor
     83      REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
     84
     85      REAL zoliq(ngrid)
     86      REAL zdz(ngrid),zrho(ngrid),ztot(ngrid), zrhol(ngrid)
     87      REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid)
     88
     89      real reffh2oliq(ngrid,nlayermx),reffh2oice(ngrid,nlayermx)
    8890 
    8991      real ttemp, ptemp, psat_tmp
    90       real tnext(ngridmx,nlayermx)
    91 
    92       real l2c(ngridmx,nlayermx)
     92      real tnext(ngrid,nlayermx)
     93
     94      real l2c(ngrid,nlayermx)
    9395      real dWtot
    9496
    … …  
    154156!     GCM -----> subroutine variables
    155157      DO k = 1, nlayermx
    156       DO i = 1, ngridmx
     158      DO i = 1, ngrid
    157159
    158160         zt(i,k)   = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
    … …  
    175177!     Initialise the outputs
    176178      DO k = 1, nlayermx
    177       DO i = 1, ngridmx
     179      DO i = 1, ngrid
    178180         d_t(i,k) = 0.0
    179181         d_q(i,k) = 0.0
    … …  
    181183      ENDDO
    182184      ENDDO
    183       DO i = 1, ngridmx
     185      DO i = 1, ngrid
    184186         zrfl(i)  = 0.0
    185187         zrfln(i) = 0.0
    … …  
    188190      ! calculate saturation mixing ratio
    189191      DO k = 1, nlayermx
    190          DO i = 1, ngridmx
     192         DO i = 1, ngrid
    191193            ttemp = zt(i,k)
    192194            ptemp = pplay(i,k)
    … …  
    199201      ! get column / layer conversion factor
    200202      DO k = 1, nlayermx
    201          DO i = 1, ngridmx
     203         DO i = 1, ngrid
    202204            l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
    203205         ENDDO
    … …  
    211213
    212214         IF (evap_prec) THEN ! note no rneb dependence!
    213             DO i = 1, ngridmx
     215            DO i = 1, ngrid
    214216               IF (zrfl(i) .GT.0.) THEN
    215217
    … …  
    241243         ENDIF
    242244
    243          DO i = 1, ngridmx
     245         DO i = 1, ngrid
    244246            zoliq(i) = 0.0
    245247         ENDDO
    … …  
    248250         if(precip_scheme.eq.1)then
    249251
    250             DO i = 1, ngridmx
     252            DO i = 1, ngrid
    251253               ttemp = zt(i,k)
    252254               IF (ttemp .ge. T_h2O_ice_liq) THEN
    … …  
    271273         elseif (precip_scheme.ge.2) then
    272274         
    273            DO i = 1, ngridmx
     275           DO i = 1, ngrid
    274276               IF (rneb(i,k).GT.0.0) THEN
    275277                  zoliq(i) = ql(i,k)
    … …  
    284286
    285287 !recalculate liquid water particle radii
    286            call h2o_cloudrad(ql,reffh2oliq,reffh2oice)
     288           call h2o_cloudrad(ngrid,ql,reffh2oliq,reffh2oice)
    287289
    288290           SELECT CASE(precip_scheme)
    … …  
    291293
    292294            DO n = 1, ninter
    293                DO i = 1, ngridmx
     295               DO i = 1, ngrid
    294296                  IF (rneb(i,k).GT.0.0) THEN
    295297                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
    … …  
    314316           
    315317            DO n = 1, ninter
    316                DO i = 1, ngridmx
     318               DO i = 1, ngrid
    317319                  IF (rneb(i,k).GT.0.0) THEN
    318320                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
    … …  
    336338
    337339            DO n = 1, ninter
    338                DO i = 1, ngridmx
     340               DO i = 1, ngrid
    339341                  IF (rneb(i,k).GT.0.0) THEN
    340342                     ! this is the ONLY place where zneb and c1 are used
    … …  
    358360
    359361            ! Change in cloud density and surface H2O values
    360             DO i = 1, ngridmx
     362            DO i = 1, ngrid
    361363               IF (rneb(i,k).GT.0.0) THEN
    362364                  d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
    … …  
    371373
    372374!     Rain or snow on the ground
    373       DO i = 1, ngridmx
     375      DO i = 1, ngrid
    374376         if(zrfl(i).lt.0.0)then
    375377            print*,'Droplets of negative rain are falling...'
    … …  
    387389!     now subroutine -----> GCM variables
    388390      DO k = 1, nlayermx
    389          DO i = 1, ngridmx
     391         DO i = 1, ngrid
    390392           
    391393            if(evap_prec)then

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

    r773 r787  
    33! to use  'getin'
    44      use ioipsl_getincom
     5      use surfdat_h
     6      use comdiurn_h
     7      use comsaison_h
     8      use comsoil_h
     9      USE comgeomfi_h
    510
    611      implicit none
    … …  
    3338#include "paramet.h"
    3439#include "dimphys.h"
    35 #include "comgeomfi.h"
    36 #include "surfdat.h"
    37 #include "comsoil.h"
    38 #include "comdiurn.h"
    3940#include "callkeys.h"
    4041#include "comcstfi.h"
    4142#include "planete.h"
    42 #include "comsaison.h"
    4343#include "control.h"
    4444#include "comvert.h"
    4545#include "netcdf.inc"
    46 #include "watercap.h"
    47 #include "fisice.h"
    4846#include "logic.h"
    49 #include "advtrac.h"
    5047#include "comgeom.h"
    5148
    … …  
    9592      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
    9693      Logical  tracerdyn
    97       integer :: nq=1 ! number of tracers
     94      integer :: nq !=1 ! number of tracers
    9895 
    9996      character*2 str2
    … …  
    110107
    111108!     added by BC
    112       REAL cloudfrac(ngridmx,nlayermx)
    113       REAL hice(ngridmx),totcloudfrac(ngridmx)
     109      REAL cloudfrac(1,nlayermx)
     110      REAL hice(1),totcloudfrac(1)
    114111      REAL qzero1D   !initial water amount on the ground
     112
     113!     added by AS to avoid the use of adv trac common
     114      character*20 nametrac(nqmx)   ! name of the tracer (no need for adv trac common)
     115
    115116
    116117c=======================================================================
    117118c INITIALISATION
    118119c=======================================================================
     120
     121      !! those are defined in surfdat_h.F90
     122      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(1))
     123      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(1))
     124      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(1))
     125      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(1))
     126      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(1))
     127      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(1))
     128      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(1))
     129      IF (.not. ALLOCATED(dryness)) ALLOCATE(dryness(1))
     130      IF (.not. ALLOCATED(watercaptag)) ALLOCATE(watercaptag(1))
     131      !! those are defined in comdiurn_h.F90
     132      IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(1))
     133      IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(1))
     134      IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(1))
     135      IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(1))
     136      !! those are defined in comsoil_h.F90
     137      IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
     138      IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
     139      IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
     140      !! those are defined in comgeomfi_h
     141      IF (.not. ALLOCATED(lati)) ALLOCATE(lati(1))
     142      IF (.not. ALLOCATED(long)) ALLOCATE(long(1))
     143      IF (.not. ALLOCATED(area)) ALLOCATE(area(1))
     144
    119145
    120146      saveprofile=.false.
    … …  
    214240          endif
    215241       
    216           ! initialize advection schemes to Van-Leer for all tracers
    217           do iq=1,nq
    218             iadv(iq)=3 ! Van-Leer
    219           enddo
    220        
    221242          do iq=1,nq
    222243          ! minimal version, just read in the tracer names, 1 per line
    223             read(90,*,iostat=ierr) tnom(iq)
     244            read(90,*,iostat=ierr) nametrac(iq)
    224245            if (ierr.ne.0) then
    225246              write(*,*) 'rcm1d: error reading tracer names...'
    … …  
    232253         i_h2o_vap=0
    233254         do iq=1,nq
    234            if (tnom(iq)=="co2_ice") then
     255           if (nametrac(iq)=="co2_ice") then
    235256             i_co2_ice=iq
    236            elseif (tnom(iq)=="h2o_ice") then
     257           elseif (nametrac(iq)=="h2o_ice") then
    237258             i_h2o_ice=iq
    238            elseif (tnom(iq)=="h2o_vap") then
     259           elseif (nametrac(iq)=="h2o_vap") then
    239260             i_h2o_vap=iq
    240261           endif
    … …  
    258279        do iq=1,nq
    259280          write(str7,'(a1,i2.2)')'q',iq
    260           tnom(iq)=str7
     281          nametrac(iq)=str7
    261282        enddo
    262283        ! actually, we'll need at least one "co2_ice" tracer
    263284        ! (for surface CO2 ice)
    264         tnom(1)="co2_ice"
     285        nametrac(1)="co2_ice"
    265286        i_co2_ice=1
    266287      endif ! of if (tracer)
    … …  
    482503c    ---------------------------
    483504
    484       DO iq=1,nqmx
     505      DO iq=1,nq
    485506        DO ilayer=1,nlayer
    486507           q(ilayer,iq) = 0.
    … …  
    488509      ENDDO
    489510     
    490       DO iq=1,nqmx
     511      DO iq=1,nq
    491512        qsurf(iq) = 0.
    492513      ENDDO
    … …  
    683704c   les variables purement physiques a "physiq"...
    684705
    685       call physdem1("startfi.nc",long,lati,nsoilmx,nqmx,
     706      call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
    686707     &              dtphys,float(day0),
    687708     &              time,tsurf,tsoil,emis,q2,qsurf,
    688709     &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
    689      &              cloudfrac,totcloudfrac,hice)
     710     &              cloudfrac,totcloudfrac,hice,nametrac)
    690711
    691712c=======================================================================
    … …  
    740761
    741762
    742       CALL physiq (1,llm,nqmx,
     763      CALL physiq (1,llm,nq,
     764     .     nametrac,
    743765     ,     firstcall,lastcall,
    744766     ,     day,time,dtphys,
    … …  
    802824c       --------------------------------------
    803825
    804         DO iq = 1, nqmx
     826        DO iq = 1, nq
    805827          DO ilayer=1,nlayer
    806828             q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)

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

    r253 r787  
    1       subroutine soil(ngrid,nsoil,firstcall,
     1      subroutine soil(ngrid,nsoil,firstcall,lastcall,
    22     &          therm_i,
    33     &          timestep,tsurf,tsoil,
    44     &          capcal,fluxgrd)
     5
     6      use comsoil_h
     7
    58      implicit none
    69
    … …  
    1720#include "dimphys.h"
    1821
    19 #include"comsoil.h"
    2022
    2123c-----------------------------------------------------------------------
    … …  
    2628      integer nsoil     ! number of soil layers
    2729      logical firstcall ! identifier for initialization call
     30      logical lastcall
    2831      real therm_i(ngrid,nsoil) ! thermal inertia
    2932      real timestep         ! time step
    … …  
    3538
    3639! local saved variables:
    37       !real,save :: mthermdiff(ngridmx,0:nsoilmx-1) ! mid-layer thermal diffusivity
    38       !real,save :: thermdiff(ngridmx,nsoilmx-1) ! inter-layer thermal diffusivity
    39       !real,save :: coefq(0:nsoilmx-1)          ! q_{k+1/2} coefficients
    40       !real,save :: coefd(ngridmx,nsoilmx-1)    ! d_k coefficients
    41       !real,save :: alph(ngridmx,nsoilmx-1)     ! alpha_k coefficients
    42       !real,save :: beta(ngridmx,nsoilmx-1)     ! beta_k coefficients
    43       !real,save :: mu
    44 
    45       real mthermdiff(ngridmx,0:nsoilmx-1) ! mid-layer thermal diffusivity
    46       real thermdiff(ngridmx,nsoilmx-1)    ! inter-layer thermal diffusivity
    47       real coefq(0:nsoilmx-1)              ! q_{k+1/2} coefficients
    48       real coefd(ngridmx,nsoilmx-1)        ! d_k coefficients
    49       real alph(ngridmx,nsoilmx-1)         ! alpha_k coefficients
    50       real beta(ngridmx,nsoilmx-1)         ! beta_k coefficients
    51       real mu
     40      real,dimension(:,:),save,allocatable :: mthermdiff ! mid-layer thermal diffusivity
     41      real,dimension(:,:),save,allocatable :: thermdiff ! inter-layer thermal diffusivity
     42      real,dimension(:),save,allocatable :: coefq ! q_{k+1/2} coefficients
     43      real,dimension(:,:),save,allocatable :: coefd ! d_k coefficients
     44      real,dimension(:,:),save,allocatable :: alph ! alpha_k coefficients
     45      real,dimension(:,:),save,allocatable :: beta ! beta_k coefficients
     46      real,save :: mu
    5247           
    53       save mthermdiff,thermdiff,coefq,coefd,alph,beta,mu
    54 
    5548! local variables:
    5649      integer ig,ik
    … …  
    6457      ! note: firstcall is set to .true. or .false. by the caller
    6558      !       and not changed by soil.F
     59
     60      ALLOCATE(mthermdiff(ngrid,0:nsoilmx-1)) ! mid-layer thermal diffusivity
     61      ALLOCATE(thermdiff(ngrid,nsoilmx-1))    ! inter-layer thermal diffusivity
     62      ALLOCATE(coefq(0:nsoilmx-1))              ! q_{k+1/2} coefficients
     63      ALLOCATE(coefd(ngrid,nsoilmx-1))        ! d_k coefficients
     64      ALLOCATE(alph(ngrid,nsoilmx-1))         ! alpha_k coefficients
     65      ALLOCATE(beta(ngrid,nsoilmx-1))         ! beta_k coefficients
     66
    6667! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
    6768      do ig=1,ngrid
    … …  
    186187      enddo
    187188
     189      if (lastcall) then
     190        IF ( ALLOCATED(mthermdiff)) DEALLOCATE(mthermdiff)
     191        IF ( ALLOCATED(thermdiff)) DEALLOCATE(thermdiff)
     192        IF ( ALLOCATED(coefq)) DEALLOCATE(coefq)
     193        IF ( ALLOCATED(coefd)) DEALLOCATE(coefd)
     194        IF ( ALLOCATED(alph)) DEALLOCATE(alph)
     195        IF ( ALLOCATED(beta)) DEALLOCATE(beta)
     196      endif
    188197
    189198      end

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

    r786 r787  
    11      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
     2
     3      use comsoil_h
     4
    25      implicit none
    36
    … …  
    2629#include "dimphys.h"
    2730
    28 #include "comsoil.h"
    2931#include "netcdf.inc"
    3032!======================================================================
    … …  
    3537      integer ngrid     ! # of horizontal grid points
    3638      integer nsoil     ! # of soil layers
     39      integer sta       ! # at which reading starts
    3740      real tsurf(ngrid) ! surface temperature
    3841!  output:
    39       real tsoil(ngridmx,nsoilmx)       ! soil temperature
     42      real tsoil(ngrid,nsoilmx) ! soil temperature
    4043
    4144!======================================================================
    … …  
    6770      real,parameter :: default_volcapa=1.e6
    6871
     72      integer, dimension(2) :: sta2d
     73      integer, dimension(2) :: siz2d
     74
    6975!======================================================================
     76
     77      !! this is defined in comsoil_h
     78      IF (.not.ALLOCATED(layer))
     79     . ALLOCATE(layer(nsoilmx))
     80      IF (.not.ALLOCATED(mlayer))
     81     . ALLOCATE(mlayer(0:nsoilmx-1))
     82      IF (.not.ALLOCATED(inertiedat))
     83     . ALLOCATE(inertiedat(ngrid,nsoilmx))
     84
     85      !! this is defined in dimphys.h
     86      sta = cursor
    7087
    7188! 1. Depth coordinate
    … …  
    94111          interpol=.true.
    95112        endif
     113
     114      sta2d = (/sta,1/)
     115      siz2d = (/ngrid,dimlen/)
     116
    96117
    97118! 1.2 Find out the # of dimensions <inertiedat> was defined as using
    … …  
    231252       allocate(surfinertia(ngrid))
    232253#ifdef NC_DOUBLE
    233        ierr = NF_GET_VAR_DOUBLE(nid, nvarid, surfinertia)
    234 #else
    235        ierr = NF_GET_VAR_REAL(nid, nvarid, surfinertia)
     254       ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, surfinertia)
     255#else
     256       ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, surfinertia)
    236257#endif
    237258        if (ierr.NE.NF_NOERR) then
    … …  
    258279         endif
    259280#ifdef NC_DOUBLE
    260         ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldinertiedat)
    261 #else
    262         ierr = NF_GET_VAR_REAL(nid,nvarid,oldinertiedat)
     281        ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldinertiedat)
     282#else
     283        ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldinertiedat)
    263284#endif
    264285        if (ierr.NE.NF_NOERR) then
    … …  
    268289       else ! put values in therm_i
    269290#ifdef NC_DOUBLE
    270         ierr = NF_GET_VAR_DOUBLE(nid,nvarid,inertiedat)
    271 #else
    272         ierr = NF_GET_VAR_REAL(nid,nvarid,inertiedat)
     291        ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,inertiedat)
     292#else
     293        ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,inertiedat)
    273294#endif
    274295        if (ierr.NE.NF_NOERR) then
    … …  
    300321         endif
    301322#ifdef NC_DOUBLE
    302         ierr=NF_GET_VAR_DOUBLE(nid,nvarid,oldtsoil)
    303 #else
    304         ierr=NF_GET_VAR_REAL(nid,nvarid,oldtsoil)
     323        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldtsoil)
     324#else
     325        ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldtsoil)
    305326#endif
    306327        if (ierr.ne.NF_NOERR) then
    … …  
    310331       else ! put values in tsoil
    311332#ifdef NC_DOUBLE
    312         ierr=NF_GET_VAR_DOUBLE(nid,nvarid,tsoil)
    313 #else
    314         ierr=NF_GET_VAR_REAL(nid,nvarid,tsoil)
     333        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,tsoil)
     334#else
     335        ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,tsoil)
    315336#endif
    316337        if (ierr.ne.NF_NOERR) then

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

    r253 r787  
    1       SUBROUTINE surface_nature(obliquit,qsurf,qsurfliquid
     1      SUBROUTINE surface_nature(ngrid,nq,obliquit,qsurf,qsurfliquid
    22     &   ,qsurfsnow,rnat,oceanarea)
     3
     4      USE surfdat_h
     5      USE comsoil_h
     6      USE comgeomfi_h
     7      USE tracer_h
     8
    39      IMPLICIT none
    410
    … …  
    3137#include "comcstfi.h"
    3238#include "callkeys.h"
    33 #include "tracer.h"
    34 #include "fisice.h"
    35 #include "comgeomfi.h"
    36 #include "surfdat.h"
    37 #include "comsoil.h"
    3839
    39         REAL qsurf(ngridmx,nqmx),ps(ngridmx)
    40         REAL qsurfliquid(ngridmx)
    41         REAL qsurfsnow(ngridmx)
     40        integer ngrid,nq
     41
     42        REAL qsurf(ngrid,nq),ps(ngrid)
     43        REAL qsurfliquid(ngrid)
     44        REAL qsurfsnow(ngrid)
    4245        INTEGER iq, ig
    43         INTEGER rnat(ngridmx)
     46        INTEGER rnat(ngrid)
    4447        REAL oceanarea
    4548        REAL obliquit
    4649
    47         do ig=1,ngridmx
     50        do ig=1,ngrid
    4851           rnat(ig)=1
    4952           dryness(ig)=1        !(coefficient for evaporation)
    … …  
    5659       
    5760        oceanarea=0.
    58         do ig=1,ngridmx
     61        do ig=1,ngrid
    5962           if (rnat(ig).eq.0)then
    6063              oceanarea=oceanarea+area(ig)

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

    r135 r787  
    1       SUBROUTINE surfini(ngrid,qsurf,psolaralb)
     1      SUBROUTINE surfini(ngrid,nq,qsurf,psolaralb)
     2
     3      USE surfdat_h
     4      USE tracer_h
     5
    26      IMPLICIT NONE
    37c=======================================================================
    … …  
    1115#include "dimensions.h"
    1216#include "dimphys.h"
    13 #include "surfdat.h"
    1417#include "callkeys.h"
    15 #include "tracer.h"
    1618c
    1719      INTEGER ngrid,ig,icap
     20      INTEGER nq
    1821      REAL  piceco2(ngrid),psolaralb(ngrid)
    19       REAL qsurf(ngrid,nqmx) !tracer on surface (kg/m2)
     22      REAL qsurf(ngrid,nq) !tracer on surface (kg/m2)
    2023
    2124      EXTERNAL ISMIN,ISMAX
    … …  
    4245!      DO ig=1,ngrid
    4346c        IF (water) THEN
    44 c          if (qsurf(ig,nqmx).gt.0.005) then
     47c          if (qsurf(ig,nq).gt.0.005) then
    4548c             psolaralb(ig,1) = 0.4
    4649c             psolaralb(ig,2) = 0.4

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

    r726 r787  
    11c=======================================================================
    2       SUBROUTINE tabfi(nid,Lmodif,tab0,day_ini,lmax,p_rad,
     2      SUBROUTINE tabfi(ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
    33     .                 p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    44c=======================================================================
    … …  
    4646      use ioipsl_getincom
    4747
     48      USE surfdat_h
     49      use comsoil_h
     50      USE comgeomfi_h
     51
    4852      implicit none
    4953 
    … …  
    5155#include "dimphys.h"
    5256#include "comcstfi.h"
    53 #include "comgeomfi.h"
    5457#include "planete.h"
    55 #include "surfdat.h"
    56 #include "comsoil.h"
    5758#include "netcdf.inc"
    5859#include "callkeys.h"
    … …  
    6162c   Declarations
    6263c-----------------------------------------------------------------------
     64
     65      INTEGER ngrid
    6366
    6467c Arguments
    … …  
    104107c  Initialization of some physical constants
    105108c informations on physics grid
    106       if(ngridmx.ne.tab_cntrl(tab0+1)) then
    107          print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngridmx'
    108          print*,tab_cntrl(tab0+1),ngridmx
     109      if(ngrid.ne.tab_cntrl(tab0+1)) then
     110         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
     111         print*,tab_cntrl(tab0+1),ngrid
    109112      endif
    110113      lmax = nint(tab_cntrl(tab0+2))
    … …  
    168171      write(*,*) 'Reading tab_cntrl when calling tabfi before changes'
    169172      write(*,*) '*****************************************************'
    170       write(*,5) '(1)      = ngridmx?',tab_cntrl(tab0+1),float(ngridmx)
     173      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
    171174      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
    172175      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)
    … …  
    490493      write(*,*) 'Reading tab_cntrl when calling tabfi AFTER changes'
    491494      write(*,*) '*****************************************************'
    492       write(*,5) '(1)      = ngridmx?',tab_cntrl(tab0+1),float(ngridmx)
     495      write(*,5) '(1)        = ngrid?',tab_cntrl(tab0+1),float(ngrid)
    493496      write(*,5) '(2)            lmax',tab_cntrl(tab0+2),float(lmax)
    494497      write(*,5) '(3)         day_ini',tab_cntrl(tab0+3),float(day_ini)

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

    r728 r787  
    1       subroutine totalcloudfrac(rneb,totalrneb,pplev,pq,tau)
     1      subroutine totalcloudfrac(ngrid,nq,rneb,totalrneb,pplev,pq,tau)
    22
    33      use watercommon_h
     4      use comdiurn_h
     5      USE comgeomfi_h
     6      USE tracer_h
    47      implicit none
    58
    … …  
    1922#include "dimphys.h"
    2023#include "comcstfi.h"
    21 #include "tracer.h"
    22 #include "fisice.h"
    23 #include "comgeomfi.h"
    24 #include "comdiurn.h"
    2524#include "callkeys.h"
    2625
     26      integer ngrid, nq
    2727
    28       real rneb(ngridmx,nlayermx)             ! cloud fraction     
    29       real pplev(ngridmx,nlayermx+1)
    30       real pq(ngridmx,nlayermx,nqmx)
    31       real tau(ngridmx,nlayermx)
     28      real rneb(ngrid,nlayermx)             ! cloud fraction     
     29      real pplev(ngrid,nlayermx+1)
     30      real pq(ngrid,nlayermx,nq)
     31      real tau(ngrid,nlayermx)
    3232
    33       real totalrneb(ngridmx)                 ! total cloud fraction
     33      real totalrneb(ngrid)                 ! total cloud fraction
    3434      real, dimension(nlayermx+1) :: masse
    3535      integer, parameter          :: recovery=7
    … …  
    5050
    5151
    52       do ig=1,ngridmx
     52      do ig=1,ngrid
    5353         totalrneb(ig) = 0.
    5454

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

    r736 r787  
    99      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol,Psat_water
    1010      use radcommon_h, only : sigma
     11      USE surfdat_h
     12      USE comgeomfi_h
     13      USE tracer_h
    1114
    1215      implicit none
    … …  
    4346#include "comcstfi.h"
    4447#include "callkeys.h"
    45 #include "surfdat.h"
    46 #include "comgeomfi.h"
    47 #include "tracer.h"
    48 
    49 #include "watercap.h"
    5048
    5149!     arguments
    … …  
    8482      integer ilev,ig,ilay,nlev
    8583
    86       REAL z4st,zdplanck(ngridmx)
    87       REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
    88       REAL zcdv(ngridmx),zcdh(ngridmx)
    89       REAL zcdv_true(ngridmx),zcdh_true(ngridmx)
    90       REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
    91       REAL zh(ngridmx,nlayermx),zt(ngridmx,nlayermx)
    92       REAL ztsrf(ngridmx)
    93       REAL z1(ngridmx),z2(ngridmx)
    94       REAL zmass(ngridmx,nlayermx)
    95       REAL zfluxv(ngridmx,nlayermx),zfluxt(ngridmx,nlayermx),zfluxq(ngridmx,nlayermx)
    96       REAL zb0(ngridmx,nlayermx)
    97       REAL zExner(ngridmx,nlayermx),zovExner(ngridmx,nlayermx)
    98       REAL zcv(ngridmx,nlayermx),zdv(ngridmx,nlayermx)  !inversion coefficient for winds
    99       REAL zct(ngridmx,nlayermx),zdt(ngridmx,nlayermx)  !inversion coefficient for temperature
    100       REAL zcq(ngridmx,nlayermx),zdq(ngridmx,nlayermx)  !inversion coefficient for tracers
     84      REAL z4st,zdplanck(ngrid)
     85      REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
     86      REAL zcdv(ngrid),zcdh(ngrid)
     87      REAL zcdv_true(ngrid),zcdh_true(ngrid)
     88      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
     89      REAL zh(ngrid,nlayermx),zt(ngrid,nlayermx)
     90      REAL ztsrf(ngrid)
     91      REAL z1(ngrid),z2(ngrid)
     92      REAL zmass(ngrid,nlayermx)
     93      REAL zfluxv(ngrid,nlayermx),zfluxt(ngrid,nlayermx),zfluxq(ngrid,nlayermx)
     94      REAL zb0(ngrid,nlayermx)
     95      REAL zExner(ngrid,nlayermx),zovExner(ngrid,nlayermx)
     96      REAL zcv(ngrid,nlayermx),zdv(ngrid,nlayermx)  !inversion coefficient for winds
     97      REAL zct(ngrid,nlayermx),zdt(ngrid,nlayermx)  !inversion coefficient for temperature
     98      REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx)  !inversion coefficient for tracers
    10199      REAL zcst1
    102100      REAL zu2!, a
    103       REAL zcq0(ngridmx),zdq0(ngridmx)
    104       REAL zx_alf1(ngridmx),zx_alf2(ngridmx)
     101      REAL zcq0(ngrid),zdq0(ngrid)
     102      REAL zx_alf1(ngrid),zx_alf2(ngrid)
    105103
    106104      LOGICAL firstcall
    … …  
    111109!     -------
    112110      INTEGER iq
    113       REAL zq(ngridmx,nlayermx,nqmx)
    114       REAL zqnoevap(ngridmx,nlayermx) !special for water case to compute where evaporated water goes.
    115       REAL pdqevap(ngridmx,nlayermx) !special for water case to compute where evaporated water goes.
    116       REAL zdmassevap(ngridmx)
    117       REAL rho(ngridmx)         ! near-surface air density
    118       REAL qsat(ngridmx),psat(ngridmx)
     111      REAL zq(ngrid,nlayermx,nq)
     112      REAL zqnoevap(ngrid,nlayermx) !special for water case to compute where evaporated water goes.
     113      REAL pdqevap(ngrid,nlayermx) !special for water case to compute where evaporated water goes.
     114      REAL zdmassevap(ngrid)
     115      REAL rho(ngrid)         ! near-surface air density
     116      REAL qsat(ngrid),psat(ngrid)
    119117      DATA firstcall/.true./
    120118      REAL kmixmin
    … …  
    122120!     Variables added for implicit latent heat inclusion
    123121!     --------------------------------------------------
    124       real dqsat(ngridmx), qsat_temp1, qsat_temp2,psat_temp
     122      real dqsat(ngrid), qsat_temp1, qsat_temp2,psat_temp
    125123
    126124      integer, save :: ivap, iliq, iliq_surf,iice_surf ! also make liq for clarity on surface...
    … …  
    241239!     ------------------------------------------------------
    242240
    243       call vdif_kc(ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
     241      call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay,pu,pv,zh,zcdv_true,pq2,zkv,zkh) !JL12 why not call vdif_kc with updated winds and temperature
    244242
    245243!     Adding eddy mixing to mimic 3D general circulation in 1D
    … …  
    262260 
    263261!JL12 calculate the flux coefficients (tables multiplied elements by elements)
    264       zfluxv(1:ngridmx,1:nlayermx)=zkv(1:ngridmx,1:nlayermx)*zb0(1:ngridmx,1:nlayermx)
     262      zfluxv(1:ngrid,1:nlayermx)=zkv(1:ngrid,1:nlayermx)*zb0(1:ngrid,1:nlayermx)
    265263     
    266264!-----------------------------------------------------------------------
    … …  
    358356!     JL12 calculate the flux coefficients (tables multiplied elements by elements)
    359357!     ---------------------------------------------------------------
    360       zfluxq(1:ngridmx,1:nlayermx)=zkh(1:ngridmx,1:nlayermx)*zb0(1:ngridmx,1:nlayermx) !JL12 we save zfluxq which doesn't need the Exner factor
    361       zfluxt(1:ngridmx,1:nlayermx)=zfluxq(1:ngridmx,1:nlayermx)*zExner(1:ngridmx,1:nlayermx)
    362 
     358      zfluxq(1:ngrid,1:nlayermx)=zkh(1:ngrid,1:nlayermx)*zb0(1:ngrid,1:nlayermx) !JL12 we save zfluxq which doesn't need the Exner factor
     359      zfluxt(1:ngrid,1:nlayermx)=zfluxq(1:ngrid,1:nlayermx)*zExner(1:ngrid,1:nlayermx)
    363360
    364361      DO ig=1,ngrid
    … …  
    404401
    405402         DO ig=1,ngrid
    406 
    407403            z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zfluxt(ig,1)*zct(ig,1)*zovExner(ig,1) &
    408404                + pfluxsrf(ig)*ptimestep + zdplanck(ig)*ptsrf(ig)
    … …  
    434430!     Calculate vertical flux from the bottom to the first layer (dust)
    435431!     -----------------------------------------------------------------
    436          do ig=1,ngridmx 
     432         do ig=1,ngrid
    437433            rho(ig) = zb0(ig,1) /ptimestep
    438434         end do
    … …  
    470466                  ENDDO
    471467               else             ! general case
    472                   DO ig=1,ngrid
     468                  do ig=1,ngrid
    473469                     z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
    474470                     zcq(ig,1)=(zmass(ig,1)*zq(ig,1,iq)+zfluxq(ig,2)*zcq(ig,2)+(-pdqsdif(ig,iq))*ptimestep)*z1(ig)
    … …  
    500496           
    501497               ! compute evaporation efficiency
    502                do ig = 1, ngrid
     498               do ig=1,ngrid
    503499                  if(rnat(ig).eq.1)then
    504500                     dryness(ig)=pqsurf(ig,iliq_surf)+pqsurf(ig,iice_surf)
    … …  
    509505
    510506               do ig=1,ngrid
    511 
    512507                ! Calculate the value of qsat at the surface (water)
    513508                call Psat_water(ptsrf(ig),pplev(ig,1),psat(ig),qsat(ig))
    … …  
    526521                  zdq(ig,nlay)=zfluxq(ig,nlay)*z1(ig)
    527522               enddo
    528            
     523         
    529524               do ilay=nlay-1,2,-1
    530525                  do ig=1,ngrid
    … …  
    658653            ENDDO
    659654
    660             DO ig=1,ngrid
     655            do ig=1,ngrid
    661656               z1(ig)=1./(zmass(ig,1)+zfluxq(ig,2)*(1.-zdq(ig,2)))
    662657               zcq(ig,1)=(zmass(ig,1)*zqnoevap(ig,1)+zfluxq(ig,2)*zcq(ig,2))*z1(ig)
    … …  
    692687      enddo
    693688     
    694       DO ig=1,ngrid  ! computing sensible heat flux (atm => surface)
     689      DO ig=1,ngrid ! computing sensible heat flux (atm => surface)
    695690         sensibFlux(ig)=cpp*zfluxt(ig,1)/ptimestep*(zt(ig,1)*zovExner(ig,1)-ztsrf(ig))
    696691      ENDDO
    … …  
    710705               enddo
    711706            enddo
    712             do ig=1,ngrid
     707            do ig=1,ngrid
    713708               zdmassevap(ig)=SUM(pdqevap(ig,:)*zmass(ig,:))*ptimestep
    714709            end do         

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

    r135 r787  
    1       SUBROUTINE vdif_kc(dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
     1      SUBROUTINE vdif_kc(ngrid,dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
    22      IMPLICIT NONE
    33c.......................................................................
    … …  
    2727c
    2828c.......................................................................
     29      INTEGER ngrid
    2930      REAL dt,g
    30       REAL zlev(ngridmx,nlayermx+1)
    31       REAL zlay(ngridmx,nlayermx)
    32       REAL u(ngridmx,nlayermx)
    33       REAL v(ngridmx,nlayermx)
    34       REAL teta(ngridmx,nlayermx)
    35       REAL cd(ngridmx)
    36       REAL q2(ngridmx,nlayermx+1)
    37       REAL km(ngridmx,nlayermx+1)
    38       REAL kn(ngridmx,nlayermx+1)
     31      REAL zlev(ngrid,nlayermx+1)
     32      REAL zlay(ngrid,nlayermx)
     33      REAL u(ngrid,nlayermx)
     34      REAL v(ngrid,nlayermx)
     35      REAL teta(ngrid,nlayermx)
     36      REAL cd(ngrid)
     37      REAL q2(ngrid,nlayermx+1)
     38      REAL km(ngrid,nlayermx+1)
     39      REAL kn(ngrid,nlayermx+1)
    3940c.......................................................................
    4041c
    … …  
    4950c
    5051c.......................................................................
    51       INTEGER nlay,nlev,ngrid
    52       REAL unsdz(ngridmx,nlayermx)
    53       REAL unsdzdec(ngridmx,nlayermx+1)
    54       REAL q(ngridmx,nlayermx+1)
     52      INTEGER nlay,nlev
     53      REAL unsdz(ngrid,nlayermx)
     54      REAL unsdzdec(ngrid,nlayermx+1)
     55      REAL q(ngrid,nlayermx+1)
    5556c.......................................................................
    5657c
    … …  
    6263c
    6364c.......................................................................
    64       REAL kmpre(ngridmx,nlayermx+1)
     65      REAL kmpre(ngrid,nlayermx+1)
    6566      REAL qcstat
    6667      REAL q2cstat
    … …  
    7071c
    7172c.......................................................................
    72       REAL long(ngridmx,nlayermx+1)
     73      REAL long(ngrid,nlayermx+1)
    7374c.......................................................................
    7475c
    … …  
    9394      REAL mcstat
    9495      REAL m2cstat
    95       REAL m(ngridmx,nlayermx+1)
    96       REAL mpre(ngridmx,nlayermx+1)
    97       REAL m2(ngridmx,nlayermx+1)
    98       REAL n2(ngridmx,nlayermx+1)
     96      REAL m(ngrid,nlayermx+1)
     97      REAL mpre(ngrid,nlayermx+1)
     98      REAL m2(ngrid,nlayermx+1)
     99      REAL n2(ngrid,nlayermx+1)
    99100c.......................................................................
    100101c
    … …  
    118119      LOGICAL gnsup
    119120      REAL gm
    120 c      REAL ri(ngridmx,nlayermx+1)
    121       REAL sn(ngridmx,nlayermx+1)
    122       REAL snq2(ngridmx,nlayermx+1)
    123       REAL sm(ngridmx,nlayermx+1)
    124       REAL smq2(ngridmx,nlayermx+1)
     121c      REAL ri(ngrid,nlayermx+1)
     122      REAL sn(ngrid,nlayermx+1)
     123      REAL snq2(ngrid,nlayermx+1)
     124      REAL sm(ngrid,nlayermx+1)
     125      REAL smq2(ngrid,nlayermx+1)
    125126c.......................................................................
    126127c
    … …  
    183184      PARAMETER (nlay=nlayermx)
    184185      PARAMETER (nlev=nlayermx+1)
    185       PARAMETER (ngrid=ngridmx)
    186186c
    187187      PARAMETER (
    … …  
    209209c
    210210      DO ilev=1,nlev
    211                                                       DO igrid=1,ngrid
     211                                                 DO igrid=1,ngrid
    212212        q2(igrid,ilev)=amax1(q2(igrid,ilev),q2min)
    213213        q(igrid,ilev)=sqrt(q2(igrid,ilev))
    214                                                       ENDDO
     214                                                 ENDDO
    215215      ENDDO
    216216c
    217                                                       DO igrid=1,ngrid
     217                                                 DO igrid=1,ngrid
    218218      tmp1=cd(igrid)*(u(igrid,1)**2+v(igrid,1)**2)
    219219      q2(igrid,1)=b1**(2.E+0/3.E+0)*tmp1
    220220      q2(igrid,1)=amax1(q2(igrid,1),q2min)
    221221      q(igrid,1)=sqrt(q2(igrid,1))
    222                                                       ENDDO
     222                                                 ENDDO
    223223c
    224224c.......................................................................
    … …  
    237237c
    238238      DO ilay=1,nlay
    239                                                       DO igrid=1,ngrid
     239                                                 DO igrid=1,ngrid
    240240        unsdz(igrid,ilay)=1.E+0/(zlev(igrid,ilay+1)-zlev(igrid,ilay))
    241                                                       ENDDO
     241                                                 ENDDO
    242242      ENDDO
    243                                                       DO igrid=1,ngrid
     243                                                 DO igrid=1,ngrid
    244244      unsdzdec(igrid,1)=1.E+0/(zlay(igrid,1)-zlev(igrid,1))
    245                                                       ENDDO
     245                                                 ENDDO
    246246      DO ilay=2,nlay
    247                                                       DO igrid=1,ngrid
     247                                                 DO igrid=1,ngrid
    248248        unsdzdec(igrid,ilay)=1.E+0/(zlay(igrid,ilay)-zlay(igrid,ilay-1))
    249                                                       ENDDO
     249                                                 ENDDO
    250250      ENDDO
    251                                                       DO igrid=1,ngrid
     251                                                 DO igrid=1,ngrid
    252252      unsdzdec(igrid,nlay+1)=1.E+0/(zlev(igrid,nlay+1)-zlay(igrid,nlay))
    253                                                       ENDDO
     253                                                 ENDDO
    254254c
    255255c.......................................................................
    … …  
    257257c.......................................................................
    258258c
    259                                                       DO igrid=1,ngrid
     259                                                 DO igrid=1,ngrid
    260260      m2(igrid,1)=(unsdzdec(igrid,1)
    261261     &                   *u(igrid,1))**2
    … …  
    264264      m(igrid,1)=sqrt(m2(igrid,1))
    265265      mpre(igrid,1)=m(igrid,1)
    266                                                       ENDDO
     266                                                 ENDDO
    267267c
    268268c-----------------------------------------------------------------------
    269269      DO ilev=2,nlev-1
    270                                                       DO igrid=1,ngrid
     270                                                 DO igrid=1,ngrid
    271271c-----------------------------------------------------------------------
    272272c
    … …  
    295295c
    296296c-----------------------------------------------------------------------
    297                                                       ENDDO
     297                                                 ENDDO
    298298      ENDDO
    299299c-----------------------------------------------------------------------
    300300c
    301                                                       DO igrid=1,ngrid
     301                                                 DO igrid=1,ngrid
    302302      m2(igrid,nlev)=m2(igrid,nlev-1)
    303303      m(igrid,nlev)=m(igrid,nlev-1)
    304304      mpre(igrid,nlev)=m(igrid,nlev)
    305                                                       ENDDO
     305                                                 ENDDO
    306306c
    307307c.......................................................................
    … …  
    311311c-----------------------------------------------------------------------
    312312      DO ilev=2,nlev-1
    313                                                       DO igrid=1,ngrid
     313                                                 DO igrid=1,ngrid
    314314c-----------------------------------------------------------------------
    315315c
    … …  
    375375c
    376376c-----------------------------------------------------------------------
    377                                                       ENDDO
     377                                                  ENDDO
    378378      ENDDO
    379379c-----------------------------------------------------------------------
    … …  
    383383c.......................................................................
    384384c
    385                                                       DO igrid=1,ngrid
     385                                                 DO igrid=1,ngrid
    386386      kn(igrid,1)=knmin
    387387      km(igrid,1)=kmmin
    388388      kmpre(igrid,1)=km(igrid,1)
    389                                                       ENDDO
     389                                                 ENDDO
    390390c
    391391c-----------------------------------------------------------------------
    392392      DO ilev=2,nlev-1
    393                                                       DO igrid=1,ngrid
     393                                                 DO igrid=1,ngrid
    394394c-----------------------------------------------------------------------
    395395c
    … …  
    401401c
    402402c-----------------------------------------------------------------------
    403                                                       ENDDO
     403                                                 ENDDO
    404404      ENDDO
    405405c-----------------------------------------------------------------------
    406406c
    407                                                       DO igrid=1,ngrid
     407                                                 DO igrid=1,ngrid
    408408      kn(igrid,nlev)=kn(igrid,nlev-1)
    409409      km(igrid,nlev)=km(igrid,nlev-1)
    410410      kmpre(igrid,nlev)=km(igrid,nlev)
    411                                                       ENDDO
     411                                                 ENDDO
    412412c
    413413c.......................................................................
    … …  
    418418      DO 10001 ilev=2,nlev-1
    419419c---->
    420       DO 10002 igrid=1,ngrid 
     420      DO 10002 igrid=1,ngrid
    421421c
    422422c.......................................................................
    … …  
    539539c
    540540c
    541                                                       DO igrid=1,ngrid
     541                                                 DO igrid=1,ngrid
    542542      kn(igrid,1)=knmin
    543543      km(igrid,1)=kmmin
    … …  
    546546      kn(igrid,nlev)=kn(igrid,nlev-1)
    547547      km(igrid,nlev)=km(igrid,nlev-1)
    548                                                       ENDDO
     548                                                 ENDDO
    549549
    550550c

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

    r650 r787  
    99      use watercommon_h, only : RLVTT, T_h2O_ice_liq, RCPD, mx_eau_sol
    1010      use radcommon_h, only : sigma
     11      USE surfdat_h
     12      USE comgeomfi_h
     13      USE tracer_h
    1114
    1215      implicit none
    … …  
    3942#include "comcstfi.h"
    4043#include "callkeys.h"
    41 #include "surfdat.h"
    42 #include "comgeomfi.h"
    43 #include "tracer.h"
    44 
    45 #include "watercap.h"
    4644
    4745!     arguments
    … …  
    7876      integer ilev,ig,ilay,nlev
    7977
    80       REAL z4st,zdplanck(ngridmx)
    81       REAL zkv(ngridmx,nlayermx+1),zkh(ngridmx,nlayermx+1)
    82       REAL zcdv(ngridmx),zcdh(ngridmx)
    83       REAL zcdv_true(ngridmx),zcdh_true(ngridmx)
    84       REAL zu(ngridmx,nlayermx),zv(ngridmx,nlayermx)
    85       REAL zh(ngridmx,nlayermx)
    86       REAL ztsrf2(ngridmx)
    87       REAL z1(ngridmx),z2(ngridmx)
    88       REAL za(ngridmx,nlayermx),zb(ngridmx,nlayermx)
    89       REAL zb0(ngridmx,nlayermx)
    90       REAL zc(ngridmx,nlayermx),zd(ngridmx,nlayermx)
     78      REAL z4st,zdplanck(ngrid)
     79      REAL zkv(ngrid,nlayermx+1),zkh(ngrid,nlayermx+1)
     80      REAL zcdv(ngrid),zcdh(ngrid)
     81      REAL zcdv_true(ngrid),zcdh_true(ngrid)
     82      REAL zu(ngrid,nlayermx),zv(ngrid,nlayermx)
     83      REAL zh(ngrid,nlayermx)
     84      REAL ztsrf2(ngrid)
     85      REAL z1(ngrid),z2(ngrid)
     86      REAL za(ngrid,nlayermx),zb(ngrid,nlayermx)
     87      REAL zb0(ngrid,nlayermx)
     88      REAL zc(ngrid,nlayermx),zd(ngrid,nlayermx)
    9189      REAL zcst1
    9290      REAL zu2!, a
    93       REAL zcq(ngridmx,nlayermx),zdq(ngridmx,nlayermx)
    94       REAL evap(ngridmx)
    95       REAL zcq0(ngridmx),zdq0(ngridmx)
    96       REAL zx_alf1(ngridmx),zx_alf2(ngridmx)
     91      REAL zcq(ngrid,nlayermx),zdq(ngrid,nlayermx)
     92      REAL evap(ngrid)
     93      REAL zcq0(ngrid),zdq0(ngrid)
     94      REAL zx_alf1(ngrid),zx_alf2(ngrid)
    9795
    9896      LOGICAL firstcall
    … …  
    103101!     variables added for CO2 condensation
    104102!     ------------------------------------
    105       REAL hh                   !, zhcond(ngridmx,nlayermx)
     103      REAL hh                   !, zhcond(ngrid,nlayermx)
    106104!     REAL latcond,tcond1mb
    107105!     REAL acond,bcond
    … …  
    112110!     -------
    113111      INTEGER iq
    114       REAL zq(ngridmx,nlayermx,nqmx)
    115       REAL zq1temp(ngridmx)
    116       REAL rho(ngridmx)         ! near-surface air density
    117       REAL qsat(ngridmx)
     112      REAL zq(ngrid,nlayermx,nq)
     113      REAL zq1temp(ngrid)
     114      REAL rho(ngrid)         ! near-surface air density
     115      REAL qsat(ngrid)
    118116      DATA firstcall/.true./
    119117      REAL kmixmin
    … …  
    121119!     Variables added for implicit latent heat inclusion
    122120!     --------------------------------------------------
    123       real latconst, dqsat(ngridmx), qsat_temp1, qsat_temp2
    124       real z1_Tdry(ngridmx), z2_Tdry(ngridmx)
    125       real z1_T(ngridmx), z2_T(ngridmx)
    126       real zb_T(ngridmx)
    127       real zc_T(ngridmx,nlayermx)
    128       real zd_T(ngridmx,nlayermx)
    129       real lat1(ngridmx), lat2(ngridmx)
     121      real latconst, dqsat(ngrid), qsat_temp1, qsat_temp2
     122      real z1_Tdry(ngrid), z2_Tdry(ngrid)
     123      real z1_T(ngrid), z2_T(ngrid)
     124      real zb_T(ngrid)
     125      real zc_T(ngrid,nlayermx)
     126      real zd_T(ngrid,nlayermx)
     127      real lat1(ngrid), lat2(ngrid)
    130128      real surfh2otot
    131129      logical surffluxdiag
    … …  
    152150
    153151      IF (firstcall) THEN
    154 !         IF(ngrid.NE.ngridmx) THEN
    155 !            PRINT*,'STOP dans vdifc'
    156 !            PRINT*,'probleme de dimensions :'
    157 !            PRINT*,'ngrid  =',ngrid
    158 !            PRINT*,'ngridmx  =',ngridmx
    159 !            STOP
    160 !         ENDIF
    161152!     To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal)
    162153!     bcond=1./tcond1mb
    … …  
    254245!     ------------------------------------------------------
    255246
    256       call vdif_kc(ptimestep,g,pzlev,pzlay
     247      call vdif_kc(ngrid,ptimestep,g,pzlev,pzlay
    257248     &     ,pu,pv,ph,zcdv_true
    258249     &     ,pq2,zkv,zkh)
    … …  
    442433!     Calculate vertical flux from the bottom to the first layer (dust)
    443434!     -----------------------------------------------------------------
    444          do ig=1,ngridmx 
     435         do ig=1,ngrid 
    445436            rho(ig) = zb0(ig,1) /ptimestep
    446437         end do

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

    r135 r787  
    3535c
    3636
    37       real dzq(ngridmx,llm),dzqw(ngridmx,llm),adzqw(ngridmx,llm),dzqmax
     37      real dzq(ngrid,llm),dzqw(ngrid,llm),adzqw(ngrid,llm),dzqmax
    3838      real newmasse
    3939      real sigw, Mtot, MQtot
    … …  
    9090
    9191       do l = 1,llm          ! loop different than when w<0
    92         do  ij = 1,ngrid
     92        do ij=1,ngrid
    9393
    9494         if(w(ij,l).gt.0.)then
    … …  
    133133
    134134c      Surface flux up:
    135        do  ij = 1,ngrid
     135       do ij=1,ngrid
    136136         if(w(ij,1).lt.0.) wq(ij,1)=0. ! warning : not always valid
    137137       end do
    138138
    139139       do l = 1,llm-1  ! loop different than when w>0
    140         do  ij = 1,ngrid
     140        do ij=1,ngrid
    141141         if(w(ij,l+1).le.0)then
    142142

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

    r526 r787  
    3939!
    4040!=================================================================
     41
     42      USE surfdat_h
    4143 
    4244      implicit none
    … …  
    5254#include "netcdf.inc"
    5355#include "temps.h"
    54 #include "surfdat.h"
    5556
    5657! Arguments on input:
    … …  
    9596      character(len=20),save :: nom_def(n_nom_def_max)
    9697      logical,save :: firstcall=.true.
    97      
    98 #ifndef MESOSCALE
     98
    9999!***************************************************************
    100100!Sortie des variables au rythme voulu
    … …  
    108108!***************************************************************
    109109
     110#ifndef NOWRITEDIAGFI
     111
    110112! At very first call, check if there is a "diagfi.def" to use and read it
    111113! -----------------------------------------------------------------------
    … …  
    113115         firstcall=.false.
    114116
    115   !      Open diagfi.def definition file if there is one:
     117!      Open diagfi.def definition file if there is one:
    116118         open(99,file="diagfi.def",status='old',form='formatted',
    117119     s   iostat=ierr2)
    … …  
    192194      endif ! if (firstnom.eq.'1234567890')
    193195
    194       if (ngridmx.eq.1) then
     196      if (ngrid.eq.1) then
    195197        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
    196198        ! are undefined; so set them to 1
    … …  
    448450
    449451#endif
     452
    450453      end

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

    r135 r787  
    1010! (yielding the sought time series of the variable)
    1111
     12use comsoil_h
     13
    1214implicit none
    1315
    … …  
    1618#include"paramet.h"
    1719#include"control.h"
    18 #include"comsoil.h"
    1920#include"netcdf.inc"
    2021
    … …  
    7879 
    7980  ! Define dimensions and axis attributes
    80   call iniwritesoil(nid)
     81  call iniwritesoil(nid,ngrid)
    8182 
    8283  ! set zitau to -1 to be compatible with zitau incrementation step below

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

    r533 r787  
    4343! Addition by RW (2010) to allow OLR to be saved in .nc format
    4444      use radinc_h, only : L_NSPECTI
     45      USE surfdat_h
    4546
    4647      implicit none
    … …  
    5657#include "netcdf.inc"
    5758#include "temps.h"
    58 #include "surfdat.h"
    5959#include "callkeys.h"
    6060

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

    r533 r787  
    4343! Addition by RW (2010) to allow OSR to be saved in .nc format
    4444      use radinc_h, only : L_NSPECTV
     45      USE surfdat_h
    4546
    4647      implicit none
    … …  
    5657#include "netcdf.inc"
    5758#include "temps.h"
    58 #include "surfdat.h"
    5959#include "callkeys.h"
    6060
Note: See TracChangeset for help on using the changeset viewer.