Changeset 1529 for trunk/LMDZ.GENERIC/libf/phystd

Timestamp:

Apr 5, 2016, 10:51:51 AM (9 years ago)

Author:

emillour

Message:

Generic GCM: Towards a cleaner separation between physics and dynamics

• Got rid of references to "dimensions.h" from physics packages: use nbp_lon (=iim), nbp_lat (=jjp1) and nbp_lev from module mod_grid_phy_lmdz (in phy_common) instead.
• Removed module "comhdiff_mod.F90", as it is only used by module surf_heat_transp_mod.F90, moved module variables there.
• Addedin "surf_heat_transp_mod" local versions of some arrays and routines (from dyn3d) required to compute gradient, divergence, etc. on the global dynamics grid. As before, the slab ocean only works in serial.

EM

Location:

trunk/LMDZ.GENERIC/libf/phystd

Files:

• aeroptproperties.F90 (modified) (1 diff)
• callcorrk.F90 (modified) (13 diffs)
• comhdiff_mod.F90 (deleted)
• dyn1d/diverg.F (deleted)
• dyn1d/gr_dyn_fi.F (deleted)
• dyn1d/gr_fi_dyn.F (deleted)
• dyn1d/gr_u_scal.F (deleted)
• dyn1d/gr_v_scal.F (deleted)
• dyn1d/grad.F (deleted)
• dyn1d/rcm1d.F (modified) (1 diff)
• inifis_mod.F90 (modified) (5 diffs)
• inistats.F (modified) (4 diffs)
• initracer.F (modified) (1 diff)
• iniwrite.F (modified) (13 diffs)
• iniwrite_specIR.F (modified) (12 diffs)
• iniwrite_specVI.F (modified) (11 diffs)
• iniwritesoil.F90 (modified) (8 diffs)
• mass_redistribution.F90 (modified) (6 diffs)
• mkstat.F90 (modified) (4 diffs)
• physiq.F90 (modified) (2 diffs)
• radcommon_h.F90 (modified) (3 diffs)
• radii_mod.F90 (modified) (8 diffs)
• radinc_h.F90 (modified) (3 diffs)
• soil.F (modified) (1 diff)
• suaer_corrk.F90 (modified) (1 diff)
• surf_heat_transp_mod.F90 (modified) (4 diffs)
• vdif_kc.F (modified) (1 diff)
• vlz_fi.F (modified) (1 diff)
• writediagfi.F (modified) (17 diffs)
• writediagsoil.F90 (modified) (12 diffs)
• writediagspecIR.F (modified) (11 diffs)
• writediagspecVI.F (modified) (11 diffs)
• wstats.F90 (modified) (14 diffs)

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

@@ -31 +31 @@
 !     Varying nueff section removed by R. Wordsworth for simplicity
 !     ==============================================================
-
-!#include "dimensions.h"
-!#include "dimphys.h"
 
 !     Local variables 

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

@@ -19 +19 @@
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval,    &
-                              kastprof,strictboundcorrk,specOLR,CLFvarying
+      use callkeys_mod, only: varactive,diurnal,tracer,water,nosurf,varfixed,satval,        &
+                              kastprof,strictboundcorrk,specOLR,CLFvarying
 
       implicit none
@@ -144 +144 @@
 
       ! Local aerosol optical properties for each column on RADIATIVE grid.
-      real*8,save ::  QXVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8,save ::  QSVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8,save ::  GVAER(L_LEVELS+1,L_NSPECTV,naerkind)
-      real*8,save ::  QXIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-      real*8,save ::  QSIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-      real*8,save ::  GIAER(L_LEVELS+1,L_NSPECTI,naerkind)
-
-      !REAL :: QREFvis3d(ngrid,nlayer,naerkind)
-      !REAL :: QREFir3d(ngrid,nlayer,naerkind)
-      !save QREFvis3d, QREFir3d 
+      real*8,save,allocatable ::  QXVAER(:,:,:)
+      real*8,save,allocatable ::  QSVAER(:,:,:)
+      real*8,save,allocatable ::  GVAER(:,:,:)
+      real*8,save,allocatable ::  QXIAER(:,:,:)
+      real*8,save,allocatable ::  QSIAER(:,:,:)
+      real*8,save,allocatable ::  GIAER(:,:,:)
+
       real, dimension(:,:,:), save, allocatable :: QREFvis3d
       real, dimension(:,:,:), save, allocatable :: QREFir3d
@@ -188 +185 @@
 
 
-      qxvaer(:,:,:)=0.0
-      qsvaer(:,:,:)=0.0
-      gvaer(:,:,:) =0.0
-
-      qxiaer(:,:,:)=0.0
-      qsiaer(:,:,:)=0.0
-      giaer(:,:,:) =0.0
-
       if(firstcall) then
+
+        ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
+        if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS+1,L_NSPECTV,naerkind))
+        if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS+1,L_NSPECTV,naerkind))
+        if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS+1,L_NSPECTV,naerkind))
+        if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS+1,L_NSPECTI,naerkind))
+        if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS+1,L_NSPECTI,naerkind))
+        if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS+1,L_NSPECTI,naerkind))
 
          !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
@@ -212 +209 @@
          endif
          call su_aer_radii(ngrid,nlayer,reffrad,nueffrad)
-         
+         
          
 !--------------------------------------------------
@@ -264 +261 @@
       end if ! of if (firstcall)
 
-
 !=======================================================================
 !          I.b  Initialization on every call   
 !=======================================================================
  
- 
+      qxvaer(:,:,:)=0.0
+      qsvaer(:,:,:)=0.0
+      gvaer(:,:,:) =0.0
+
+      qxiaer(:,:,:)=0.0
+      qsiaer(:,:,:)=0.0
+      giaer(:,:,:) =0.0
+
 !--------------------------------------------------
 !     Effective radius and variance of the aerosols
@@ -277 +280 @@
 
          if ((iaer.eq.iaero_co2).and.tracer.and.(igcm_co2_ice.gt.0)) then ! Treat condensed co2 particles.
-            call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
+            call co2_reffrad(ngrid,nlayer,nq,pq,reffrad(1,1,iaero_co2))
             print*,'Max. CO2 ice particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
             print*,'Min. CO2 ice particle size = ',minval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
-         end if
+         end if
          
          if ((iaer.eq.iaero_h2o).and.water) then ! Treat condensed water particles. To be generalized for other aerosols ...
-            call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
+            call h2o_reffrad(ngrid,nlayer,pq(1,1,igcm_h2o_ice),pt, &
                              reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
             print*,'Max. H2O cloud particle size = ',maxval(reffrad(1:ngrid,1:nlayer,iaer))/1.e-6,' um'
@@ -290 +293 @@
          
          if(iaer.eq.iaero_dust)then
-            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
+            call dust_reffrad(ngrid,nlayer,reffrad(1,1,iaero_dust))
             print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          
          if(iaer.eq.iaero_h2so4)then
-            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
+            call h2so4_reffrad(ngrid,nlayer,reffrad(1,1,iaero_h2so4))
             print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
          endif
          
           if(iaer.eq.iaero_back2lay)then
-            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
+            call back2lay_reffrad(ngrid,reffrad(1,1,iaero_back2lay),nlayer,pplev)
          endif
          
@@ -321 +324 @@
            reffrad,QREFvis3d,QREFir3d,                             & 
            tau_col,cloudfrac,totcloudfrac,clearsky)                
-          
+          
 
 
@@ -475 +478 @@
          DO nw=1,L_NSPECTV ! Short Wave loop.
             albv(nw)=albedo(ig,nw)
-         ENDDO
+         ENDDO
 
          if (nosurf) then ! Case with no surface.
             DO nw=1,L_NSPECTV
                if(albv(nw).gt.0.0) then
-                  print*,'For open lower boundary in callcorrk must'
+                  print*,'For open lower boundary in callcorrk must'
                   print*,'have spectral surface band albedos all set to zero!'
                   call abort
-               endif
-            ENDDO         
+               endif
+            ENDDO         
          endif
 
@@ -770 +773 @@
 
          Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
-         glat_ig=glat(ig)
+         glat_ig=glat(ig)
 
 !-----------------------------------------------------------------------
@@ -798 +801 @@
          else ! During the night, fluxes = 0.
             nfluxtopv       = 0.0d0
-            fluxtopvdn      = 0.0d0
+            fluxtopvdn      = 0.0d0
             nfluxoutv_nu(:) = 0.0d0
             nfluxgndv_nu(:) = 0.0d0
@@ -813 +816 @@
             surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV))      
             if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive.
-               DO nw=1,L_NSPECTV                  
-                  albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
+               DO nw=1,L_NSPECTV                  
+                  albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw)
                ENDDO
-               albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
+               albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux
                albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV))
             else
                albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
             endif
-         else
-            albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
-         endif
+         else
+            albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0.
+         endif
 
 
@@ -956 +959 @@
         IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
 !$OMP END MASTER
-!$OMP BARRIER   
+!$OMP BARRIER
         IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad)
         IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad)

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

@@ -907 +907 @@
       end                       !rcm1d
  
-c***********************************************************************
-c***********************************************************************
-c     Subroutines Bidons utilise seulement en 3D, mais
-c     necessaire a la compilation de rcm1d en 1D
-
-!      subroutine gr_fi_dyn
-!      RETURN
-!      END
- 
-c***********************************************************************
-c***********************************************************************
-
-!#include "../dyn3d/disvert.F"
-!#include "../dyn3d/abort_gcm.F"
-!#include "../dyn3d/diverg.F"
-!#include "../dyn3d/grad.F"
-!#include "../dyn3d/gr_u_scal.F"
-!#include "../dyn3d/gr_v_scal.F"
-!#include "../dyn3d/gr_dyn_fi.F"
-
+

trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

@@ -9 +9 @@
              prad,pg,pr,pcpp)
 
-  use radinc_h, only : naerkind
+  use radinc_h, only: ini_radinc_h, naerkind
+  use radcommon_h, only: ini_radcommon_h
   use datafile_mod, only: datadir
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
@@ -21 +22 @@
   use planetwide_mod, only: planetwide_sumval
   use callkeys_mod
+  use mod_phys_lmdz_para, only : is_parallel
 
 !=======================================================================
@@ -74 +76 @@
   REAL SSUM
  
-  CHARACTER ch1*12
-  CHARACTER ch80*80
-
-  logical chem, h2o
-
-  real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
-
   ! initialize constants in comcstfi_mod
   rad=prad
@@ -333 +328 @@
      call getin_p("ok_slab_ocean",ok_slab_ocean)
      write(*,*) "ok_slab_ocean = ",ok_slab_ocean
+     ! Sanity check: for now slab oncean only works in serial mode
+     if (ok_slab_ocean.and.is_parallel) then
+       write(*,*) " Error: slab ocean should only be used in serial mode!"
+       call abort
+     endif
 
      write(*,*) "Use slab-sea-ice ?"
@@ -718 +718 @@
   ENDDO
 
+  ! initialize variables in radinc_h
+  call ini_radinc_h(nlayer)
+ 
+  ! allocate "radcommon_h" arrays
+  call ini_radcommon_h()
+
   ! allocate "comsoil_h" arrays
   call ini_comsoil_h(ngrid)

trunk/LMDZ.GENERIC/libf/phystd/inistats.F

@@ -1 +1 @@
       subroutine inistats(ierr)
 
+      use statto_mod, only: istats,istime
       use mod_phys_lmdz_para, only : is_master
-      use statto_mod, only: istats,istime
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff,pseudoalt,presnivs
-      USE comconst_mod, ONLY: daysec,dtphys,pi
+      USE comconst_mod, ONLY: pi
+      USE time_phylmdz_mod, ONLY: daysec,dtphys
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       implicit none
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
       integer,intent(out) :: ierr
       integer :: nid
       integer :: l,nsteppd
-      real, dimension(llm) ::  sig_s
+      real, dimension(nbp_lev) ::  sig_s
+      real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
       integer :: idim_lat,idim_lon,idim_llm,idim_llmp1,idim_time
       real, dimension(istime) :: lt
@@ -39 +40 @@
       write (*,*) 
 
-      do l= 1, llm
+      do l= 1, nbp_lev
          sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
          pseudoalt(l)=-10.*log(presnivs(l)/preff)   
       enddo
+      
+      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+      !add extra redundant point (180 degrees, since lon_reg starts at -180
+      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
 
       if (is_master) then
@@ -53 +58 @@
       endif
 
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_lat)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_lon)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
-      ierr = NF_DEF_DIM (nid, "llmp1", llm+1, idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_lat)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_lon)
+      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
+      ierr = NF_DEF_DIM (nid, "llmp1", nbp_lev+1, idim_llmp1)
       ierr = NF_DEF_DIM (nid, "Time", NF_UNLIMITED, idim_time)
 
@@ -68 +73 @@
      &            "degrees_north",1,idim_lat,nvarid,ierr)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
 #else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
 #endif
       call def_var_stats(nid,"longitude","East longitude",
      &            "degrees_east",1,idim_lon,nvarid,ierr)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
 #else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 

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

@@ -22 +22 @@
 c            Ehouarn Millour (oct. 2008) identify tracers by their names
 c=======================================================================
-
-
-#include "dimensions.h"
 
       integer :: ngrid,nq

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

@@ -1 @@
-      SUBROUTINE iniwrite(nid,idayref,phis)
+      SUBROUTINE iniwrite(nid,idayref,phis,area)
 
       use comsoil_h, only: mlayer, nsoilmx
-      use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
-      use time_phylmdz_mod, only: daysec, dtphys
+      USE comcstfi_mod, only: g, mugaz, omeg, rad, rcp, pi 
       USE comvert_mod, ONLY: ap,bp,aps,bps,pseudoalt
       USE logic_mod, ONLY: fxyhypb,ysinus
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
+      USE time_phylmdz_mod, ONLY: daysec, dtphys
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
       IMPLICIT NONE
 
@@ -26 +28 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 c   Arguments:
@@ -36 +35 @@
       integer,intent(in) :: nid        ! NetCDF file ID
       INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
-      real,intent(in) :: phis(ip1jmp1) ! surface geopotential
+      real,intent(in) :: phis(nbp_lon+1,nbp_lat) ! surface geopotential
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -44 +44 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
-
-      integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
-      integer :: idim_rlatu,idim_rlatv,idim_llmp1,idim_llm
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
+
+      integer :: nvarid,idim_index,idim_rlonv
+      integer :: idim_rlatu,idim_llmp1,idim_llm
       integer :: idim_nsoilmx ! "subsurface_layers" dimension ID #
       integer, dimension(2) :: id  
@@ -54 +55 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = real(iim)
-      tab_cntrl(2)  = real(jjm)
-      tab_cntrl(3)  = real(llm)
+      tab_cntrl(1)  = real(nbp_lon)
+      tab_cntrl(2)  = real(nbp_lat-1)
+      tab_cntrl(3)  = real(nbp_lev)
       tab_cntrl(4)  = real(idayref)
       tab_cntrl(5)  = rad
@@ -99 +100 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
-      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
-      ierr = NF_DEF_DIM (nid, "interlayer", (llm+1), idim_llmp1)
-      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
+!      ierr = NF_DEF_DIM (nid, "rlonu", iip1, idim_rlonu)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
+!      ierr = NF_DEF_DIM (nid, "rlatv", jjm, idim_rlatv)
+      ierr = NF_DEF_DIM (nid, "interlayer", (nbp_lev+1), idim_llmp1)
+      ierr = NF_DEF_DIM (nid, "altitude", nbp_lev, idim_llm)
       ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
 c
@@ -129 +130 @@
 c --------------------------
 c  longitudes and latitudes
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
-     .                       "Longitudes at u nodes")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
-#endif
+!
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "rlonu", NF_DOUBLE, 1, idim_rlonu,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "rlonu", NF_FLOAT, 1, idim_rlonu,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 21,
+!     .                       "Longitudes at u nodes")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonu/pi*180)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonu/pi*180)
+!#endif
 c
 c --------------------------
@@ -156 +158 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+c
+c --------------------------
+      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+      !add extra redundant point (180 degrees, since lon_reg starts at -180
+      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -173 +179 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -199 +205 @@
 c
 c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
-     .                       "Latitudes at v nodes")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
-#endif
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "rlatv", NF_DOUBLE, 1, idim_rlatv,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "rlatv", NF_FLOAT, 1, idim_rlatv,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
+!     .                       "Latitudes at v nodes")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatv/pi*180)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatv/pi*180)
+!#endif
 c
 c --------------------------
@@ -274 +280 @@
 c  Mesh area and conversion coefficients cov. <-> contra. <--> natural
 
-      id(1)=idim_rlonu
-      id(2)=idim_rlatu
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
-#endif
-c
-      id(1)=idim_rlonv
-      id(2)=idim_rlatv
-c
-c --------------------------
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
-     .             "Conversion coefficients cov <--> natural")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
-#endif
+!      id(1)=idim_rlonu
+!      id(2)=idim_rlatu
+c
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "cu", NF_DOUBLE, 2, id,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "cu", NF_FLOAT, 2, id,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+!     .             "Conversion coefficients cov <--> natural")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cu)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,cu)
+!#endif
+c
+!      id(1)=idim_rlonv
+!      id(2)=idim_rlatv
+c
+c --------------------------
+!      ierr = NF_REDEF (nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_DEF_VAR (nid, "cv", NF_DOUBLE, 2, id,nvarid)
+!#else
+!      ierr = NF_DEF_VAR (nid, "cv", NF_FLOAT, 2, id,nvarid)
+!#endif
+!      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 40,
+!     .             "Conversion coefficients cov <--> natural")
+!      ierr = NF_ENDDEF(nid)
+!#ifdef NC_DOUBLE
+!      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cv)
+!#else
+!      ierr = NF_PUT_VAR_REAL (nid,nvarid,cv)
+!#endif
 c
       id(1)=idim_rlonv
@@ -325 +331 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 c
@@ -350 +356 @@
 c
 
-      write(*,*)'iniwrite: iim,jjm,llm,idayref',iim,jjm,llm,idayref
+      write(*,*)'iniwrite: nbp_lon,nbp_lat,nbp_lev,idayref',
+     & nbp_lon,nbp_lat,nbp_lev,idayref
       write(*,*)'iniwrite: rad,omeg,g,mugaz,rcp',
-     s rad,omeg,g,mugaz,rcp
+     & rad,omeg,g,mugaz,rcp
       write(*,*)'iniwrite: daysec,dtphys',daysec,dtphys
 

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

@@ -1 @@
-      SUBROUTINE iniwrite_specIR(nid,idayref)
+      SUBROUTINE iniwrite_specIR(nid,idayref,area)
 
       use radinc_h, only: L_NSPECTI
@@ -8 +8 @@
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
 
       implicit none
@@ -27 +29 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
-!#include"dimphys.h"
+      include "netcdf.inc"
 
 c   Arguments:
@@ -38 +36 @@
       integer,intent(in) :: nid        ! NetCDF file ID
       INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -45 +44 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
 
       integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
@@ -57 +57 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = FLOAT(iim)
-      tab_cntrl(2)  = FLOAT(jjm)
-      tab_cntrl(3)  = FLOAT(llm)
+      tab_cntrl(1)  = FLOAT(nbp_lon)
+      tab_cntrl(2)  = FLOAT(nbp_lat-1)
+      tab_cntrl(3)  = FLOAT(nbp_lev)
       tab_cntrl(4)  = FLOAT(idayref)
       tab_cntrl(5)  = rad
@@ -102 +102 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
       ierr = NF_DEF_DIM (nid, "IR Wavenumber",L_NSPECTI,idim_bandsIR)
 
@@ -141 +141 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+c
+c --------------------------
+      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+      !add extra redundant point (180 degrees, since lon_reg starts at -180
+      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+      
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -158 +162 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -172 +176 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_DOUBLE,1,
-     .                          idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"IR Wavenumber",NF_FLOAT,1,
-     .                          idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 34,
@@ -196 +200 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"IR Bandwidth",NF_DOUBLE,1,
-     .                          idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"IR Bandwidth",NF_FLOAT,1,
-     .                          idim_bandsIR,nvarid)
+     .                                idim_bandsIR,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 25,
@@ -230 +234 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 

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

@@ -1 @@
-      SUBROUTINE iniwrite_specVI(nid,idayref)
+      SUBROUTINE iniwrite_specVI(nid,idayref,area)
 
       use radinc_h, only: L_NSPECTV
       use radcommon_h, only: WNOV,DWNV
-      use comsoil_h
       use comcstfi_mod, only: rad, omeg, g, mugaz, rcp, pi
       use time_phylmdz_mod, only: daysec, dtphys
@@ -9 +8 @@
       USE serre_mod, ONLY: clon,clat,grossismx,grossismy,dzoomx,dzoomy
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
+      USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
 
       implicit none
@@ -28 +29 @@
 c   -------------
 
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom.h"
-#include "netcdf.inc"
-!#include"dimphys.h"
+      include "netcdf.inc"
 
 c   Arguments:
 c   ----------
 
-      integer nid        ! NetCDF file ID
-      INTEGER*4 idayref  ! date (initial date for this run)
+      integer,intent(in) :: nid        ! NetCDF file ID
+      INTEGER*4,INTENT(IN) :: idayref  ! date (initial date for this run)
+      real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 c   Local:
@@ -46 +44 @@
       REAL tab_cntrl(length) ! run parameters are stored in this array
       INTEGER ierr
+      REAl :: lon_reg_ext(nbp_lon+1) ! extended longitudes
 
       integer :: nvarid,idim_index,idim_rlonu,idim_rlonv
@@ -58 +57 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = FLOAT(iim)
-      tab_cntrl(2)  = FLOAT(jjm)
-      tab_cntrl(3)  = FLOAT(llm)
+      tab_cntrl(1)  = FLOAT(nbp_lon)
+      tab_cntrl(2)  = FLOAT(nbp_lat-1)
+      tab_cntrl(3)  = FLOAT(nbp_lev)
       tab_cntrl(4)  = FLOAT(idayref)
       tab_cntrl(5)  = rad
@@ -103 +102 @@
 
       ierr = NF_DEF_DIM (nid, "index", length, idim_index)
-      ierr = NF_DEF_DIM (nid, "latitude", jjp1, idim_rlatu)
-      ierr = NF_DEF_DIM (nid, "longitude", iip1, idim_rlonv)
+      ierr = NF_DEF_DIM (nid, "latitude", nbp_lat, idim_rlatu)
+      ierr = NF_DEF_DIM (nid, "longitude", nbp_lon+1, idim_rlonv)
       ierr = NF_DEF_DIM (nid, "VI Wavenumber",L_NSPECTV,idim_bandsVI)
 
@@ -141 +140 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
-#endif
-c
-c --------------------------
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lat_reg/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lat_reg/pi*180)
+#endif
+c
+c --------------------------
+      lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+      !add extra redundant point (180 degrees, since lon_reg starts at -180
+      lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+
       ierr = NF_REDEF (nid)
 #ifdef NC_DOUBLE
@@ -158 +161 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlonv/pi*180)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lon_reg_ext/pi*180)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,lon_reg_ext/pi*180)
 #endif
 c
@@ -172 +175 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_DOUBLE,1,
-     .                          idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"VI Wavenumber",NF_FLOAT,1,
-     .                          idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 33,
@@ -197 +200 @@
 #ifdef NC_DOUBLE
       ierr=NF_DEF_VAR(nid,"VI Bandwidth",NF_DOUBLE,1,
-     .                          idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #else
       ierr=NF_DEF_VAR(nid,"VI Bandwidth",NF_FLOAT,1,
-     .                          idim_bandsVI,nvarid)
+     .                                idim_bandsVI,nvarid)
 #endif
       ierr=NF_PUT_ATT_TEXT (nid,nvarid,"long_name", 24,
@@ -231 +234 @@
       ierr = NF_ENDDEF(nid)
 #ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aire)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,aire)
+      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
+#else
+      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
 #endif
 

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

@@ -1 @@
-subroutine iniwritesoil(nid,ngrid)
+subroutine iniwritesoil(nid,ngrid,inertia,area)
 
 ! initialization routine for 'writediagoil'. Here we create/define
@@ -5 +5 @@
 ! (time-independent) parameters.
 
-use comsoil_h, only: mlayer, inertiedat, nsoilmx
-use comcstfi_mod, only: pi
+use comsoil_h, only: mlayer, nsoilmx
+USE comcstfi_mod, only: pi
+USE regular_lonlat_mod, ONLY: lon_reg, lat_reg
+use mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
 
 implicit none
 
-#include"dimensions.h"
-#include"paramet.h"
-#include"comgeom.h"
-#include"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
 integer,intent(in) :: ngrid
 integer,intent(in) :: nid ! NetCDF output file ID
+real,intent(in) :: inertia(nbp_lon+1,nbp_lat,nsoilmx)
+real,intent(in) :: area(nbp_lon+1,nbp_lat) ! mesh area (m2)
 
 ! Local variables:
@@ -30 +31 @@
 integer,dimension(3) :: dimids ! to store IDs of dimensions of a variable
 character(len=60) :: text ! to store some text
-real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
+real,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
 integer :: i,j,l,ig0
+real :: lon_reg_ext(nbp_lon+1) ! extended longitudes
 
 ! 1. Define the dimensions
@@ -38 +40 @@
 
 ! Define the dimensions
-ierr=NF_DEF_DIM(nid,"longitude",iip1,idim_rlonv)
-! iip1 known from paramet.h
+ierr=NF_DEF_DIM(nid,"longitude",nbp_lon+1,idim_rlonv)
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define longitude dimension"
 endif
-ierr=NF_DEF_DIM(nid,"latitude",jjp1,idim_rlatu)
-! jjp1 known from paramet.h
+ierr=NF_DEF_DIM(nid,"latitude",nbp_lat,idim_rlatu)
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define latitude dimension"
 endif
 ierr=NF_DEF_DIM(nid,"depth",nsoilmx,idim_depth)
-! nsoilmx known from dimphys.h
+! nsoilmx known from comsoil_h
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not define depth dimension"
@@ -81 +81 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
+lon_reg_ext(1:nbp_lon)=lon_reg(1:nbp_lon)
+!add extra redundant point (180 degrees, since lon_reg starts at -180
+lon_reg_ext(nbp_lon+1)=-lon_reg_ext(1)
+
 ! Write longitude to file
 ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,rlonv*(180./pi))
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,rlonv*(180./pi))
-#endif
-! Note: rlonv is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,lon_reg_ext*(180./pi))
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,lon_reg_ext*(180./pi))
+#endif
+! Note: rlonv is known from comgeom.h and pi from comcstfi.h
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write longitude variable"
@@ -117 +121 @@
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,rlatu*(180./pi))
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,rlatu*(180./pi))
-#endif
-! Note: rlatu is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,lat_reg*(180./pi))
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,lat_reg*(180./pi))
+#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write longitude variable"
@@ -209 +212 @@
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,aire)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,aire)
-#endif
-! Note: aire is known from comgeom.h
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,area)
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,area)
+#endif
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write area variable"
@@ -240 +242 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
-! Recast data along 'dynamics' grid
-! Note: inertiedat is known from comsoil_h
-
-do l=1,nsoilmx
-  ! handle the poles
-  do i=1,iip1
-    data3(i,1,l)=inertiedat(1,l)
-    data3(i,jjp1,l)=inertiedat(ngrid,l)
-  enddo
-  ! rest of the grid
-  do j=2,jjm
-    ig0=1+(j-2)*iim
-    do i=1,iim
-      data3(i,j,l)=inertiedat(ig0+i,l)
-    enddo
-    data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude
-  enddo
-enddo ! of do l=1,nsoilmx
-
 ! Write data2 to file
 ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
 ! Write
 #ifdef NC_DOUBLE
-ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
-#else
-ierr=NF_PUT_VAR_REAL(nid,varid,data3)
+ierr=NF_PUT_VAR_DOUBLE(nid,varid,inertia)
+#else
+ierr=NF_PUT_VAR_REAL(nid,varid,inertia)
 #endif
 if (ierr.ne.NF_NOERR) then

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

@@ -1 @@
-      SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
+SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep,   &
                        rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs,     &
                        pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr,  &
@@ -224 +224 @@
 !        Van Leer scheme:
          w(1:nlayer+1)=-zmflux(1:nlayer+1)*ptimestep
-         call vl1d(zzt,2.,zzmass,w,ztm) 
-         call vl1d(zzu ,2.,zzmass,w,zum) 
-         call vl1d(zzv ,2.,zzmass,w,zvm) 
+         call vl1d(nlayer,zzt,2.,zzmass,w,ztm) 
+         call vl1d(nlayer,zzu,2.,zzmass,w,zum) 
+         call vl1d(nlayer,zzv,2.,zzmass,w,zvm) 
          do iq=1,nq
            zq1(1:nlayer)=zzq(1:nlayer,iq)
@@ -232 +232 @@
 !               print*,iq
 !               print*,zq1
-           call vl1d(zq1,2.,zzmass,w,zqm1)
+           call vl1d(nlayer,zq1,2.,zzmass,w,zqm1)
            do l=2,nlayer
               zzq(l,iq)=zq1(l)
@@ -281 +281 @@
       END DO  ! loop on ig 
 
-      return
-      end
-
-
+CONTAINS
 
 ! *****************************************************************
-      SUBROUTINE vl1d(q,pente_max,zzmass,w,qm)
+      SUBROUTINE vl1d(llm,q,pente_max,zzmass,w,qm)
 !
 !    
@@ -299 +296 @@
 !
 !   --------------------------------------------------------------------
+      
       IMPLICIT NONE
-
-#include "dimensions.h"
 
 !   Arguments:
 !   ----------
+      integer,intent(in) :: llm
       real zzmass(llm),pente_max
       REAL q(llm),qm(llm+1)
@@ -406 +403 @@
 !         end if
 
-      return
-      end
+       END SUBROUTINE vl1d
+
+END SUBROUTINE mass_redistribution

trunk/LMDZ.GENERIC/libf/phystd/mkstat.F90

@@ -10 +10 @@
 !  Yann W. july 2003
 
+use statto_mod, only: istime,count
 use mod_phys_lmdz_para, only : is_master
-use statto_mod, only: istime,count
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
 
 implicit none
 
-#include "dimensions.h"
-#include "netcdf.inc"
+include "netcdf.inc"
 
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
 integer :: ierr,nid,nbvar,i,ndims,lt,nvarid
 integer, dimension(4) :: id,varid,start,size
 integer, dimension(5) :: dimids
 character (len=50) :: name,nameout,units,title
-real, dimension(iip1,jjp1,llm) :: sum3d,square3d,mean3d,sd3d
-real, dimension(iip1,jjp1) :: sum2d,square2d,mean2d,sd2d
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: sum3d,square3d,mean3d,sd3d
+real, dimension(nbp_lon+1,nbp_lat) :: sum2d,square2d,mean2d,sd2d
 real, dimension(istime) :: time
-real, dimension(jjp1) :: lat
-real, dimension(iip1) :: lon
-real, dimension(llm) :: alt
+real, dimension(nbp_lat) :: lat
+real, dimension(nbp_lon+1) :: lon
+real, dimension(nbp_lev) :: alt
 logical :: lcopy=.true.
 !integer :: latid,lonid,altid,timeid
@@ -35 +33 @@
 !integer, dimension(4) :: dimout
 
+! Incrementation of count for the last step, which is not done in wstats
+count(istime)=count(istime)+1
+
 if (is_master) then
 ! only the master needs do this
-
-! Incrementation of count for the last step, which is not done in wstats
-count(istime)=count(istime)+1
 
 ierr = NF_OPEN("stats.nc",NF_WRITE,nid)
@@ -107 +105 @@
 !      dimout(4)=timeid
 
-      size=(/iip1,jjp1,llm,1/)
+      size=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
       do lt=1,istime
          start=(/1,1,1,lt/)
@@ -137 +135 @@
 !      dimout(3)=timeid
 
-      size=(/iip1,jjp1,1,0/)
+      size=(/nbp_lon+1,nbp_lat,1,0/)
       do lt=1,istime
          start=(/1,1,lt,0/)

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

@@ -26 +26 @@
       use time_phylmdz_mod, only: ecritphy, iphysiq, nday
       use phyredem, only: physdem0, physdem1
-      use slab_ice_h
+      use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
+                            noceanmx
       use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
+                                ini_surf_heat_transp_mod, &
                                 ocean_slab_get_vars,ocean_slab_final
       use surf_heat_transp_mod,only :init_masquv
@@ -558 +560 @@
                                 sea_ice, pctsrf_sic)
 
+           call ini_surf_heat_transp_mod()
+           
            knindex(:) = 0
 

trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90

@@ -1 @@
-      module radcommon_h
-
-      use radinc_h
+module radcommon_h
+      use radinc_h, only: L_NSPECTI, L_NSPECTV, L_NGAUSS, NTstar, NTstop, &
+                          naerkind, nsizemax
       implicit none
 
@@ -121 +121 @@
       REAL :: omegaREFir(naerkind,nsizemax)
 
-      REAL tstellar ! Stellar brightness temperature (SW)
+      REAL,SAVE :: tstellar ! Stellar brightness temperature (SW)
 
-      real*8 planckir(L_NSPECTI,NTstop-NTstar+1)
+      real*8,save :: planckir(L_NSPECTI,NTstop-NTstar+1)
 
-      real*8 PTOP, TAUREF(L_LEVELS+1)
+      real*8,save :: PTOP
+      real*8,save,allocatable :: TAUREF(:)
 
-      real*8, parameter :: UBARI = 0.5D0
+      real*8,parameter :: UBARI = 0.5D0
 
-      real*8 gweight(L_NGAUSS)
+      real*8,save :: gweight(L_NGAUSS)
 !$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFvis,omegaREFir,&     ! gweight read by master in sugas_corrk
                 !$OMP tstellar,planckir,PTOP,TAUREF)
@@ -144 +145 @@
       real*8, parameter :: grav = 6.672E-11
 
-      real*8 Cmk
-      save Cmk
-      real*8 glat_ig
-      save glat_ig
+      real*8,save :: Cmk
+      real*8,save :: glat_ig
 !$OMP THREADPRIVATE(Cmk,glat_ig)
 
       ! extinction of incoming sunlight (Saturn's rings, eclipses, etc...)
-      REAL, DIMENSION(:), ALLOCATABLE :: eclipse
+      REAL, DIMENSION(:), ALLOCATABLE ,SAVE :: eclipse
 
       !Latitude-dependent gravity
-      REAL, DIMENSION(:), ALLOCATABLE :: glat
+      REAL, DIMENSION(:), ALLOCATABLE , SAVE :: glat
 !$OMP THREADPRIVATE(glat,eclipse)
 
-      end module radcommon_h
+contains
+
+      subroutine ini_radcommon_h
+      use radinc_h, only: L_LEVELS
+      implicit none
+      
+        allocate(TAUREF(L_LEVELS+1))
+      
+      end subroutine ini_radcommon_h
+
+end module radcommon_h

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

@@ -8 +8 @@
 !     water cloud optical properties
 
-      use callkeys_mod, only: radfixed,Nmix_co2,                        &
-                pres_bottom_tropo,pres_top_tropo,size_tropo,    &
-                pres_bottom_strato,size_strato
+      use callkeys_mod, only: radfixed,Nmix_co2,                    &
+                pres_bottom_tropo,pres_top_tropo,size_tropo,        &
+                pres_bottom_strato,size_strato
       
       real, save ::  rad_h2o
@@ -21 +21 @@
 
 
-      contains
+contains
 
 
@@ -38 +38 @@
       use ioipsl_getin_p_mod, only: getin_p
       use radinc_h, only: naerkind
-      use aerosol_mod
-!      USE tracer_h
-      Implicit none
-
-!      include "dimensions.h"
-!      include "dimphys.h"
+      use aerosol_mod, only: iaero_back2lay, iaero_co2, iaero_dust, &
+                             iaero_h2o, iaero_h2so4
+      Implicit none
 
       integer,intent(in) :: ngrid
@@ -81 +78 @@
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
             endif
-            
-            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
+            
+            if(iaer.eq.iaero_back2lay)then ! Two-layer aerosols
                reffrad(1:ngrid,1:nlayer,iaer) = 2.e-6
                nueffrad(1:ngrid,1:nlayer,iaer) = 0.1 
@@ -101 +98 @@
          if (radfixed) then
 
-            write(*,*)"radius of H2O water particles:"
+            write(*,*)"radius of H2O water particles:"
             rad_h2o=13. ! default value
             call getin_p("rad_h2o",rad_h2o)
             write(*,*)" rad_h2o = ",rad_h2o
 
-            write(*,*)"radius of H2O ice particles:"
+            write(*,*)"radius of H2O ice particles:"
             rad_h2o_ice=35. ! default value
             call getin_p("rad_h2o_ice",rad_h2o_ice)
             write(*,*)" rad_h2o_ice = ",rad_h2o_ice
 
-         else
+         else
 
             write(*,*)"Number mixing ratio of H2O water particles:"
@@ -122 +119 @@
             call getin_p("Nmix_h2o_ice",Nmix_h2o_ice)
             write(*,*)" Nmix_h2o_ice = ",Nmix_h2o_ice
-         endif
+         endif
 
       print*,'exit su_aer_radii'
@@ -173 +170 @@
                zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
                zfice = MIN(MAX(zfice,0.0),1.0)
-               zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
-               zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
+               zrad_liq  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o*pi*rhowater) )
+               zrad_ice  = CBRT( 3*pq(ig,l)/(4*Nmix_h2o_ice*pi*rhowaterice) )
                nueffrad(ig,l) = coef_chaud * (1.-zfice) + coef_froid * zfice
                zrad = zrad_liq * (1.-zfice) + zrad_ice * zfice
 
-               reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
+               reffrad(ig,l) = min(max(zrad,1.e-6),1000.e-6)
                enddo
             enddo      
@@ -213 +210 @@
 
       if (radfixed) then
-         reffliq(1:ngrid,1:nlayer)= rad_h2o
-         reffice(1:ngrid,1:nlayer)= rad_h2o_ice
+         reffliq(1:ngrid,1:nlayer)= rad_h2o
+         reffice(1:ngrid,1:nlayer)= rad_h2o_ice
       else
-         do k=1,nlayer
-           do i=1,ngrid
-             reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
-             reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
-           
-             reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
-             reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
-           enddo
-         enddo
+         do k=1,nlayer
+           do i=1,ngrid
+             reffliq(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o*pi*rhowater))
+             reffliq(i,k) = min(max(reffliq(i,k),1.e-6),1000.e-6)
+           
+             reffice(i,k) = CBRT(3*pql(i,k)/(4*Nmix_h2o_ice*pi*rhowaterice))
+             reffice(i,k) = min(max(reffice(i,k),1.e-6),1000.e-6)
+           enddo
+         enddo
       endif
 

trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90

@@ -1 @@
-      module radinc_h
+module radinc_h
 
-      implicit none
+  implicit none
 
-#include "dimensions.h"
-#include "bands.h"
-#include "scatterers.h"
+  include "bands.h"
+  include "scatterers.h"
 
 !======================================================================
@@ -57 +56 @@
 !----------------------------------------------------------------------
 
-      integer, parameter :: L_NLAYRAD  = llm
-      integer, parameter :: L_LEVELS   = 2*(llm-1)+3
-      integer, parameter :: L_NLEVRAD  = llm+1
+      integer,save :: L_NLAYRAD  ! = nbp_lev ! set by ini_radinc_h
+      integer,save :: L_LEVELS   ! = 2*(nbp_lev-1)+3 ! set by ini_radinc_h
+      integer,save :: L_NLEVRAD  ! = nbp_lev+1 ! set by ini_radinc_h
+!$OMP THREADPRIVATE(L_NLAYRAD,L_LEVELS,L_NLEVRAD)
 
       ! These are set in sugas_corrk
@@ -96 +96 @@
 !$OMP THREADPRIVATE(banddir)
 
-      end module radinc_h
+contains
+
+  subroutine ini_radinc_h(nbp_lev)
+  ! Initialize module variables
+  implicit none
+  integer,intent(in) :: nbp_lev
+  
+  L_NLAYRAD = nbp_lev
+  L_LEVELS = 2*(nbp_lev-1)+3
+  L_NLEVRAD = nbp_lev+1
+  
+  end subroutine
+
+end module radinc_h

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

@@ -20 +20 @@
 !        heat capacity are commons in comsoil.h
 !-----------------------------------------------------------------------
-
-!      include "dimensions.h"
-!      include "dimphys.h"
-
 
 c-----------------------------------------------------------------------

trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90

@@ -2 +2 @@
 
       ! inputs
-      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,iim,jjm,naerkind
+      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
       use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
       use datafile_mod, only: datadir, aerdir

trunk/LMDZ.GENERIC/libf/phystd/surf_heat_transp_mod.F90

@@ -4 +4 @@
 MODULE surf_heat_transp_mod
 
-
+IMPLICIT NONE  
+
+  ! Variables copied over from dyn3d dynamics:
+  REAL,SAVE,ALLOCATABLE :: fext(:)
+  REAL,SAVE,ALLOCATABLE :: unsairez(:)
+  REAL,SAVE,ALLOCATABLE :: unsaire(:)
+  REAL,SAVE,ALLOCATABLE :: cu(:)
+  REAL,SAVE,ALLOCATABLE :: cv(:)
+  REAL,SAVE,ALLOCATABLE :: cuvsurcv(:)
+  REAL,SAVE,ALLOCATABLE :: cvusurcu(:)
+  REAL,SAVE,ALLOCATABLE :: aire(:)
+  REAL,SAVE :: apoln
+  REAL,SAVE :: apols
+  REAL,SAVE,ALLOCATABLE :: aireu(:)
+  REAL,SAVE,ALLOCATABLE :: airev(:) 
+
+  LOGICAL,SAVE :: alpha_var
+  LOGICAL,SAVE :: slab_upstream
+  REAL,SAVE,ALLOCATABLE :: zmasqu(:)
+  REAL,SAVE,ALLOCATABLE :: zmasqv(:)
+  REAL,SAVE,ALLOCATABLE :: unsfv(:)
+  REAL,SAVE,ALLOCATABLE :: unsev(:)
+  REAL,SAVE,ALLOCATABLE :: unsfu(:)
+  REAL,SAVE,ALLOCATABLE :: unseu(:)
+
+  ! Routines usable only by routines within this module:
+  PRIVATE :: gr_fi_dyn, gr_dyn_fi
+  ! Routines from dyn3d, valid on global dynamics grid only:
+  PRIVATE :: grad,diverg,gr_u_scal,gr_v_scal
+  
 CONTAINS
-
-      SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
-
-      USE comhdiff_mod, ONLY: zmasqu,zmasqv
-
+  
+  SUBROUTINE ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez_,fext_,unsaire_,&
+                                  cu_,cuvsurcv_,cv_,cvusurcu_, &
+                                  aire_,apoln_,apols_, &
+                                  aireu_,airev_)
+    USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
+    IMPLICIT NONE
+    ! Transfer some variables from dyn3d dynamics
+    INTEGER,INTENT(IN) :: ip1jm
+    INTEGER,INTENT(IN) :: ip1jmp1
+    REAL,INTENT(IN) :: unsairez_(ip1jm)
+    REAL,INTENT(IN) :: fext_(ip1jm)
+    REAL,INTENT(IN) :: unsaire_(ip1jmp1)
+    REAL,INTENT(IN) :: cu_(ip1jmp1)
+    REAL,INTENT(IN) :: cuvsurcv_(ip1jm)
+    REAL,INTENT(IN) :: cv_(ip1jm)
+    REAL,INTENT(IN) :: cvusurcu_(ip1jmp1)
+    REAL,INTENT(IN) :: aire_(ip1jmp1)
+    REAL,INTENT(IN) :: apoln_
+    REAL,INTENT(IN) :: apols_
+    REAL,INTENT(IN) :: aireu_(ip1jmp1)
+    REAL,INTENT(IN) :: airev_(ip1jm)
+    
+    ! Sanity check
+    if ((ip1jm.ne.((nbp_lon+1)*(nbp_lat-1))).or. &
+        (ip1jmp1.ne.((nbp_lon+1)*nbp_lat))) then
+      write(*,*) "ini_surf_heat_transp Error: wrong array sizes"
+      stop
+    endif
+    
+    allocate(unsairez(ip1jm))
+    unsairez(:)=unsairez_(:)
+    allocate(fext(ip1jm))
+    fext(:)=fext_(:)
+    allocate(unsaire(ip1jmp1))
+    unsaire(:)=unsaire_(:)
+    allocate(cu(ip1jmp1))
+    cu(:)=cu_(:)
+    allocate(cuvsurcv(ip1jm))
+    cuvsurcv(:)=cuvsurcv_(:)
+    allocate(cv(ip1jm))
+    cv(:)=cv_(:)
+    allocate(cvusurcu(ip1jmp1))
+    cvusurcu(:)=cvusurcu_(:)
+    allocate(aire(ip1jmp1))
+    aire(:)=aire_(:)
+    apoln=apoln_
+    apols=apols_
+    allocate(aireu(ip1jmp1))
+    aireu(:)=aireu_(:)
+    allocate(airev(ip1jm))
+    airev(:)=airev_(:)
+    
+  END SUBROUTINE ini_surf_heat_transp
+  
+  SUBROUTINE ini_surf_heat_transp_mod
+    USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
+    IMPLICIT NONE
+    INTEGER :: ip1jm, ip1jmp1
+    
+    ip1jm=(nbp_lon+1)*(nbp_lat-1)
+    ip1jmp1=(nbp_lon+1)*nbp_lat
+  
+    allocate(zmasqu(ip1jmp1))
+    allocate(zmasqv(ip1jm))
+    allocate(unsfv(ip1jm))
+    allocate(unsev(ip1jm))
+    allocate(unsfu(ip1jmp1))
+    allocate(unseu(ip1jmp1))
+
+  END SUBROUTINE ini_surf_heat_transp_mod
+
+  SUBROUTINE divgrad_phy(ngrid,nlevs,temp,delta)
+
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
       IMPLICIT NONE
 
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-#include "comgeom.h"
-      
       INTEGER,INTENT(IN) :: ngrid, nlevs
       REAL,INTENT(IN) :: temp(ngrid,nlevs)
       REAL,INTENT(OUT) :: delta(ngrid,nlevs)
-      REAL delta_2d(ip1jmp1,nlevs)
+      REAL delta_2d((nbp_lon+1)*nbp_lat,nlevs)
       INTEGER :: ll
-      REAL ghx(ip1jmp1,nlevs), ghy(ip1jm,nlevs)
-
+      REAL ghx((nbp_lon+1)*nbp_lat,nlevs)
+      REAL ghy((nbp_lon+1)*(nbp_lat-1),nlevs)
+      INTEGER :: iip1,jjp1
+      
+      iip1=nbp_lon+1
+      jjp1=nbp_lat
 
       CALL gr_fi_dyn(nlevs,ngrid,iip1,jjp1,temp,delta_2d)
@@ -42 +140 @@
 
 
-      SUBROUTINE init_masquv(ngrid,zmasq)
-
-      USE comhdiff_mod, ONLY: zmasqu,zmasqv,unsfu,unsfv,unseu,unsev
-     
+  SUBROUTINE init_masquv(ngrid,zmasq)
+
+      USE mod_grid_phy_lmdz, only: nbp_lon, nbp_lat
       IMPLICIT NONE
-
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-#include "comgeom.h"
 
 
       INTEGER,INTENT(IN) :: ngrid
       REAL zmasq(ngrid)
-      REAL zmasq_2d(ip1jmp1)
-      REAL ff(ip1jm)
+      REAL zmasq_2d((nbp_lon+1)*nbp_lat)
+      REAL ff((nbp_lon+1)*(nbp_lat-1))
       REAL eps
       INTEGER i
-
+      INTEGER :: iim,iip1,jjp1,ip1jm,ip1jmp1
+      
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      jjp1=nbp_lat
+      ip1jm=(nbp_lon+1)*(nbp_lat-1)
+      ip1jmp1=(nbp_lon+1)*nbp_lat
 
 ! Masques u,v
@@ -104 +202 @@
 
 
-      SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
-  
-          use slab_ice_h 
-          USE comhdiff_mod, ONLY: zmasqu,zmasqv,unsfu,unsfv,unseu,unsev
+  SUBROUTINE slab_ekman2(ngrid,tx_phy,ty_phy,ts_phy,dt_phy)
+  
+      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+      USE slab_ice_h, ONLY: noceanmx
 
       IMPLICIT NONE
       
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-#include "comgeom.h"
-
       INTEGER,INTENT(IN) :: ngrid
       INTEGER ij
-      REAL txv(ip1jm),fluxm(ip1jm),tyv(ip1jm)
-      REAL fluxtm(ip1jm,noceanmx),fluxtz(ip1jmp1,noceanmx)
-      REAL tyu(ip1jmp1),txu(ip1jmp1),fluxz(ip1jmp1),fluxv(ip1jmp1)
-      REAL dt(ip1jmp1,noceanmx),ts(ip1jmp1,noceanmx)
+      REAL txv((nbp_lon+1)*(nbp_lat-1)),fluxm((nbp_lon+1)*(nbp_lat-1))
+      REAL tyv((nbp_lon+1)*(nbp_lat-1))
+      REAL fluxtm((nbp_lon+1)*(nbp_lat-1),noceanmx)
+      REAL fluxtz((nbp_lon+1)*nbp_lat,noceanmx)
+      REAL tyu((nbp_lon+1)*nbp_lat),txu((nbp_lon+1)*nbp_lat)
+      REAL fluxz((nbp_lon+1)*nbp_lat),fluxv((nbp_lon+1)*nbp_lat)
+      REAL dt((nbp_lon+1)*nbp_lat,noceanmx),ts((nbp_lon+1)*nbp_lat,noceanmx)
       REAL tx_phy(ngrid),ty_phy(ngrid)
       REAL dt_phy(ngrid,noceanmx),ts_phy(ngrid,noceanmx)
 
-
-
+      INTEGER iim,iip1,iip2,jjp1,ip1jm,ip1jmi1,ip1jmp1
+
+      iim=nbp_lon
+      iip1=nbp_lon+1
+      iip2=nbp_lon+2
+      jjp1=nbp_lon
+      ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+      ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+      ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
 
 ! Passage taux,y sur grilles 2d
@@ -226 +329 @@
   END SUBROUTINE slab_ekman2
 
- 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+  SUBROUTINE gr_fi_dyn(nfield,ngrid,im,jm,pfi,pdyn)
+  ! Transfer a variable on global "physics" grid to global "dynamics" grid
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
+  REAL,INTENT(IN) :: pfi(ngrid,nfield)
+  REAL,INTENT(OUT) :: pdyn(im,jm,nfield)
+  
+  INTEGER :: i,j,ifield,ig
+  
+  DO ifield=1,nfield
+    ! Handle poles
+    DO i=1,im
+      pdyn(i,1,ifield)=pfi(1,ifield)
+      pdyn(i,jm,ifield)=pfi(ngrid,ifield)
+    ENDDO
+    ! Other points
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pfi(ig,ifield),1,pdyn(1,j,ifield),1)
+      pdyn(im,j,ifield)=pdyn(1,j,ifield)
+    ENDDO
+  ENDDO ! of DO ifield=1,nfield
+  
+  END SUBROUTINE gr_fi_dyn
+  
+
+
+  SUBROUTINE gr_dyn_fi(nfield,im,jm,ngrid,pdyn,pfi)
+  ! Transfer a variable on global "dynamics" grid to global "physics" grid
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: im,jm,ngrid,nfield
+  REAL,INTENT(IN) :: pdyn(im,jm,nfield)
+  REAL,INTENT(OUT) :: pfi(ngrid,nfield)
+  
+  INTEGER j,ifield,ig
+
+  ! Sanity check:
+  IF(ngrid.NE.2+(jm-2)*(im-1)) THEN
+    WRITE(*,*) "gr_dyn_fi error, wrong sizes"
+    STOP
+  ENDIF
+  
+  ! Handle poles
+  CALL SCOPY(nfield,pdyn,im*jm,pfi,ngrid)
+  CALL SCOPY(nfield,pdyn(1,jm,1),im*jm,pfi(ngrid,1),ngrid)
+  ! Other points
+  DO ifield=1,nfield
+    DO j=2,jm-1
+      ig=2+(j-2)*(im-1)
+      CALL SCOPY(im-1,pdyn(1,j,ifield),1,pfi(ig,ifield),1)
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_dyn_fi
+
+
+
+  SUBROUTINE  grad(klevel,pg,pgx,pgy)
+  ! compute the covariant components x,y of the gradient of pg
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: pg((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgx((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(OUT) :: pgy((nbp_lon+1)*(nbp_lat-1),klevel)
+
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,ip1jm,ip1jmp1
+  
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  
+  DO l=1,klevel
+    DO ij=1,ip1jmp1-1
+      pgx(ij,l)=pg(ij+1,l)-pg(ij,l)
+    ENDDO
+    ! correction for pgx(ip1,j,l) ...
+    ! ... pgx(iip1,j,l)=pgx(1,j,l) ...
+    DO ij=iip1,ip1jmp1,iip1
+      pgx(ij,l)=pgx(ij-iim,l)
+    ENDDO
+    DO ij=1,ip1jm
+      pgy(ij,l)=pg(ij,l)-pg(ij+iip1,l)
+    ENDDO
+  ENDDO
+  
+  END SUBROUTINE grad
+
+
+
+  SUBROUTINE diverg(klevel,x,y,div)
+  ! compute the divergence of a vector field of components
+  ! x,y. y and y being covriant components
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: klevel
+  REAL,INTENT(IN) :: x((nbp_lon+1)*nbp_lat,klevel)
+  REAL,INTENT(IN) :: y((nbp_lon+1)*(nbp_lat-1),klevel)
+  REAL,INTENT(OUT) :: div((nbp_lon+1)*nbp_lat,klevel)
+  
+  INTEGER :: l,ij
+  INTEGER :: iim,iip1,iip2,ip1jm,ip1jmp1,ip1jmi1
+  
+  REAL :: aiy1(nbp_lon+1),aiy2(nbp_lon+1)
+  REAL :: sumypn,sumyps
+  REAL,EXTERNAL :: SSUM
+  
+  iim=nbp_lon
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  ip1jmi1=(nbp_lon+1)*(nbp_lat-1)-(nbp_lon+1) ! = ip1jm - iip1
+  
+  DO l=1,klevel
+    DO ij=iip2,ip1jm-1
+      div(ij+1,l)= &
+        cvusurcu(ij+1)*x(ij+1,l)-cvusurcu(ij)*x(ij,l)+ &
+        cuvsurcv(ij-iim)*y(ij-iim,l)-cuvsurcv(ij+1)*y(ij+1,l) 
+    ENDDO
+    ! correction for div(1,j,l) ...
+    ! ... div(1,j,l)= div(iip1,j,l) ...
+    DO ij=iip2,ip1jm,iip1
+      div(ij,l)=div(ij+iim,l)
+    ENDDO
+    ! at the poles
+    DO ij=1,iim
+      aiy1(ij)=cuvsurcv(ij)*y(ij,l)
+      aiy2(ij)=cuvsurcv(ij+ip1jmi1)*y(ij+ip1jmi1,l)
+    ENDDO
+    sumypn=SSUM(iim,aiy1,1)/apoln
+    sumyps=SSUM(iim,aiy2,1)/apols
+    DO ij=1,iip1
+      div(ij,l)=-sumypn
+      div(ij+ip1jm,l)=sumyps
+    ENDDO
+  ENDDO ! of DO l=1,klevel
+  
+  !!! CALL filtreg( div, jjp1, klevel, 2, 2, .TRUE., 1 )
+  
+  DO l=1,klevel
+    DO ij=iip2,ip1jm
+      div(ij,l)=div(ij,l)*unsaire(ij) 
+    ENDDO
+  ENDDO
+  
+  END SUBROUTINE diverg
+
+
+
+  SUBROUTINE gr_u_scal(nx,x_u,x_scal)
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: nx
+  REAL,INTENT(IN) :: x_u((nbp_lon+1)*nbp_lat,nx)
+  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+  
+  INTEGER :: l,ij
+  INTEGER :: iip1,jjp1,ip1jmp1
+  
+  iip1=nbp_lon+1
+  jjp1=nbp_lat
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  
+  DO l=1,nx
+    DO ij=ip1jmp1,2,-1
+      x_scal(ij,l)= &
+                   (aireu(ij)*x_u(ij,l)+aireu(ij-1)*x_u(ij-1,l)) &
+                  /(aireu(ij)+aireu(ij-1))
+    ENDDO
+  ENDDO
+  
+  CALL SCOPY(nx*jjp1,x_scal(iip1,1),iip1,x_scal(1,1),iip1)
+
+  END SUBROUTINE gr_u_scal
+
+
+
+  SUBROUTINE gr_v_scal(nx,x_v,x_scal)
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  IMPLICIT NONE
+  
+  INTEGER,INTENT(IN) :: nx
+  REAL,INTENT(IN) :: x_v((nbp_lon+1)*(nbp_lat-1),nx)
+  REAL,INTENT(OUT) :: x_scal((nbp_lon+1)*nbp_lat,nx)
+  
+  INTEGER :: l,ij
+  INTEGER :: iip1,iip2,ip1jm,ip1jmp1
+  
+  iip1=nbp_lon+1
+  iip2=nbp_lon+2
+  ip1jm=(nbp_lon+1)*(nbp_lat-1) ! = iip1*jjm
+  ip1jmp1=(nbp_lon+1)*nbp_lat ! = iip1*jjp1
+  
+  DO l=1,nx
+    DO ij=iip2,ip1jm
+      x_scal(ij,l)= &
+                   (airev(ij-iip1)*x_v(ij-iip1,l)+airev(ij)*x_v(ij,l)) &
+                  /(airev(ij-iip1)+airev(ij))
+    ENDDO
+    DO ij=1,iip1
+      x_scal(ij,l)=0.
+    ENDDO
+    DO ij=ip1jm+1,ip1jmp1
+      x_scal(ij,l)=0.
+    ENDDO
+  ENDDO
+
+  END SUBROUTINE gr_v_scal
+  
 END MODULE surf_heat_transp_mod
 

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

@@ -1 +1 @@
       SUBROUTINE vdif_kc(ngrid,nlay,dt,g,zlev,zlay,u,v,teta,cd,q2,km,kn)
       IMPLICIT NONE
-c.......................................................................
-!#include "dimensions.h"
-!#include "dimphys.h"
 c.......................................................................
 c

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

@@ -15 +15 @@
 c   --------------------------------------------------------------------
       IMPLICIT NONE
-c
-!#include "dimensions.h"
-!#include "dimphys.h"
-
 c
 c

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

@@ -40 +40 @@
 !=================================================================
       use surfdat_h, only: phisfi
+      use comgeomphy, only: airephy
       use time_phylmdz_mod, only: ecritphy, day_step, iphysiq, day_ini
       USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
      &                               is_master, gather
-      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat, nbp_lev
       implicit none
 
 ! Commons
-      include "dimensions.h"
-      include "paramet.h"
-      include "comgeom.h"
       include "netcdf.inc"
 
@@ -56 +55 @@
       character (len=*),intent(in) :: nom,titre,unite
       integer,intent(in) :: dim
-      real,intent(in) :: px(ngrid,llm)
+      real,intent(in) :: px(ngrid,nbp_lev)
 
 ! Local variables:
 
-      real*4 dx3(iip1,jjp1,llm) ! to store a 3D data set
-      real*4 dx2(iip1,jjp1)     ! to store a 2D (surface) data set
-      real*4 dx1(llm)           ! to store a 1D (column) data set
+      real*4 dx3(nbp_lon+1,nbp_lat,nbp_lev) ! to store a 3D data set
+      real*4 dx2(nbp_lon+1,nbp_lat)     ! to store a 2D (surface) data set
+      real*4 dx1(nbp_lev)           ! to store a 1D (column) data set
       real*4 dx0
 
@@ -68 +67 @@
 !$OMP THREADPRIVATE(date)
 
-      REAL phis(ip1jmp1)
+      REAL phis((nbp_lon+1),nbp_lat)
+      REAL area((nbp_lon+1),nbp_lat)
 
       integer irythme
       integer ierr,ierr2
-      integer iq
-      integer i,j,l,zmax , ig0
+      integer i,j,l, ig0
 
       integer,save :: zitau=0
@@ -102 +101 @@
 #ifdef CPP_PARA
 ! Added to work in parallel mode
-      real dx3_glop(klon_glo,llm)
-      real dx3_glo(iim,jjp1,llm) ! to store a global 3D data set
+      real dx3_glop(klon_glo,nbp_lev)
+      real dx3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set
       real dx2_glop(klon_glo)
-      real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+      real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
       real px2(ngrid)
-!      real dx1_glo(llm)          ! to store a 1D (column) data set
+!      real dx1_glo(nbp_lev)          ! to store a 1D (column) data set
 !      real dx0_glo
       real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
 #else
       real phisfi_glo(ngrid) ! surface geopotential on global physics grid
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -117 +118 @@
 !Sortie des variables au rythme voulu
 
-      irythme = ecritphy ! sortie au rythme de ecritphy
+      irythme = int(ecritphy) ! sortie au rythme de ecritphy
 !     irythme = iconser  ! sortie au rythme des variables de controle
 !     irythme = iphysiq  ! sortie a tous les pas physique
@@ -188 +189 @@
           ! Gather phisfi() geopotential on physics grid
           call Gather(phisfi,phisfi_glo)
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
 #else
          phisfi_glo(:)=phisfi(:)
+         areafi_glo(:)=airephy(:)
 #endif
 
@@ -216 +220 @@
          ierr = NF_ENDDEF(nid)
 
+         ! Build phis() and area()
+         do i=1,nbp_lon+1 ! poles
+           phis(i,1)=phisfi_glo(1)
+           phis(i,nbp_lat)=phisfi_glo(klon_glo)
+           ! for area, divide at the poles by nbp_lon
+           area(i,1)=areafi_glo(1)/nbp_lon
+           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+         enddo
+         do j=2,nbp_lat-1
+           ig0= 1+(j-2)*nbp_lon
+           do i=1,nbp_lon
+              phis(i,j)=phisfi_glo(ig0+i)
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           phis(nbp_lon+1,j)=phis(1,j)
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+         
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
-         call iniwrite(nid,day_ini,phis)
+         call iniwrite(nid,day_ini,phis,area)
 
          endif ! of if (is_master)
@@ -234 +256 @@
 
       if (ngrid.eq.1) then
-        ! in testphys1d, for the 1d version of the GCM, iphysiq 
-        ! should be most likely 1 (because no dyn!)
+        ! in testphys1d, for the 1d version of the GCM, iphysiq and irythme
+        ! are undefined; so set them to 1
         iphysiq=1 
+        irythme=1
         ! NB: 
       endif
@@ -296 +319 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -304 +327 @@
 !         Passage variable physique -->  variable dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-           DO l=1,llm
-             DO i=1,iip1
+           DO l=1,nbp_lev
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -344 +367 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
-           edges(3)=llm
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
+           edges(3)=nbp_lev
            edges(4)=1
 !#ifdef NC_DOUBLE
@@ -380 +403 @@
             call Grid1Dto2D_glo(dx2_glop,dx2_glo)
             ! copy dx2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=dx2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -390 +413 @@
 !         recast (copy) variable from physics grid to dynamics grid
 
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx2(i,1)=px(1,1)
-                dx2(i,jjp1)=px(ngrid,1)
+                dx2(i,nbp_lat)=px(ngrid,1)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx2(i,j)=px(ig0+i,1)
                 ENDDO
-                dx2(iip1,j)=dx2(1,j)
+                dx2(nbp_lon+1,j)=dx2(1,j)
              ENDDO
 #endif
@@ -426 +449 @@
            corner(2)=1
            corner(3)=ntime
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=1
 
@@ -457 +480 @@
 !         Passage variable physique -->  physique dynamique
 !         recast (copy) variable from physics grid to dynamics grid
-          do l=1,llm
+          do l=1,nbp_lev
             dx1(l)=px(1,l)
           enddo
@@ -479 +502 @@
            corner(2)=ntime
            
-           edges(1)=llm
+           edges(1)=nbp_lev
            edges(2)=1
 !#ifdef NC_DOUBLE
@@ -543 +566 @@
 
 #endif
+! of #ifndef MESOSCALE
       end

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

@@ -12 +12 @@
 ! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
 
-use comsoil_h, only: nsoilmx
+use comsoil_h, only: nsoilmx, inertiedat
+use comgeomphy, only: airephy
 use time_phylmdz_mod, only: ecritphy, day_step, iphysiq
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
+                              nbp_lon, nbp_lat
 
 implicit none
 
-#include"dimensions.h"
-#include"paramet.h"
-#include"netcdf.inc"
+include"netcdf.inc"
 
 ! Arguments:
@@ -31 +31 @@
 integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
 real,dimension(ngrid,nsoilmx),intent(in) :: px ! variable
-! Note: nsoilmx is a parameter set in 'comsoil_h'
 
 ! Local variables:
-real*4,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
-! Note iip1,jjp1 known from paramet.h; nsoilmx known from comsoil_h
-real*4,dimension(iip1,jjp1) :: data2 ! to store 2D data
+real*4,dimension(nbp_lon+1,nbp_lat,nsoilmx) :: data3 ! to store 3D data
+real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
 real*4 :: data0 ! to store 0D data
 integer :: i,j,l ! for loops
@@ -42 +40 @@
 
 real*4,save :: date ! time counter (in elapsed days)
+
+real :: inertia((nbp_lon+1),nbp_lat,nsoilmx)
+real :: area((nbp_lon+1),nbp_lat)
+
+real :: inertiafi_glo(klon_glo,nsoilmx)
+real :: areafi_glo(klon_glo)
+
 integer,save :: isample ! sample rate at which data is to be written to output
 integer,save :: ntime=0 ! counter to internally store time steps
@@ -60 +65 @@
 ! Added to work in parallel mode
 real dx3_glop(klon_glo,nsoilmx)
-real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
+real dx3_glo(nbp_lon,nbp_lat,nsoilmx) ! to store a global 3D data set
 real dx2_glop(klon_glo)
-real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
+real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
 real px2(ngrid)
 #endif
@@ -81 +86 @@
   
   ! Set output sample rate
-  isample=ecritphy ! same as for diagfi outputs
+  isample=int(ecritphy) ! same as for diagfi outputs
   ! Note ecritphy is known from control.h
   
@@ -91 +96 @@
     stop
    endif
+
+#ifdef CPP_PARA
+   ! Gather inertiedat() soil thermal inertia on physics grid
+   call Gather(inertiedat,inertiafi_glo)
+   ! Gather airephy() mesh area on physics grid
+   call Gather(airephy,areafi_glo)
+#else
+         inertiafi_glo(:,:)=inertiedat(:,:)
+         areafi_glo(:)=airephy(:)
+#endif
+
+   ! build inertia() and area()
+   do i=1,nbp_lon+1 ! poles
+     inertia(i,1,1:nsoilmx)=inertiafi_glo(1,1:nsoilmx)
+     inertia(i,nbp_lat,1:nsoilmx)=inertiafi_glo(klon_glo,1:nsoilmx)
+     ! for area, divide at the poles by nbp_lon
+     area(i,1)=areafi_glo(1)/nbp_lon
+     area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+   enddo
+   do j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     do i=1,nbp_lon
+        inertia(i,j,1:nsoilmx)=inertiafi_glo(ig0+i,1:nsoilmx)
+        area(i,j)=areafi_glo(ig0+i)
+     enddo
+     ! handle redundant point in longitude
+     inertia(nbp_lon+1,j,1:nsoilmx)=inertia(1,j,1:nsoilmx)
+     area(nbp_lon+1,j)=area(1,j)
+   enddo
+   
+   ! write "header" of file (longitudes, latitudes, geopotential, ...)
+   call iniwritesoil(nid,ngrid,inertia,area)
+
   endif ! of if (is_master)
-
-  ! Define dimensions and axis attributes
-  call iniwritesoil(nid,ngrid)
   
   ! set zitau to -1 to be compatible with zitau incrementation step below
@@ -148 +183 @@
     call Grid1Dto2D_glo(dx3_glop,dx3_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-    data3(iip1,:,:)=data3(1,:,:)
+    data3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+    data3(nbp_lon+1,:,:)=data3(1,:,:)
   endif
 !$OMP END MASTER
@@ -156 +191 @@
   do l=1,nsoilmx
     ! handle the poles
-    do i=1,iip1
+    do i=1,nbp_lon+1
       data3(i,1,l)=px(1,l)
-      data3(i,jjp1,l)=px(ngrid,l)
+      data3(i,nbp_lat,l)=px(ngrid,l)
     enddo
     ! rest of the grid
-    do j=2,jjm
-      ig0=1+(j-2)*iim
-      do i=1,iim
+    do j=2,nbp_lat-1
+      ig0=1+(j-2)*nbp_lon
+      do i=1,nbp_lon
         data3(i,j,l)=px(ig0+i,l)
       enddo
-      data3(iip1,j,l)=data3(1,j,l) ! extra (modulo) longitude
+      data3(nbp_lon+1,j,l)=data3(1,j,l) ! extra (modulo) longitude
     enddo
   enddo
@@ -195 +230 @@
   corners(4)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=nsoilmx
   edges(4)=1
@@ -223 +258 @@
     call Grid1Dto2D_glo(dx2_glop,dx2_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    data2(1:iim,:)=dx2_glo(1:iim,:)
-    data2(iip1,:)=data2(1,:)
+    data2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
+    data2(nbp_lon+1,:)=data2(1,:)
   endif
 !$OMP END MASTER
@@ -230 +265 @@
 #else
   ! handle the poles
-  do i=1,iip1
+  do i=1,nbp_lon+1
     data2(i,1)=px(1,1)
-    data2(i,jjp1)=px(ngrid,1)
+    data2(i,nbp_lat)=px(ngrid,1)
   enddo
   ! rest of the grid
-  do j=2,jjm
-    ig0=1+(j-2)*iim
-    do i=1,iim
+  do j=2,nbp_lat-1
+    ig0=1+(j-2)*nbp_lon
+    do i=1,nbp_lon
       data2(i,j)=px(ig0+i,1)
     enddo
-    data2(iip1,j)=data2(1,j) ! extra (modulo) longitude
+    data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
   enddo
 #endif
@@ -266 +301 @@
   corners(3)=ntime
   
-  edges(1)=iip1
-  edges(2)=jjp1
+  edges(1)=nbp_lon+1
+  edges(2)=nbp_lat
   edges(3)=1
   

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

@@ -43 +43 @@
 ! Addition by RW (2010) to allow OLR to be saved in .nc format
       use radinc_h, only : L_NSPECTI
-!      USE surfdat_h, only : phisfi
-#ifdef CPP_PARA
+      use comgeomphy, only: airephy
       use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-#endif
+      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat
       use time_phylmdz_mod, only: ecritphy, iphysiq, day_step, day_ini
       use callkeys_mod, only: iradia
@@ -53 +52 @@
       implicit none
 
-! Commons
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-!#include "control.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 ! Arguments on input:
@@ -73 +66 @@
 !      real dx0
 
-      real date
-
-!      REAL phis(ip1jmp1)
-
       integer irythme
       integer ierr
@@ -82 +71 @@
       integer i,j,l,zmax , ig0
 
-      integer zitau
-      character firstnom*20
-      SAVE firstnom
-      SAVE zitau
-      SAVE date
-      data firstnom /'1234567890'/
-      data zitau /0/
+      integer,save :: zitau=0
+      character(len=20),save :: firstnom='1234567890'
+      real,save :: date
 !$OMP THREADPRIVATE(firstnom,zitau,date)
 
@@ -100 +85 @@
       integer, dimension(4) :: edges,corner
 
+      real area((nbp_lon+1),nbp_lat)
 ! added by RDW for OLR output
-       real dx3(iip1,jjp1,L_NSPECTI) ! to store the data set
+       real dx3(nbp_lon+1,nbp_lat,L_NSPECTI) ! to store the data set
 
 #ifdef CPP_PARA
 ! Added to work in parallel mode
       real dx3_glop(klon_glo,L_NSPECTI)
-      real dx3_glo(iim,jjp1,L_NSPECTI) ! to store a global 3D data set
-#else
-      logical,parameter :: is_master=.true.
-      logical,parameter :: is_mpi_root=.true.
+      real dx3_glo(nbp_lon,nbp_lat,L_NSPECTI) ! to store a global 3D data set
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
+#else
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -139 +125 @@
          endif
 
+#ifdef CPP_PARA
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
+#else
+         areafi_glo(:)=airephy(:)
+#endif
          ! Create the NetCDF file
          if (is_master) then
@@ -159 +151 @@
          ierr = NF_ENDDEF(nid)
 
-!         call gr_fi_dyn(1,size(phisfi_glo),iip1,jjp1,phisfi_glo,phis)
+         ! Build area()
+         do i=1,nbp_lon+1 ! poles
+           ! divide at the poles by nbp_lon
+           area(i,1)=areafi_glo(1)/nbp_lon
+           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+         enddo
+         do j=2,nbp_lat-1
+           ig0= 1+(j-2)*nbp_lon
+           do i=1,nbp_lon
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+
          ! write "header" of file (longitudes, latitudes, area, ...)
-         call iniwrite_specIR(nid,day_ini)
+         call iniwrite_specIR(nid,day_ini,area)
          endif ! of if (is_master)
 
@@ -215 +221 @@
              endif
 
-             write(6,*)'WRITEDIAGSPEC: date= ', date
+             write(6,*)'WRITEDIAGSPECIR: date= ', date
            endif ! of if (is_master)
         end if ! of if (nom.eq.firstnom)
@@ -233 +239 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -240 +246 @@
 #else
            DO l=1,L_NSPECTI
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -279 +285 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=L_NSPECTI
            edges(4)=1

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

@@ -43 +43 @@
 ! Addition by RW (2010) to allow OSR to be saved in .nc format
       use radinc_h, only : L_NSPECTV
-!      USE surfdat_h
-#ifdef CPP_PARA
+      use comgeomphy, only: airephy
       use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-#endif
+      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo,
+     &                              nbp_lon, nbp_lat
       use time_phylmdz_mod, only: ecritphy, iphysiq, day_step, day_ini
       use callkeys_mod, only: iradia
@@ -53 +52 @@
       implicit none
 
-! Commons
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "paramet.h"
-!#include "control.h"
-#include "comgeom.h"
-#include "netcdf.inc"
+      include "netcdf.inc"
 
 ! Arguments on input:
@@ -73 +66 @@
 !      real dx0
 
-      real date
-
-!      REAL phis(ip1jmp1)
-
       integer irythme
       integer ierr
@@ -82 +71 @@
       integer i,j,l,zmax , ig0
 
-      integer zitau
-      character firstnom*20
-      SAVE firstnom
-      SAVE zitau
-      SAVE date
-      data firstnom /'1234567890'/
-      data zitau /0/
+      integer,save :: zitau=0
+      character(len=20),save :: firstnom='1234567890'
+      real,save :: date
 !$OMP THREADPRIVATE(firstnom,zitau,date)
 
@@ -100 +85 @@
       integer, dimension(4) :: edges,corner
 
+      real area((nbp_lon+1),nbp_lat)
 ! added by RDW for OSR output
-       real dx3(iip1,jjp1,L_NSPECTV) ! to store the data set
+       real dx3(nbp_lon+1,nbp_lat,L_NSPECTV) ! to store the data set
 
 #ifdef CPP_PARA
 ! Added to work in parallel mode
       real dx3_glop(klon_glo,L_NSPECTV)
-      real dx3_glo(iim,jjp1,L_NSPECTV) ! to store a global 3D data set
-#else
-      logical,parameter :: is_master=.true.
-      logical,parameter :: is_mpi_root=.true.
+      real dx3_glo(nbp_lon,nbp_lat,L_NSPECTV) ! to store a global 3D data set
+      real areafi_glo(klon_glo) ! mesh area on global physics grid
+#else
+      real areafi_glo(ngrid) ! mesh area on global physics grid
 #endif
 
@@ -138 +124 @@
          endif
 
+#ifdef CPP_PARA
+          ! Gather airephy() mesh area on physics grid
+          call Gather(airephy,areafi_glo)
+#else
+         areafi_glo(:)=airephy(:)
+#endif
          ! Create the NetCDF file
          if (is_master) then
@@ -158 +150 @@
          ierr = NF_ENDDEF(nid)
 
+         ! Build area()
+         do i=1,nbp_lon+1 ! poles
+           ! divide at the poles by nbp_lon
+           area(i,1)=areafi_glo(1)/nbp_lon
+           area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
+         enddo
+         do j=2,nbp_lat-1
+           ig0= 1+(j-2)*nbp_lon
+           do i=1,nbp_lon
+              area(i,j)=areafi_glo(ig0+i)
+           enddo
+           ! handle redundant point in longitude
+           area(nbp_lon+1,j)=area(1,j)
+         enddo
+
          ! write "header" of file (longitudes, latitudes, geopotential, ...)
-!         call gr_fi_dyn(1,ngrid,iip1,jjp1,phisfi,phis)
-!         call iniwrite(nid,day_ini,phis)
-         call iniwrite_specVI(nid,day_ini)
+         call iniwrite_specVI(nid,day_ini,area)
          endif ! of if (is_master)
 
@@ -215 +220 @@
              endif
 
-             write(6,*)'WRITEDIAGSPEC: date= ', date
+             write(6,*)'WRITEDIAGSPECVI: date= ', date
            endif ! of if (is_master)
         end if ! of if (nom.eq.firstnom)
@@ -233 +238 @@
             call Grid1Dto2D_glo(dx3_glop,dx3_glo)
             ! copy dx3_glo() to dx3(:) and add redundant longitude
-            dx3(1:iim,:,:)=dx3_glo(1:iim,:,:)
-            dx3(iip1,:,:)=dx3(1,:,:)
+            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
+            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
           endif
 !$OMP END MASTER
@@ -240 +245 @@
 #else
            DO l=1,L_NSPECTV
-             DO i=1,iip1
+             DO i=1,nbp_lon+1
                 dx3(i,1,l)=px(1,l)
-                dx3(i,jjp1,l)=px(ngrid,l)
+                dx3(i,nbp_lat,l)=px(ngrid,l)
              ENDDO
-             DO j=2,jjm
-                ig0= 1+(j-2)*iim
-                DO i=1,iim
+             DO j=2,nbp_lat-1
+                ig0= 1+(j-2)*nbp_lon
+                DO i=1,nbp_lon
                    dx3(i,j,l)=px(ig0+i,l)
                 ENDDO
-                dx3(iip1,j,l)=dx3(1,j,l)
+                dx3(nbp_lon+1,j,l)=dx3(1,j,l)
              ENDDO
            ENDDO
@@ -279 +284 @@
            corner(4)=ntime
 
-           edges(1)=iip1
-           edges(2)=jjp1
+           edges(1)=nbp_lon+1
+           edges(2)=nbp_lat
            edges(3)=L_NSPECTV
            edges(4)=1

trunk/LMDZ.GENERIC/libf/phystd/wstats.F90

@@ -1 +1 @@
 subroutine wstats(ngrid,nom,titre,unite,dim,px)
 
+use statto_mod, only: istats,istime,count
 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
-use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-use statto_mod, only: istats,istime,count
-
+use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo, &
+                              nbp_lon, nbp_lat, nbp_lev
 implicit none
 
-#include "dimensions.h"
-!#include "dimphys.h"
-#include "netcdf.inc"
+include "netcdf.inc"
 
 integer,intent(in) :: ngrid
 character (len=*),intent(in) :: nom,titre,unite
 integer,intent(in) :: dim
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
-real,intent(in) :: px(ngrid,llm)
-real, dimension(iip1,jjp1,llm) :: mean3d,sd3d,dx3
-real, dimension(iip1,jjp1) :: mean2d,sd2d,dx2
+real,intent(in) :: px(ngrid,nbp_lev)
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: mean3d,sd3d,dx3
+real, dimension(nbp_lon+1,nbp_lat) :: mean2d,sd2d,dx2
 character (len=50) :: namebis
 character (len=50), save :: firstvar
@@ -34 +30 @@
 ! Added to work in parallel mode
 #ifdef CPP_PARA
-real px3_glop(klon_glo,llm) ! to store a 3D data set on global physics grid
-real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px3_glop(klon_glo,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
 real px2_glop(klon_glo) ! to store a 2D data set on global physics grid
-real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
 real px2(ngrid)
-real px3(ngrid,llm)
+real px3(ngrid,nbp_lev)
 #else
 ! When not running in parallel mode:
-real px3_glop(ngrid,llm) ! to store a 3D data set on global physics grid
-real px3_glo(iim,jjp1,llm) ! to store a global 3D data set on global lonxlat grid
+real px3_glop(ngrid,nbp_lev) ! to store a 3D data set on global physics grid
+real px3_glo(nbp_lon,nbp_lat,nbp_lev) ! to store a global 3D data set on global lonxlat grid
 real px2_glop(ngrid) ! to store a 2D data set on global physics grid
-real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
+real px2_glo(nbp_lon,nbp_lat) ! to store a 2D data set on global lonxlat grid
 #endif
 
@@ -67 +63 @@
 #ifdef CPP_PARA
  if (dim.eq.3) then
-  px3(1:ngrid,1:llm)=px(1:ngrid,1:llm)
+  px3(1:ngrid,1:nbp_lev)=px(1:ngrid,1:nbp_lev)
   ! Gather fieds on a "global" (without redundant longitude) array
   call Gather(px3,px3_glop)
@@ -74 +70 @@
     call Grid1Dto2D_glo(px3_glop,px3_glo)
     ! copy dx3_glo() to dx3(:) and add redundant longitude
-    dx3(1:iim,:,:)=px3_glo(1:iim,:,:)
-    dx3(iip1,:,:)=dx3(1,:,:)
+    dx3(1:nbp_lon,:,:)=px3_glo(1:nbp_lon,:,:)
+    dx3(nbp_lon+1,:,:)=dx3(1,:,:)
   endif
 !$OMP END MASTER
@@ -87 +83 @@
             call Grid1Dto2D_glo(px2_glop,px2_glo)
             ! copy px2_glo() to dx2(:) and add redundant longitude
-            dx2(1:iim,:)=px2_glo(1:iim,:)
-            dx2(iip1,:)=dx2(1,:)
+            dx2(1:nbp_lon,:)=px2_glo(1:nbp_lon,:)
+            dx2(nbp_lon+1,:)=dx2(1,:)
           endif
 !$OMP END MASTER
@@ -95 +91 @@
 #else
   if (dim.eq.3) then
-    px3_glop(:,1:llm)=px(:,1:llm)
+    px3_glop(:,1:nbp_lev)=px(:,1:nbp_lev)
 !  Passage variable physique -->  variable dynamique
-    DO l=1,llm
-      DO i=1,iim
+    DO l=1,nbp_lev
+      DO i=1,nbp_lon
          px3_glo(i,1,l)=px(1,l)
-         px3_glo(i,jjp1,l)=px(ngrid,l)
+         px3_glo(i,nbp_lat,l)=px(ngrid,l)
       ENDDO
-      DO j=2,jjm
-         ig0= 1+(j-2)*iim
-         DO i=1,iim
+      DO j=2,nbp_lat-1
+         ig0= 1+(j-2)*nbp_lon
+         DO i=1,nbp_lon
             px3_glo(i,j,l)=px(ig0+i,l)
          ENDDO
@@ -112 +108 @@
     px2_glop(:)=px(:,1)
 !    Passage variable physique -->  physique dynamique
-   DO i=1,iim
+   DO i=1,nbp_lon
      px2_glo(i,1)=px(1,1)
-     px2_glo(i,jjp1)=px(ngrid,1)
+     px2_glo(i,nbp_lat)=px(ngrid,1)
    ENDDO
-   DO j=2,jjm
-     ig0= 1+(j-2)*iim
-     DO i=1,iim
+   DO j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     DO i=1,nbp_lon
         px2_glo(i,j)=px(ig0+i,1)
      ENDDO
@@ -198 +194 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/iip1,jjp1,llm,1/)
+      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
       mean3d(:,:,:)=0
       sd3d(:,:,:)=0
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lev,1,0/)
       mean2d(:,:)=0
       sd2d(:,:)=0
@@ -211 +207 @@
    if (dim.eq.3) then
       start=(/1,1,1,indx/)
-      sizes=(/iip1,jjp1,llm,1/)
+      sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean3d)
@@ -226 +222 @@
    else if (dim.eq.2) then
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lat,1,0/)
 #ifdef NC_DOUBLE
       ierr = NF_GET_VARA_DOUBLE(nid,meanid,start,sizes,mean2d)
@@ -245 +241 @@
 
 if (dim.eq.3) then
-  dx3(1:iim,:,:)=px3_glo(:,:,:)
-  dx3(iip1,:,:)=dx3(1,:,:)
+  dx3(1:nbp_lon,:,:)=px3_glo(:,:,:)
+  dx3(nbp_lon+1,:,:)=dx3(1,:,:)
 else ! dim.eq.2
-  dx2(1:iim,:)=px2_glo(:,:)
-  dx2(iip1,:)=dx2(1,:)
+  dx2(1:nbp_lon,:)=px2_glo(:,:)
+  dx2(nbp_lon+1,:)=dx2(1,:)
 endif
 
@@ -290 +286 @@
 !======================================================
 subroutine inivar(nid,varid,ngrid,dim,indx,px,ierr)
+use mod_grid_phy_lmdz, only : nbp_lon, nbp_lat, nbp_lev
 
 implicit none
 
-include "dimensions.h"
-!include "dimphys.h"
 include "netcdf.inc"
 
 integer, intent(in) :: nid,varid,dim,indx,ngrid
-real, dimension(ngrid,llm), intent(in) :: px
+real, dimension(ngrid,nbp_lev), intent(in) :: px
 integer, intent(out) :: ierr
-
-integer,parameter :: iip1=iim+1
-integer,parameter :: jjp1=jjm+1
 
 integer :: l,i,j,ig0
 integer, dimension(4) :: start,sizes
-real, dimension(iip1,jjp1,llm) :: dx3
-real, dimension(iip1,jjp1) :: dx2
+real, dimension(nbp_lon+1,nbp_lat,nbp_lev) :: dx3
+real, dimension(nbp_lon+1,nbp_lat) :: dx2
 
 if (dim.eq.3) then
 
    start=(/1,1,1,indx/)
-   sizes=(/iip1,jjp1,llm,1/)
+   sizes=(/nbp_lon+1,nbp_lat,nbp_lev,1/)
 
 !  Passage variable physique -->  variable dynamique
 
-   DO l=1,llm
-      DO i=1,iip1
+   DO l=1,nbp_lev
+      DO i=1,nbp_lon+1
          dx3(i,1,l)=px(1,l)
-         dx3(i,jjp1,l)=px(ngrid,l)
+         dx3(i,nbp_lat,l)=px(ngrid,l)
       ENDDO
-      DO j=2,jjm
-         ig0= 1+(j-2)*iim
-         DO i=1,iim
+      DO j=2,nbp_lat-1
+         ig0= 1+(j-2)*nbp_lon
+         DO i=1,nbp_lon
             dx3(i,j,l)=px(ig0+i,l)
          ENDDO
-         dx3(iip1,j,l)=dx3(1,j,l)
+         dx3(nbp_lon+1,j,l)=dx3(1,j,l)
       ENDDO
    ENDDO
@@ -339 +331 @@
 
       start=(/1,1,indx,0/)
-      sizes=(/iip1,jjp1,1,0/)
+      sizes=(/nbp_lon+1,nbp_lat,1,0/)
 
 !    Passage variable physique -->  physique dynamique
 
-  DO i=1,iip1
+  DO i=1,nbp_lon+1
      dx2(i,1)=px(1,1)
-     dx2(i,jjp1)=px(ngrid,1)
+     dx2(i,nbp_lat)=px(ngrid,1)
   ENDDO
-  DO j=2,jjm
-     ig0= 1+(j-2)*iim
-     DO i=1,iim
+  DO j=2,nbp_lat-1
+     ig0= 1+(j-2)*nbp_lon
+     DO i=1,nbp_lon
         dx2(i,j)=px(ig0+i,1)
      ENDDO
-     dx2(iip1,j)=dx2(1,j)
+     dx2(nbp_lon+1,j)=dx2(1,j)
   ENDDO
 
@@ -380 +372 @@
 implicit none
 
-#include "netcdf.inc"
+include "netcdf.inc"
 
 integer,intent(in) :: nid ! NetCDF file ID