Changeset 1216

Timestamp:

Apr 3, 2014, 9:09:47 AM (11 years ago)

Author:

emillour

Message:

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

• Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
• In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
• Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
• Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
• added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
• added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

EM

Location:

trunk/LMDZ.GENERIC

Files:

• README (modified) (1 diff)
• libf/dyn3d/calfis.F (modified) (10 diffs)
• libf/dyn3d/control_mod.F90 (moved) (moved from trunk/LMDZ.GENERIC/libf/dyn3d/control.h) (1 diff)
• libf/dyn3d/defrun_new.F (modified) (1 diff)
• libf/dyn3d/dynetat0.F (modified) (5 diffs)
• libf/dyn3d/dynredem.F (modified) (7 diffs)
• libf/dyn3d/gcm.F (modified) (19 diffs)
• libf/dyn3d/infotrac.F90 (moved) (moved from trunk/LMDZ.GENERIC/libf/dyn3d/iniadvtrac.F) (8 diffs)
• libf/dyn3d/iniconst.F (modified) (2 diffs)
• libf/dyn3d/inidissip.F (modified) (2 diffs)
• libf/dyn3d/logic.h (modified) (1 diff)
• libf/dyn3d/mod_const_mpi.F90 (added)
• libf/dyn3d/test_period.F (modified) (3 diffs)
• libf/grid/dimension/makdim (modified) (3 diffs)
• libf/phystd/callkeys.h (modified) (1 diff)
• libf/phystd/comgeomphy.F90 (added)
• libf/phystd/comsoil_h.F90 (modified) (1 diff)
• libf/phystd/condense_cloud.F90 (modified) (2 diffs)
• libf/phystd/dimphy.F90 (added)
• libf/phystd/dimphys.h (modified) (1 diff)
• libf/phystd/ini_archive.F (modified) (2 diffs)
• libf/phystd/inifis.F (modified) (5 diffs)
• libf/phystd/iniphysiq.F90 (added)
• libf/phystd/inistats.F (modified) (9 diffs)
• libf/phystd/init_phys_lmdz.F90 (added)
• libf/phystd/iniwrite.F (modified) (4 diffs)
• libf/phystd/iniwritesoil.F90 (modified) (8 diffs)
• libf/phystd/iostart.F90 (added)
• libf/phystd/kcm1d.F90 (modified) (5 diffs)
• libf/phystd/lect_start_archive.F (modified) (16 diffs)
• libf/phystd/mkstat.F90 (modified) (3 diffs)
• libf/phystd/mod_grid_phy_lmdz.F90 (added)
• libf/phystd/mod_phys_lmdz_mpi_data.F90 (added)
• libf/phystd/mod_phys_lmdz_mpi_transfert.F90 (added)
• libf/phystd/mod_phys_lmdz_omp_data.F90 (added)
• libf/phystd/mod_phys_lmdz_omp_transfert.F90 (added)
• libf/phystd/mod_phys_lmdz_para.F90 (added)
• libf/phystd/mod_phys_lmdz_transfert_para.F90 (added)
• libf/phystd/newstart.F (modified) (23 diffs)
• libf/phystd/phyetat0.F90 (moved) (moved from trunk/LMDZ.GENERIC/libf/phystd/phyetat0.F) (4 diffs)
• libf/phystd/phyredem.F90 (moved) (moved from trunk/LMDZ.GENERIC/libf/phystd/physdem1.F) (1 diff)
• libf/phystd/physiq.F90 (modified) (12 diffs)
• libf/phystd/planetwide_mod.F90 (added)
• libf/phystd/rcm1d.F (modified) (21 diffs)
• libf/phystd/soil.F (modified) (6 diffs)
• libf/phystd/soil_settings.F (modified) (15 diffs)
• libf/phystd/start2archive.F (modified) (11 diffs)
• libf/phystd/stelang.F (modified) (2 diffs)
• libf/phystd/surfini.F (modified) (4 diffs)
• libf/phystd/tabfi.F (modified) (6 diffs)
• libf/phystd/write_archive.F (modified) (1 diff)
• libf/phystd/writediagfi.F (modified) (10 diffs)
• libf/phystd/writediagsoil.F90 (modified) (8 diffs)
• libf/phystd/writediagspecIR.F (modified) (3 diffs)
• libf/phystd/writediagspecVI.F (modified) (3 diffs)
• libf/phystd/wstats.F90 (modified) (4 diffs)
• makegcm_g95 (modified) (7 diffs)
• makegcm_gfortran (modified) (9 diffs)
• makegcm_ifort (modified) (7 diffs)
• makegcm_pgf90 (modified) (7 diffs)

trunk/LMDZ.GENERIC/README

@@ -1001 +1001 @@
 == 05/12/2013 == JL
 - corrected sugascorrk to work in the two band gray aproximation with -b 1x1 and NGAUSS=2
+
+== 03/04/2014 == EM
+Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
+in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
+in the near future)
+- Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce
+  full recomputation of the model) option
+- In dyn3d: converted control.h to module control_mod.F90 and converted
+  iniadvtrac.F to module infotrac.F90
+- Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
+- Rearanged input/outputs routines everywhere to handle serial/MPI cases.
+  physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write
+  routines for startfi files are gathered in module iostart.F90
+- added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90,
+  dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90,
+  mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90,
+  mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90,
+  mod_phys_lmdz_transfert_para.F90 in phymars
+  and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
+- added created generic routines 'planetwide_maxval' and 'planetwide_minval',
+  in module "planetwide_mod", that enable obtaining the max and min of a field
+  over the whole planet. This should be further imroved with computation of
+  means (possibly area weighed), etc.

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

@@ -6 +6 @@
 c    Auteur :  P. Le Van, F. Hourdin 
 c   .........
-
+      USE infotrac, ONLY: tname, nqtot
       IMPLICIT NONE
 c=======================================================================
@@ -71 +71 @@
 #include "comvert.h"
 #include "comgeom2.h"
-#include "control.h"
-
-#include "advtrac.h"
+!#include "control.h"
+
+!#include "advtrac.h"
 !! this is to get tnom (tracers name)
 
@@ -85 +85 @@
       REAL pteta(iip1,jjp1,llm)
       REAL pmasse(iip1,jjp1,llm)
-      REAL pq(iip1,jjp1,llm,nqmx)
+      REAL pq(iip1,jjp1,llm,nqtot)
       REAL pphis(iip1,jjp1)
       REAL pphi(iip1,jjp1,llm)
@@ -92 +92 @@
       REAL pducov(iip1,jjp1,llm)
       REAL pdteta(iip1,jjp1,llm)
-      REAL pdq(iip1,jjp1,llm,nqmx)
+      REAL pdq(iip1,jjp1,llm,nqtot)
 c
       REAL pw(iip1,jjp1,llm)
@@ -103 +103 @@
       REAL pdufi(iip1,jjp1,llm)
       REAL pdhfi(iip1,jjp1,llm)
-      REAL pdqfi(iip1,jjp1,llm,nqmx)
+      REAL pdqfi(iip1,jjp1,llm,nqtot)
       REAL pdpsfi(iip1,jjp1)
       logical tracer
@@ -116 +116 @@
 c
       REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
-      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
+      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
 c
       REAL zvervel(ngridmx,llm)
 c
       REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
-      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
+      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
       REAL zdpsrf(ngridmx)
 c
@@ -170 +170 @@
 
 c
-      IF (firstcal) THEN
-         latfi(1)=rlatu(1)
-         lonfi(1)=0.
-         DO j=2,jjm
-            DO i=1,iim
-               latfi((j-2)*iim+1+i)= rlatu(j)
-               lonfi((j-2)*iim+1+i)= rlonv(i)
-            ENDDO
-         ENDDO
-         latfi(ngridmx)= rlatu(jjp1)
-         lonfi(ngridmx)= 0.
-
-         ! build airefi(), mesh area on physics grid
-         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
-         ! Poles are single points on physics grid
-         airefi(1)=airefi(1)*iim
-         airefi(ngridmx)=airefi(ngridmx)*iim
-
-         CALL inifis(ngridmx,llm,day_ini,daysec,dtphys,
-     .                latfi,lonfi,airefi,rad,g,r,cpp)
-      ENDIF
+!      IF (firstcal) THEN
+!         latfi(1)=rlatu(1)
+!         lonfi(1)=0.
+!         DO j=2,jjm
+!            DO i=1,iim
+!               latfi((j-2)*iim+1+i)= rlatu(j)
+!               lonfi((j-2)*iim+1+i)= rlonv(i)
+!            ENDDO
+!         ENDDO
+!         latfi(ngridmx)= rlatu(jjp1)
+!         lonfi(ngridmx)= 0.
+!
+!         ! build airefi(), mesh area on physics grid
+!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
+!         ! Poles are single points on physics grid
+!         airefi(1)=airefi(1)*iim
+!         airefi(ngridmx)=airefi(ngridmx)*iim
+!
+!         CALL inifis(ngridmx,llm,day_ini,daysec,dtphys,
+!     .                latfi,lonfi,airefi,rad,g,r,cpp)
+!      ENDIF
 
 c
@@ -278 +278 @@
 c   43.bis Taceurs (en kg/kg)
 c   --------------------------
-      DO iq=1,nqmx
+      DO iq=1,nqtot
          DO l=1,llm
             zqfi(1,l,iq) = pq(1,1,l,iq)
@@ -425 +425 @@
 
       CALL physiq (ngridmx,llm,nq,
-     .     tnom,
+     .     tname,
      ,     debut,lafin,
      ,     rday_ecri,heure,dtphys,
@@ -467 +467 @@
 
 
-c   62. humidite specifique
+c   62. traceurs
 c   ---------------------
 
-      DO iq=1,nqmx
+      DO iq=1,nqtot
          DO l=1,llm
             DO i=1,iip1

trunk/LMDZ.GENERIC/libf/dyn3d/control_mod.F90

@@ -1 @@
-!-----------------------------------------------------------------------
-! INCLUDE 'control.h'
-! For Fortran 77/Fortran 90 compliance always use line continuation
-! symbols '&' in columns 73 and 6
-!
 
-      COMMON/control/nday,day_step,                                     &
-     &              iperiod,iconser,idissip,iphysiq ,                   &
-     &              periodav,ecritphy,anneeref
+module control_mod
 
-      INTEGER   nday,day_step,iperiod,iconser,                          &
-     &          idissip,iphysiq,anneeref
-      REAL periodav, ecritphy
+  implicit none
 
-!-----------------------------------------------------------------------
+  integer,save :: nday ! # of days to run
+  integer,save :: day_step ! # of dynamical time steps per day
+  integer,save :: iperiod  ! make a Matsuno step before avery iperiod-1 LF steps
+  integer,save :: iconser !
+  integer,save :: idissip ! apply dissipation every idissip dynamical step
+  integer,save :: iphysiq ! call physics every iphysiq dynamical steps
+  integer,save :: anneeref ! reference year # ! not used
+  real,save :: periodav
+  integer,save :: ecritphy ! output data in "diagfi.nc" every ecritphy dynamical steps 
+
+end module control_mod

trunk/LMDZ.GENERIC/libf/dyn3d/defrun_new.F

@@ -38 +38 @@
       USE ioipsl_getincom
       use sponge_mod,only: callsponge,nsponge,mode_sponge,tetasponge
+      use control_mod,only: nday, day_step, iperiod, anneeref,
+     &                      iconser, idissip, iphysiq, ecritphy
       IMPLICIT NONE
 
 #include "dimensions.h"
 #include "paramet.h"
-#include "control.h"
+!#include "control.h"
 #include "logic.h"
 #include "serre.h"

trunk/LMDZ.GENERIC/libf/dyn3d/dynetat0.F

@@ -1 +1 @@
       SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
      .                    teta,q,masse,ps,phis,time)
+      use infotrac, only: tname, nqtot
       IMPLICIT NONE
 
@@ -30 +31 @@
 #include "serre.h"
 #include "logic.h"
-#include"advtrac.h"
+!#include"advtrac.h"
 
 c   Arguments:
@@ -328 +329 @@
 !           WRITE(str3(2:3),'(i2.2)') iq
 !           ierr =  NF_INQ_VARID (nid, str3, nvarid)
-! NB: tracers are now read in using their name ('tnom' from advtrac.h)
-!           write(*,*) "  loading tracer:",trim(tnom(iq))
-           ierr=NF_INQ_VARID(nid,tnom(iq),nvarid)
+! NB: tracers are now read in using their name ('tname' from infotrac)
+!           write(*,*) "  loading tracer:",trim(tname(iq))
+           ierr=NF_INQ_VARID(nid,tname(iq),nvarid)
            IF (ierr .NE. NF_NOERR) THEN
 !              PRINT*, "dynetat0: Le champ <"//str3//"> est absent"
-              PRINT*, "dynetat0: Le champ <"//trim(tnom(iq))//
+              PRINT*, "dynetat0: Le champ <"//trim(tname(iq))//
      &                "> est absent"
               PRINT*, "          Il est donc initialise a zero"
@@ -346 +347 @@
              IF (ierr .NE. NF_NOERR) THEN
 !                 PRINT*, "dynetat0: Lecture echouee pour "//str3
-               PRINT*, "dynetat0: Lecture echouee pour "//trim(tnom(iq))
+               PRINT*,"dynetat0: Lecture echouee pour "//trim(tname(iq))
                CALL abort
              ENDIF
@@ -354 +355 @@
 c        case when new tracer are added in addition to old ones
              write(*,*)'tracers 1 to ', nqold,'were already present'
-             write(*,*)'tracers ', nqold+1,' to ', nqmx,'are new'
+             write(*,*)'tracers ', nqold+1,' to ', nqtot,'are new'
              write(*,*)' and initialized to zero'
-             q(:,:,:,nqold+1:nqmx)=0.0
+             q(:,:,:,nqold+1:nqtot)=0.0
 !             yes=' '
 !            do while ((yes.ne.'y').and.(yes.ne.'n'))

trunk/LMDZ.GENERIC/libf/dyn3d/dynredem.F

@@ -1 +1 @@
       SUBROUTINE dynredem0(fichnom,idayref,anneeref,phis,nq)
+      use infotrac, only: tname
       IMPLICIT NONE
 c=======================================================================
@@ -16 +17 @@
 #include "netcdf.inc"
 #include "serre.h"
-#include "advtrac.h"
+!#include "advtrac.h"
 c   Arguments:
 c   ----------
@@ -902 +903 @@
 !               ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,
 !     .                          "Traceurs "//str3)
-             txt="Traceur "//trim(tnom(iq))
-#ifdef NC_DOUBLE
-               ierr=NF_DEF_VAR(nid,tnom(iq),NF_DOUBLE,4,dims4,nvarid)
-#else
-               ierr=NF_DEF_VAR(nid,tnom(iq),NF_FLOAT,4,dims4,nvarid)
+             txt="Traceur "//trim(tname(iq))
+#ifdef NC_DOUBLE
+               ierr=NF_DEF_VAR(nid,tname(iq),NF_DOUBLE,4,dims4,nvarid)
+#else
+               ierr=NF_DEF_VAR(nid,tname(iq),NF_FLOAT,4,dims4,nvarid)
 #endif
                ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",
@@ -954 +955 @@
       SUBROUTINE dynredem1(fichnom,time,
      .                     vcov,ucov,teta,q,nq,masse,ps)
+      use infotrac, only: nqtot, tname
       IMPLICIT NONE
 c=================================================================
@@ -963 +965 @@
 #include "comvert.h"
 #include "comgeom.h"
-#include"advtrac.h"
+!#include"advtrac.h"
 
       INTEGER nq, l
@@ -969 +971 @@
       REAL teta(ip1jmp1,llm)                   
       REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
-      REAL q(iip1,jjp1,llm,nqmx)
+      REAL q(iip1,jjp1,llm,nqtot)
       REAL q3d(iip1,jjp1,llm) !temporary variable
       CHARACTER*(*) fichnom
@@ -1052 +1054 @@
 !            WRITE(str3(2:3),'(i2.2)') iq
 !            ierr = NF_INQ_VARID(nid, str3, nvarid)
-            ierr=NF_INQ_VARID(nid,tnom(iq),nvarid)
+            ierr=NF_INQ_VARID(nid,tname(iq),nvarid)
             IF (ierr .NE. NF_NOERR) THEN
 !               PRINT*, "Variable "//str3//" n est pas definie"
-              PRINT*, "Variable "//trim(tnom(iq))//" n est pas definie"
+              PRINT*,"Variable "//trim(tname(iq))//" n est pas definie"
               CALL abort
             ENDIF

trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F

@@ -1 +1 @@
       PROGRAM gcm
 
+      use infotrac, only: iniadvtrac, nqtot, iadv
       use sponge_mod,only: callsponge,mode_sponge,sponge
+      use control_mod, only: nday, day_step, iperiod, iphysiq,
+     &                       iconser, ecritphy, idissip
+      use comgeomphy, only: initcomgeomphy
       IMPLICIT NONE
 
@@ -42 +46 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
+!#include "control.h"
 #include "ener.h"
 #include "netcdf.inc"
 #include "serre.h"
 #include "tracstoke.h"
-#include"advtrac.h"
+!#include"advtrac.h"
 
       INTEGER*4  iday ! jour julien
@@ -56 +60 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       real, dimension(ip1jmp1,llm) :: teta   ! temperature potentielle 
-      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
+      REAL,allocatable :: q(:,:,:)           ! champs advectes
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL pext(ip1jmp1)                     ! pression  extensive
@@ -79 +83 @@
 c   tendances dynamiques
       REAL dv(ip1jm,llm),du(ip1jmp1,llm)
-      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
+      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
+      REAL,ALLOCATABLE :: dq(:,:,:)
 
 c   tendances de la dissipation
@@ -87 +92 @@
 c   tendances physiques
       REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
-      REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
+      REAL dhfi(ip1jmp1,llm),dpfi(ip1jmp1)
+      REAL,ALLOCATABLE :: dqfi(:,:,:)
 
 c   variables pour le fichier histoire
@@ -123 +129 @@
       LOGICAL tracer
           data tracer/.true./
-      INTEGER nq
+!      INTEGER nq
 
 C Calendrier
@@ -150 +156 @@
       REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
 
+c-----------------------------------------------------------------------
+c    variables pour l'initialisation de la physique :
+c    ------------------------------------------------
+      INTEGER ngridmx
+      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
+      REAL zcufi(ngridmx),zcvfi(ngridmx)
+      REAL latfi(ngridmx),lonfi(ngridmx)
+      REAL airefi(ngridmx)
+      SAVE latfi, lonfi, airefi
+      INTEGER i,j
 
 c-----------------------------------------------------------------------
@@ -159 +175 @@
 
 c-----------------------------------------------------------------------
-c  Initialize tracers using iniadvtrac (Ehouarn, oct 2008)
-      CALL iniadvtrac(nq,numvanle)
-
-
-      CALL dynetat0("start.nc",nqmx,vcov,ucov,
+      CALL defrun_new( 99, .TRUE. )
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! FH 2008/05/02
+! A nettoyer. On ne veut qu'une ou deux routines d'interface 
+! dynamique -> physique pour l'initialisation
+!#ifdef CPP_PHYS
+      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+      call initcomgeomphy
+!#endif
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+! Initialize tracers
+      CALL iniadvtrac(nqtot,numvanle)
+! Allocation de la tableau q : champs advectes   
+      allocate(q(ip1jmp1,llm,nqtot))
+      allocate(dq(ip1jmp1,llm,nqtot))
+      allocate(dqfi(ip1jmp1,llm,nqtot))
+
+      CALL dynetat0("start.nc",nqtot,vcov,ucov,
      .              teta,q,masse,ps,phis, time_0)
-
-      CALL defrun_new( 99, .TRUE. )
 
 c  on recalcule eventuellement le pas de temps
@@ -196 +225 @@
      *                tetagdiv, tetagrot , tetatemp              )
 c
+
+c-----------------------------------------------------------------------
+c   Initialisation de la physique :
+c   -------------------------------
+
+!      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
+         latfi(1)=rlatu(1)
+         lonfi(1)=0.
+         zcufi(1) = cu(1)
+         zcvfi(1) = cv(1)
+         DO j=2,jjm
+            DO i=1,iim
+               latfi((j-2)*iim+1+i)= rlatu(j)
+               lonfi((j-2)*iim+1+i)= rlonv(i)
+               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
+               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
+            ENDDO
+         ENDDO
+         latfi(ngridmx)= rlatu(jjp1)
+         lonfi(ngridmx)= 0.
+         zcufi(ngridmx) = cu(ip1jm+1)
+         zcvfi(ngridmx) = cv(ip1jm-iim)
+
+         ! build airefi(), mesh area on physics grid
+         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
+         ! Poles are single points on physics grid
+         airefi(1)=airefi(1)*iim
+         airefi(ngridmx)=airefi(ngridmx)*iim
+
+! Initialisation de la physique: pose probleme quand on tourne 
+! SANS physique, car iniphysiq.F est dans le repertoire phy[]... 
+! Il faut une cle CPP_PHYS
+!#ifdef CPP_PHYS
+!         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
+         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys,
+     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
+     &                1)
+!     &                iflag_phys)
+!#endif
+!         call_iniphys=.false.
+!      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
 
       CALL pression ( ip1jmp1, ap, bp, ps, p       )
@@ -229 +299 @@
      . 'c''est a dire du jour',i7,3x,'au jour',i7//)
 
-      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
+      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqtot)
 
       ecripar = .TRUE.
@@ -237 +307 @@
 
 c   Quelques initialisations pour les traceurs
-      call initial0(ijp1llm*nqmx,dq)
+      call initial0(ijp1llm*nqtot,dq)
 c     istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
 c     istphy=istdyn/iphysiq
@@ -328 +398 @@
        IF( forward. OR . leapf )  THEN
 
-        DO iq = 1, nqmx
+        DO iq = 1, nqtot
 c
          IF ( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
             CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
 
-         ELSE IF( iq.EQ. nqmx )   THEN
+         ELSE IF( iq.EQ. nqtot )   THEN
 c
             iapp_tracvl = 5
@@ -341 +411 @@
 c
 
-            CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
+            CALL vanleer(numvanle,iapp_tracvl,nqtot,q,pbaru,pbarv,
      *                      p, masse, dq,  iadv(1), teta, pk     )
 
@@ -413 +483 @@
 
 
-        CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
+        CALL calfis( nqtot, lafin ,rdayvrai,rday_ecri,time  ,
      $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
      $     du,dv,dteta,dq,w, dufi,dvfi,dhfi,dqfi,dpfi,tracer)
@@ -421 +491 @@
 c      ------------------------------
 !        if(1.eq.2)then
-          CALL addfi( nqmx, dtphys, leapf, forward   ,
+          CALL addfi( nqtot, dtphys, leapf, forward   ,
      $                  ucov, vcov, teta , q   ,ps , masse,
      $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
@@ -556 +626 @@
 c                 iav=0
 c              ENDIF
-c              CALL writedynav(histaveid, nqmx, itau,vcov ,
+c              CALL writedynav(histaveid, nqtot, itau,vcov ,
 c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
 c           ENDIF
@@ -569 +639 @@
        CALL test_period ( ucov,vcov,teta,q,p,phis )
        CALL dynredem1("restart.nc",0.0,
-     .                     vcov,ucov,teta,q,nqmx,masse,ps)
+     .                     vcov,ucov,teta,q,nqtot,masse,ps)
 
               CLOSE(99)
@@ -636 +706 @@
                   iav=0
                ENDIF
-c              CALL writedynav(histaveid, nqmx, itau,vcov ,
+c              CALL writedynav(histaveid, nqtot, itau,vcov ,
 c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
 
@@ -644 +714 @@
                  IF(itau.EQ.itaufin)
      . CALL dynredem1("restart.nc",0.0,
-     .                     vcov,ucov,teta,q,nqmx,masse,ps)
+     .                     vcov,ucov,teta,q,nqtot,masse,ps)
 
                  forward = .TRUE.

trunk/LMDZ.GENERIC/libf/dyn3d/infotrac.F90

@@ +1 @@
+MODULE infotrac
+
+IMPLICIT NONE
+! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
+  INTEGER, SAVE :: nqtot
+  INTEGER,allocatable :: iadv(:)   ! tracer advection scheme number
+  CHARACTER(len=20),allocatable ::  tname(:) ! tracer name
+
+CONTAINS
+
       subroutine iniadvtrac(nq,numvanle)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -8 +18 @@
       IMPLICIT NONE
 
-#include "dimensions.h"
-#include "advtrac.h"
-#include "control.h"
+!#include "dimensions.h"
+!#include "advtrac.h"
+!#include "control.h"
 
 ! routine arguments:
@@ -20 +30 @@
       INTEGER :: iq
       INTEGER :: ierr
-
-
-      if (nqmx > 0) then
+      CHARACTER(len=3) :: qname
 
 ! Look for file traceur.def
-      OPEN(90,file='traceur.def',form='formatted',status='old',
-     &        iostat=ierr)
+      OPEN(90,file='traceur.def',form='formatted',status='old', &
+              iostat=ierr)
       IF (ierr.eq.0) THEN
         write(*,*) "iniadvtrac: Reading file traceur.def"
@@ -35 +43 @@
           write(*,*) "   (first line of traceur.def) "
           stop
-        else
-          ! check that the number of tracers is indeed nqmx
-          if (nq.ne.nqmx) then
-            write(*,*) "iniadvtrac: error, wrong number of tracers:"
-            write(*,*) "nq=",nq," whereas nqmx=",nqmx
-            stop
-          endif
         endif
+        
+        ! allocate arrays:
+        allocate(iadv(nq))
+        allocate(tname(nq))
         
         ! initialize advection schemes to Van-Leer for all tracers
@@ -49 +54 @@
         enddo
         
-
-
-!     MODIFICATION TO TEST OTHER SCHEMES BY RDW
-!        do iq=1,nq
-!           iadv(iq)=1
-!        enddo
-!        print*,'IADV SET TO 1 IN iniadvtrac!!!!'
-
         do iq=1,nq
         ! minimal version, just read in the tracer names, 1 per line
-          read(90,*,iostat=ierr) tnom(iq)
+          read(90,*,iostat=ierr) tname(iq)
           if (ierr.ne.0) then
             write(*,*) 'iniadvtrac: error reading tracer names...'
@@ -65 +62 @@
           endif
         enddo !of do iq=1,nq
-      ELSE
-        write(*,*) "iniadvtrac: can't find file traceur.def..."
-        stop
+        close(90) ! done reading tracer names, close file
       ENDIF ! of IF (ierr.eq.0)
 
-c  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
-c  ...................................................................
-c
-c     iadv = 1    shema  transport type "humidite specifique LMD"  
-c     iadv = 2    shema   amont
-c     iadv = 3    shema  Van-leer
-c     iadv = 4    schema  Van-leer + humidite specifique
-c                        Modif F.Codron
-c
-c
-      DO  iq = 1, nqmx
-       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
-     * ,' pour le traceur no ', iq
-       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'
-
-     * ,' traceur no ', iq
-       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
-     * ,'le traceur no ', iq
+!  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
+!  ...................................................................
+!
+!     iadv = 1    shema  transport type "humidite specifique LMD"  
+!     iadv = 2    shema   amont
+!     iadv = 3    shema  Van-leer
+!     iadv = 4    schema  Van-leer + humidite specifique
+!                        Modif F.Codron
+! 
+!
+      DO  iq = 1, nq-1
+       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'&
+       ,' pour le traceur no ', iq
+       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'  &
+       ,' traceur no ', iq
+       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour' &
+       ,'le traceur no ', iq
 
        IF( iadv(iq).EQ.4 )  THEN
-         PRINT *,' Le shema  Van-Leer + humidite specifique ',
-     * ' est  uniquement pour la vapeur d eau .'
+         PRINT *,' Le shema  Van-Leer + humidite specifique ',          &
+       ' est  uniquement pour la vapeur d eau .'
          PRINT *,' Corriger iadv( ',iq, ')  et repasser ! '
          CALL ABORT
@@ -97 +91 @@
 
        IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 )   THEN
-        PRINT *,' Erreur dans le choix de iadv (nqmx).Corriger et '
-     * ,' repasser car  iadv(iq) = ', iadv(iq)
+        PRINT *,' Erreur dans le choix de iadv (nqtot).Corriger et '    &
+       ,' repasser car  iadv(iq) = ', iadv(iq)
          CALL ABORT
        ENDIF
       ENDDO
 
-!!!! AS: compiler complains about iadv(nqmx) when there is nqmx=0
-!!!! AS: so I commented those lines and changed nqmx-1 for nqmx above
-!       IF( iadv(nqmx).EQ.1 ) PRINT *,' Choix du shema humidite '
-!     * ,'specifique pour la vapeur d''eau'
-!       IF( iadv(nqmx).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'
-!     * ,' vapeur d''eau '
-!       IF( iadv(nqmx).EQ.3 ) PRINT *,' Choix du shema  Van-Leer '
-!     * ,' pour la vapeur d''eau'
-!       IF( iadv(nqmx).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + '
-!     * ,' humidite specifique pour la vapeur d''eau'
+       IF( iadv(nq).EQ.1 ) PRINT *,' Choix du shema humidite '          &
+       ,'specifique pour la vapeur d''eau'
+       IF( iadv(nq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'  &
+       ,' vapeur d''eau '
+       IF( iadv(nq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer '         &
+       ,' pour la vapeur d''eau'
+       IF( iadv(nq).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + '       &
+       ,' humidite specifique pour la vapeur d''eau'
 !
-!c
-!!       IF( (iadv(nqmx).LE.0).OR.(iadv(nqmx).GT.4) )   THEN
-!!     MODIFICATION TO TEST WITHOUT TRACER ADVECTION BY RDW
-!       IF( (iadv(nqmx).LT.0).OR.(iadv(nqmx).GT.4) )   THEN
-!        PRINT *,' Erreur dans le choix de iadv (nqmx).Corriger et '
-!     * ,' repasser car  iadv(nqmx) = ', iadv(nqmx)
-!         CALL ABORT
-!       ENDIF
+       IF( (iadv(nq).LE.0).OR.(iadv(nq).GT.4) )   THEN
+        PRINT *,' Erreur dans le choix de iadv (nqtot).Corriger et '    &
+       ,' repasser car  iadv(nqtot) = ', iadv(nqtot)
+         CALL ABORT
+       ENDIF
 
       first = .TRUE.
-      numvanle = nqmx + 1
-      DO  iq = 1, nqmx
+      numvanle = nq + 1
+      DO  iq = 1, nq
         IF(((iadv(iq).EQ.3).OR.(iadv(iq).EQ.4)).AND.first ) THEN
           numvanle = iq
@@ -131 +120 @@
         ENDIF
       ENDDO
-c
-      DO  iq = 1, nqmx
+!
+      DO  iq = 1, nq
 
       IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
-          PRINT *,' Il y a discontinuite dans le choix du shema de ',
-     *    'Van-leer pour les traceurs . Corriger et repasser . '
+          PRINT *,' Il y a discontinuite dans le choix du shema de ',   &
+          'Van-leer pour les traceurs . Corriger et repasser . '
            CALL ABORT
       ENDIF
 
       ENDDO
-c
-      endif ! of if nqmx > 0
+!
+      end subroutine iniadvtrac
 
-      end
+END MODULE infotrac

trunk/LMDZ.GENERIC/libf/dyn3d/iniconst.F

@@ -1 +1 @@
       SUBROUTINE iniconst
 
+      use control_mod, only: iphysiq, idissip
       IMPLICIT NONE
 c
@@ -13 +14 @@
 #include "comconst.h"
 #include "temps.h"
-#include "control.h"
+!#include "control.h"
 #include "comvert.h"
 

trunk/LMDZ.GENERIC/libf/dyn3d/inidissip.F

@@ -8 +8 @@
 c   -------------
 
+      use control_mod, only: idissip, iperiod
       IMPLICIT NONE
 #include "dimensions.h"
@@ -14 +15 @@
 #include "comconst.h"
 #include "comvert.h"
-#include "control.h"
+!#include "control.h"
 
       LOGICAL lstardis

trunk/LMDZ.GENERIC/libf/dyn3d/logic.h

@@ -3 +3 @@
 
       COMMON/logic/ purmats,physic,forward,leapf,apphys,grireg,
-     *  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,hybrid,autozlevs
+     *  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,
+     &  hybrid,autozlevs
 
       LOGICAL purmats,physic,forward,leapf,apphys,grireg,statcl,conser,

trunk/LMDZ.GENERIC/libf/dyn3d/test_period.F

@@ -1 +1 @@
       SUBROUTINE test_period ( ucov, vcov, teta, q, p, phis )
+      
+      USE infotrac, ONLY: nqtot
+      IMPLICIT NONE
 c
 c     Auteur : P. Le Van  
@@ -14 +17 @@
 c
       REAL ucov(ip1jmp1,llm), vcov(ip1jm,llm), teta(ip1jmp1,llm) ,
-     ,      q(ip1jmp1,llm,nqmx), p(ip1jmp1,llmp1), phis(ip1jmp1)
+     ,      q(ip1jmp1,llm,nqtot), p(ip1jmp1,llmp1), phis(ip1jmp1)
 c
 c   .....  Variables  locales  .....
@@ -51 +54 @@
       
 c
-      DO nq =1, nqmx
+      DO nq =1, nqtot
         DO l =1, llm
           DO ij = 1, ip1jmp1, iip1

trunk/LMDZ.GENERIC/libf/grid/dimension/makdim

@@ -1 +1 @@
 #!/bin/bash
 
-nqmx=$1
-shift
+#nqmx=$1
+#shift
 for i in $* ; do
    list=$list.$i
 done
-
-fichdim=dimensions${list}.t${nqmx}
+fichdim=dimensions${list}
+#fichdim=dimensions${list}.t${nqmx}
 
 echo $fichdim
@@ -18 +18 @@
    jm=$2
    lm=$3
-   n2=$1
+#   n2=$1
    ndm=1
 
@@ -51 +51 @@
 !   INCLUDE 'dimensions.h'
 !
-!   dimensions.h contient les dimensions du modele
-!   ndm est tel que iim=2**ndm
-!   nqmx est la dimension de la variable traceur q
+!   dimensions.h contains the model dimensions
+!   
 !-----------------------------------------------------------------------
 
-      INTEGER, parameter :: iim= $im
+      INTEGER, parameter :: iim=$im
       INTEGER, parameter :: jjm=$jm
       INTEGER, parameter :: llm=$lm
       INTEGER, parameter :: ndm=$ndm
-
-      integer, parameter :: nqmx=$nqmx
 
 !-----------------------------------------------------------------------

trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

@@ -3 +3 @@
 ! symbols '&' in columns 73 and 6
 !
-      COMMON/callkeys/callrad,corrk,calldifv,UseTurbDiff,calladj        &
+! Group commons according to their type for minimal performance impact
+
+      COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj      &
      &   , co2cond,callsoil                                             &
-     &   , season,diurnal,tlocked,rings_shadow,iradia,lwrite            &
-     &   , iaervar,iddist,topdustref,callstats,calleofdump              &
+     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
+     &   , callstats,calleofdump                                        &
      &   , enertest                                                     &
      &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
-     &   , Nmix_co2, radfixed, dusttau                                  &
+     &   , radfixed                                                     &
      &   , meanOLR, specOLR                                             &
-     &   , kastprof, Tstrat                                             &
-     &   , nosurf, intheat, flatten, oblate                             &     
-     &   , newtonian, tau_relax, testradtimes                           &
+     &   , kastprof                                                     &
+     &   , nosurf, oblate                                               &     
+     &   , newtonian, testradtimes                                      &
      &   , check_cpp_match, force_cpp                                   &
      &   , rayleigh                                                     &
-     &   , stelbbody, stelTbb                                           &
-     &   , tplanet                                                      &
-     &   , obs_tau_col_tropo                                            &
-     &   , obs_tau_col_strato                                           &
-     &   , pres_bottom_tropo                                            &
-     &   , pres_top_tropo                                               &
-     &   , pres_bottom_strato                                           &
-     &   , pres_top_strato                                              &
-     &   , size_tropo                                                   &
-     &   , size_strato                                                  &
-     &   , startype, Fat1AU                                             &
+     &   , stelbbody                                                    &
      &   , nearco2cond                                                  &
-     &   , tracer, mass_redistrib, varactive, varfixed, satval          &
+     &   , tracer, mass_redistrib, varactive, varfixed                  &
      &   , sedimentation,water,watercond,waterrain                      &
      &   , aeroco2,aeroh2o,aeroh2so4,aeroback2lay                       &
      &   , aerofixco2,aerofixh2o                                        &
-     &   , hydrology, sourceevol, icetstep, albedosnow                  &
-     &   , maxicethick, Tsaldiff                                        &
-     &   , CLFfixval, CLFvarying                                        &
-     &   , n2mixratio                                                   &
-     &   , co2supsat                                                    &
-     &   , cloudlvl                                                     &
-     &   , pceil                                                        &
-     &   , Rmean, J2, MassPlanet                                        &
+     &   , hydrology, sourceevol                                        &
+     &   , CLFvarying                                                   &
      &   , strictboundcorrk                                                                                     
 
+      COMMON/callkeys_i/iaervar,iddist,iradia,startype
+      
+      COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb,     &
+     &                  Tstrat,tplanet,obs_tau_col_tropo,               &
+     &                  obs_tau_col_strato,pres_bottom_tropo,           &
+     &                  pres_top_tropo,pres_bottom_strato,              &
+     &                  pres_top_strato,size_tropo,size_strato,satval,  &
+     &                  CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
+     &                  maxicethick,Tsaldiff,tau_relax,cloudlvl,        &
+     &                  icetstep,intheat,flatten,Rmean,J2,MassPlanet
+      
       logical callrad,corrk,calldifv,UseTurbDiff                        &
      &   , calladj,co2cond,callsoil                                     &

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

@@ +1 @@
+module comsoil_h
 
-       module comsoil_h
+implicit none
+! nsoilmx : number of subterranean layers
+!integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
+!integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
+  integer, parameter :: nsoilmx = 18 
 
-       implicit none
+  real,save,allocatable,dimension(:) :: layer      ! soil layer depths
+  real,save,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
+  real,save,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia
+  real,save :: volcapa    ! soil volumetric heat capacity
+       ! NB: volcapa is read fromn control(35) from physicq start file
+       !     in physdem (or set via tabfi, or initialized in
+       !                 soil_settings.F)
 
-       real,allocatable,dimension(:) :: layer      ! soil layer depths
-       real,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
-       real,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia
-       real volcapa    ! soil volumetric heat capacitya
+contains
 
-       end module comsoil_h
+  subroutine ini_comsoil_h(ngrid)
+  
+  implicit none
+  integer,intent(in) :: ngrid ! number of atmospheric columns
+  
+    if (.not.allocated(layer)) allocate(layer(nsoilmx)) !soil layer depths
+    if (.not.allocated(mlayer)) allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
+    if (.not.allocated(inertiedat)) allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia
+  
+  end subroutine ini_comsoil_h
 
+end module comsoil_h
+

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

@@ -7 +7 @@
 
       use radinc_h, only : naerkind
-      use gases_h
+      use gases_h, only: gfrac, igas_co2
       use radii_mod, only : co2_reffrad
       use aerosol_mod, only : iaero_co2
-      USE surfdat_h
-      USE comgeomfi_h
-      USE tracer_h
+      USE surfdat_h, only: albedodat, albedice, emisice, emissiv
+      USE comgeomfi_h, only: lati
+      USE tracer_h, only: noms, rho_co2
 
 
@@ -453 +453 @@
 !     --------------------------------------------------------------------
       DO ig=1,ngrid
-         IF(ig.GT.ngrid/2+1) THEN
-            icap=2
+         IF(lati(ig).LT.0.) THEN
+            icap=2 ! Southern Hemisphere
          ELSE
-            icap=1
+            icap=1 ! Nortnern hemisphere
          ENDIF
 

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

@@ -7 +7 @@
 ! nlayermx : number of atmospheric layers
       integer, parameter :: nlayermx = llm 
-! nsoilmx : number of subterranean layers
+! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
       !integer, parameter :: nsoilmx = 4 ! for a test 
-      integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
+      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
       !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
 !-----------------------------------------------------------------------

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

@@ -36 +36 @@
 
       USE comsoil_h
- 
+!      use control_mod
       implicit none
 
@@ -49 +49 @@
 #include "logic.h"
 #include "serre.h"
-#include "control.h"
+!#include "control.h"
 
 #include "netcdf.inc"

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

@@ -1 @@
-      SUBROUTINE inifis(ngrid,nlayer,
+      SUBROUTINE inifis(ngrid,nlayer,nq,
      $           day_ini,pdaysec,ptimestep,
      $           plat,plon,parea,
@@ -6 +6 @@
       use radinc_h, only : naerkind
       use datafile_mod, only: datadir
-      use comdiurn_h
-
-      !! to be conservative, but why this include here?
-      use surfdat_h
-      use comsaison_h
-
-      USE comgeomfi_h
+      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
+      use comgeomfi_h, only: long, lati, area, totarea
+      use comsoil_h, only: ini_comsoil_h
+      use control_mod, only: ecritphy
 
 !=======================================================================
@@ -61 +58 @@
 
 
-      REAL prad,pg,pr,pcpp,pdaysec,ptimestep
+      REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
  
-      INTEGER ngrid,nlayer
-      REAL plat(ngrid),plon(ngrid),parea(ngrid)
+      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
+      REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
       integer day_ini
       INTEGER ig,ierr
@@ -102 +99 @@
          STOP
       ENDIF
+
+      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
+      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
+      call getin("ecritphy",ecritphy)
 
 ! --------------------------------------------------------------
@@ -691 +692 @@
          coslon(ig)=cos(plon(ig))
       ENDDO
-         
+
       pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
 
-      RETURN
+      ! allocate "comsoil_h" arrays
+      call ini_comsoil_h(ngrid)
+      
       END

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

@@ -1 +1 @@
       subroutine inistats(ierr)
 
-#ifdef CPP_PARA
       use mod_phys_lmdz_para, only : is_master
-#endif
       implicit none
 
@@ -14 +12 @@
 #include "netcdf.inc"
 
-#ifndef CPP_PARA
-      logical,parameter :: is_master=.true.
-#endif
       integer,intent(out) :: ierr
       integer :: nid
@@ -32 +27 @@
       nsteppd=nint(daysec/dtphys)
       write (*,*) 'nsteppd=',nsteppd
-      write (*,*) 'istime=',istime
-
-
-
       if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec)
-     ,   stop'Dans Instat:  1jour .ne. n pas physiques'
-
+     &   stop'Dans Instat:  1jour .ne. n pas physiques'
 
       if(mod(nsteppd,istime).ne.0)
-     ,   stop'Dans Instat:  1jour .ne. n*istime pas physiques'
+     &   stop'Dans Instat:  1jour .ne. n*istime pas physiques'
 
       istats=nsteppd/istime
@@ -57 +47 @@
       ! only the master needs do this
 
-      ierr = NF_CREATE("stats.nc",NF_CLOBBER,nid)
+      ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
       if (ierr.ne.NF_NOERR) then
          write (*,*) NF_STRERROR(ierr)
@@ -70 +60 @@
 
       ierr = NF_ENDDEF(nid)
-      call def_var(nid,"Time","Time",
-     .            "days since 0000-00-0 00:00:00",1,
-     .            idim_time,nvarid,ierr)
+      call def_var_stats(nid,"Time","Time",
+     &            "days since 0000-00-0 00:00:00",1,
+     &            idim_time,nvarid,ierr)
 ! Time is initialised later by mkstats subroutine
 
-      call def_var(nid,"latitude","latitude","degrees_north",1,
-     .            idim_lat,nvarid,ierr)
+      call def_var_stats(nid,"latitude","latitude",
+     &            "degrees_north",1,idim_lat,nvarid,ierr)
 #ifdef NC_DOUBLE
       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
@@ -82 +72 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
 #endif
-      call def_var(nid,"longitude","East longitude","degrees_east",1,
-     .            idim_lon,nvarid,ierr)
+      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)
@@ -106 +96 @@
       ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
 #endif 
-      call def_var(nid,"aps","hybrid pressure at midlayers"," ",
-     .            1,idim_llm,nvarid,ierr)
+      call def_var_stats(nid,"aps","hybrid pressure at midlayers"
+     & ," ",1,idim_llm,nvarid,ierr)
 #ifdef NC_DOUBLE
       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
@@ -114 +104 @@
 #endif
 
-      call def_var(nid,"bps","hybrid sigma at midlayers"," ",
-     .            1,idim_llm,nvarid,ierr)
+      call def_var_stats(nid,"bps","hybrid sigma at midlayers"
+     & ," ",1,idim_llm,nvarid,ierr)
 #ifdef NC_DOUBLE
       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
@@ -125 +115 @@
 
       endif ! of if (is_master)
-      end
+      end subroutine inistats
+

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

@@ -1 +1 @@
       SUBROUTINE iniwrite(nid,idayref,phis)
 
-      USE comsoil_h
-
+      use comsoil_h, only: mlayer, nsoilmx
       IMPLICIT NONE
 
@@ -23 +22 @@
 #include "dimensions.h"
 #include "paramet.h"
-!include "comconst.h"
 #include "comcstfi.h"
 #include "comvert.h"
 #include "comgeom.h"
-#include "temps.h"
 #include "ener.h"
 #include "logic.h"
 #include "netcdf.inc"
 #include "serre.h"
-#include"dimphys.h"
 
 c   Arguments:
 c   ----------
 
-      integer nid        ! NetCDF file ID
-      INTEGER*4 idayref  ! date (initial date for this run)
-      REAL phis(ip1jmp1) ! surface geopotential
+      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
 
 c   Local:
@@ -57 +53 @@
          tab_cntrl(l)=0.
       ENDDO
-      tab_cntrl(1)  = FLOAT(iim)
-      tab_cntrl(2)  = FLOAT(jjm)
-      tab_cntrl(3)  = FLOAT(llm)
-      tab_cntrl(4)  = FLOAT(idayref)
+      tab_cntrl(1)  = real(iim)
+      tab_cntrl(2)  = real(jjm)
+      tab_cntrl(3)  = real(llm)
+      tab_cntrl(4)  = real(idayref)
       tab_cntrl(5)  = rad
       tab_cntrl(6)  = omeg
@@ -252 +248 @@
 
 !-------------------------------
-! (soil) depth variable mlayer() (known from comsoil.h)
+! (soil) depth variable mlayer() (known from comsoil_h)
 !-------------------------------
       ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode

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

@@ -1 +1 @@
 subroutine iniwritesoil(nid,ngrid)
-
-use comsoil_h, only : inertiedat, mlayer
-#ifdef CPP_PARA
-use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
-use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-#endif
 
 ! initialization routine for 'writediagoil'. Here we create/define
@@ -11 +5 @@
 ! (time-independent) parameters.
 
+use comsoil_h, only: mlayer, inertiedat, nsoilmx
+
 implicit none
 
 #include"dimensions.h"
-#include"dimphys.h"
 #include"paramet.h"
 #include"comcstfi.h"
@@ -38 +33 @@
 integer :: i,j,l,ig0
 
-#ifdef CPP_PARA
-! 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
-#else
-logical,parameter :: is_master=.true.
-logical,parameter :: is_mpi_root=.true.
-#endif
-
 ! 1. Define the dimensions
-if (is_master) then
 ! Switch to NetCDF define mode
 ierr=NF_REDEF(nid)
@@ -170 +155 @@
 ierr=NF_PUT_VAR_REAL(nid,varid,mlayer)
 #endif
-! Note mlayer(0:nsoilmx-1) known from comsoil.h
+! Note mlayer(0:nsoilmx-1) known from comsoil_h
 if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write depth variable"
@@ -197 +182 @@
 ! Note no need to write time variable here; it is done in writediagsoil.
 
-endif ! of if (is_master) 
-
 ! 3. Other variables to be included
 
 ! 3.1 mesh area surrounding each horizontal point
-if (is_master) then
 ierr=NF_REDEF(nid) ! switch to NetCDF define mode
 
@@ -236 +218 @@
 endif
 
-endif ! of if (is_master)
-
-
 ! 3.2 Thermal inertia
-if (is_master) then
 ierr=NF_REDEF(nid) ! switch to NetCDF define mode
 
@@ -262 +240 @@
 ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
 
-endif !of if (is_master)
-
 ! Recast data along 'dynamics' grid
-#ifdef CPP_PARA
-  ! gather field on a "global" (without redundant longitude) array
-  call Gather(inertiedat,dx3_glop)
-!$OMP MASTER
-  if (is_mpi_root) then
-    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,:,:)
-  endif
-!$OMP END MASTER
-!$OMP BARRIER
-#else
-! Note: inertiedat is known from comsoil.h
+! Note: inertiedat is known from comsoil_h
+
 do l=1,nsoilmx
   ! handle the poles
@@ -294 +258 @@
   enddo
 enddo ! of do l=1,nsoilmx
-#endif
-
-! Write data3 to file
-if (is_master) then
- 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)
-#endif
- if (ierr.ne.NF_NOERR) then
+
+! 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)
+#endif
+if (ierr.ne.NF_NOERR) then
   write(*,*)"iniwritesoil: Error, could not write th_inertia variable"
- endif
-endif ! of if (is_master)
+endif
 
 end subroutine iniwritesoil

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

@@ -1 +1 @@
 program kcm1d
 
+  use infotrac, only: nqtot
   use radinc_h,      only: NAERKIND
   use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
   use watercommon_h, only: mH2O
-  use ioipsl_getincom 
-  use comsaison_h
+  use ioipsl_getincom, only: getin
+  use comsaison_h, only: mu0, fract, dist_star
+!  use control_mod
 
   implicit none
@@ -32 +34 @@
 #include "comcstfi.h"
 #include "planete.h"
-#include "control.h"
+!#include "control.h"
 
   ! --------------------------------------------------------------
@@ -44 +46 @@
   real plev(nlayermx+1)   ! Intermediate pressure levels [Pa]
   real temp(nlayermx)     ! temperature at the middle of the layers [K]
-  real q(nlayermx,nqmx)   ! tracer mixing ratio [kg/kg]
-  real vmr(nlayermx,nqmx) ! tracer mixing ratio [mol/mol]
-  real qsurf(nqmx)        ! tracer surface budget [kg/kg] ????
+  real,allocatable :: q(:,:)   ! tracer mixing ratio [kg/kg]
+  real,allocatable :: vmr(:,:) ! tracer mixing ratio [mol/mol]
+  real,allocatable :: qsurf(:)        ! tracer surface budget [kg/kg] ????
   real psurf,psurf_n,tsurf      
 
@@ -77 +79 @@
   logical firstcall,lastcall,global1d
 
-  character*20 nametrac(nqmx)   ! name of the tracer (no need for adv trac common)
+  character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
 
   ! --------------
@@ -219 +221 @@
            write(*,*) "   (first line of traceur.def) "
            stop
-        else
-           ! check that the number of tracers is indeed nqmx
-           if (nq.ne.nqmx) then
-              write(*,*) "kcm1d: error, wrong number of tracers:"
-              write(*,*) "nq=",nq," whereas nqmx=",nqmx
-              stop
-           endif
         endif
+        nqtot=nq
+        ! allocate arrays which depend on number of tracers
+        allocate(nametrac(nq))
+        allocate(q(nlayermx,nq))
+        allocate(vmr(nlayermx,nq))
+        allocate(qsurf(nq))
 
         do iq=1,nq

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

@@ -3 +3 @@
      &     q,qsurf,surfith,nid)
 
-      USE surfdat_h
-      USE comsoil_h
-      USE tracer_h
+!      USE surfdat_h
+      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
+      USE tracer_h, ONLY: igcm_co2_ice
+      USE infotrac, ONLY: tname, nqtot
+!      USE control_mod
 
 c=======================================================================
@@ -26 +28 @@
 #include "dimensions.h"
 #include "dimphys.h"
-#include "planete.h"
+!#include "planete.h"
 #include "paramet.h"
 #include "comconst.h"
 #include "comvert.h"
 #include "comgeom2.h"
-#include "control.h"
-#include "logic.h"
+!#include "control.h"
+!#include "logic.h"
 #include "ener.h"
 #include "temps.h"
 #include "netcdf.inc"
-#include"advtrac.h"
+!#include"advtrac.h"
 c=======================================================================
 c   Declarations
@@ -79 +81 @@
       REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
       REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
-      REAL q(iip1,jjp1,llm,nqmx),qtot(iip1,jjp1,llm)
+      REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
 
 c autre variables dynamique nouvelle grille
@@ -97 +99 @@
       REAL co2ice(ngridmx) ! CO2 ice layer
       REAL emis(ngridmx)
-      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
+      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqtot)
 c     REAL phisfi(ngridmx)
 
@@ -119 +121 @@
       real co2iceS(iip1,jjp1)
       real emisS(iip1,jjp1)
-      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqmx)
+      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
 
       real ptotal, co2icetotal
@@ -283 +285 @@
       allocate(psold(imold+1,jmold+1))
       allocate(phisold(imold+1,jmold+1))
-      allocate(qold(imold+1,jmold+1,lmold,nqmx))
+      allocate(qold(imold+1,jmold+1,lmold,nqtot))
       allocate(co2iceold(imold+1,jmold+1))
       allocate(tsurfold(imold+1,jmold+1))
@@ -296 +298 @@
       allocate(surfithold(imold+1,jmold+1))
       allocate(mlayerold(nsoilold))
-      allocate(qsurfold(imold+1,jmold+1,nqmx))
+      allocate(qsurfold(imold+1,jmold+1,nqtot))
 
       allocate(var (imold+1,jmold+1,llm))
@@ -703 +705 @@
 
 ! Surface tracers:      
-      do iq=1,nqmx
+      do iq=1,nqtot
         ! initialize all surface tracers to zero
         call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
@@ -709 +711 @@
 
 
-!      print*,'tnom=',tnom
+!      print*,'tname=',tname
 !      print*,'nid',nid
 !      print*,'nvarid',nvarid
 !      stop
 
-      DO iq=1,nqmx
-          txt=trim(tnom(iq))//"_surf"
+      DO iq=1,nqtot
+          txt=trim(tname(iq))//"_surf"
           if (txt.eq."h2o_vap") then
             ! There is no surface tracer for h2o_vap;
@@ -748 +750 @@
         ENDIF
 
-      ENDDO ! of DO iq=1,nqmx
+      ENDDO ! of DO iq=1,nqtot
 
 
@@ -894 +896 @@
 
 ! Tracers:      
-      do iq=1,nqmx
+      do iq=1,nqtot
          call initial0((jmold+1)*(imold+1)*lmold,qold(1,1,1,iq) )
       enddo
 
-      DO iq=1,nqmx
-        txt=tnom(iq)
+      DO iq=1,nqtot
+        txt=tname(iq)
         write(*,*)"lect_start_archive: loading tracer ",trim(txt)
         ierr = NF_INQ_VARID (nid,txt,nvarid)
@@ -925 +927 @@
         ENDIF
 
-      ENDDO ! of DO iq=1,nqmx
+      ENDDO ! of DO iq=1,nqtot
 
 
@@ -1254 +1256 @@
 c     write(49,*) 'ucov',vcov
 
-      if (nqmx .gt. 0) then
+      if (nqtot .gt. 0) then
 c traceurs surface
-      do iq = 1, nqmx
+      do iq = 1, nqtot
             call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
      &                  imold,jmold,iim,jjm,1,
@@ -1262 +1264 @@
       enddo
 
-      call gr_dyn_fi (nqmx,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
+      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
 
 c traceurs 3D
-      do  iq = 1, nqmx
+      do  iq = 1, nqtot
             call interp_vert(qold(1,1,1,iq),var,lmold,llm,
      &        apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
@@ -1305 +1307 @@
 
 c     Periodicite :
-      do  iq = 1, nqmx
+      do  iq = 1, nqtot
          do l=1, llm
             do j = 1, jjp1
@@ -1316 +1318 @@
 ! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
 
-      endif !! if nqmx .ne. 0
+      endif !! if nqtot .ne. 0
 
 c-----------------------------------------------------------------------

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

@@ -1 +1 @@
 subroutine mkstats(ierr)
+
  
 !
@@ -9 +10 @@
 !  Yann W. july 2003
 
-#ifdef CPP_PARA
 use mod_phys_lmdz_para, only : is_master
-#endif
 
 implicit none
@@ -35 +34 @@
 integer :: meanid,sdid
 !integer, dimension(4) :: dimout
-#ifndef CPP_PARA
-logical,parameter :: is_master=.true.
-#endif
 
 ! Incrementation of count for the last step, which is not done in wstats

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

@@ -15 +15 @@
 c=======================================================================
 
-      USE tracer_h
-      USE comsoil_h
-      USE surfdat_h
-      USE comgeomfi_h
+      use infotrac, only: iniadvtrac, nqtot, tname
+      USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
+      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
+      USE surfdat_h, ONLY: phisfi, albedodat,
+     &                     zmea, zstd, zsig, zgam, zthe
+      USE comgeomfi_h, ONLY: lati, long, area
       use datafile_mod, only: datadir
 ! to use  'getin'
       USE ioipsl_getincom, only: getin
+      use control_mod, only: day_step, iphysiq, anneeref
+      use phyredem, only: physdem0, physdem1
+      use iostart, only: open_startphy
+      use comgeomphy, only: initcomgeomphy
       implicit none
 
@@ -31 +37 @@
 #include "comvert.h"
 #include "comgeom2.h"
-#include "control.h"
+!#include "control.h"
 #include "logic.h"
 #include "ener.h"
@@ -38 +44 @@
 #include "serre.h"
 #include "netcdf.inc"
-#include "advtrac.h"
+!#include "advtrac.h"
 c=======================================================================
 c   Declarations
@@ -64 +70 @@
       REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
       REAL phis(iip1,jjp1)
-      REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
+      REAL,ALLOCATABLE :: q(:,:,:,:)               ! champs advectes
 
 c autre variables dynamique nouvelle grille
@@ -90 +96 @@
       REAL emis(ngridmx)        ! surface emissivity
       real emisread             ! added by RW
-      REAL qsurf(ngridmx,nqmx)
+      REAL,ALLOCATABLE :: qsurf(:,:)
       REAL q2(ngridmx,nlayermx+1)
 !      REAL rnaturfi(ngridmx)
@@ -176 +182 @@
       pa     = 0. ! to ensure disaster rather than minor error if we don`t
                   ! make deliberate choice of these values elsewhere.
+
+! initialize "serial/parallel" related stuff
+      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+      call initcomgeomphy
+
+! Load tracer number and names:
+      call iniadvtrac(nqtot,numvanle)
+! allocate arrays
+      allocate(q(iip1,jjp1,llm,nqtot))
+      allocate(qsurf(ngridmx,nqtot))
 
 c=======================================================================
@@ -243 +259 @@
 c  INITIALISATIONS DIVERSES
 c=======================================================================
-! Load tracer names:
-      call iniadvtrac(nq,numvanle)
-      ! tnom(:) now contains tracer names
-! JL12 we will need the tracer names to read start in dyneta0
 
 ! Initialize global tracer indexes (stored in tracer.h)
 ! ... this has to be done before phyetat0
-      call initracer(ngridmx,nqmx,tnom)
+      call initracer(ngridmx,nqtot,tname)
 
 ! Take care of arrays in common modules
@@ -320 +332 @@
         write(*,*) 'Reading file START'
         fichnom = 'start.nc'
-        CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
+        CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
      .       ps,phis,time)
 
         write(*,*) 'Reading file STARTFI'
         fichnom = 'startfi.nc'
-        CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqmx,
+        CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqtot,
      .        day_ini,time,
      .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
@@ -387 +399 @@
 c-----------------------------------------------------------------------
       if (choix_1.eq.0) then
+         ! tabfi requires that input file be first opened by open_startphy(fichnom)
+         fichnom = 'start_archive.nc'
+         call open_startphy(fichnom)
          call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
      .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
       else if (choix_1.eq.1) then
+         fichnom = 'startfi.nc'
+         call open_startphy(fichnom)
          Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
          call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
@@ -434 +451 @@
       airefi(1)=airefi(1)*iim
       airefi(ngridmx)=airefi(ngridmx)*iim
+
+! also initialize various physics flags/settings which might be needed
+!    (for instance initracer needs to know about some flags, and/or
+!      'datafile' path may be changed by user)
+      call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
+     &                latfi,lonfi,airefi,rad,g,r,cpp)
 
 c=======================================================================
@@ -574 +597 @@
 
         read(*,fmt='(a20)') modif
-        if (modif(1:1) .eq. ' ') exit ! exit loop on changes
+        if (modif(1:1) .eq. ' ')then
+         exit ! exit loop on changes
+        endif
 
         write(*,*)
@@ -774 +799 @@
          if (yes.eq.'y') then
            write(*,*) 'OK : conservation of tracer total mass'
-           DO iq =1, nqmx
+           DO iq =1, nqtot
              DO l=1,llm
                DO j=1,jjp1
@@ -827 +852 @@
          do while (yes.eq.'y')
           write(*,*) 'Which tracer name do you want to change ?'
-          do iq=1,nqmx
-            write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
+          do iq=1,nqtot
+            write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
           enddo
           write(*,'(a35,i3)')
-     &            '(enter tracer number; between 1 and ',nqmx
+     &            '(enter tracer number; between 1 and ',nqtot
           write(*,*)' or any other value to quit this option)'
           read(*,*) iq
-          if ((iq.ge.1).and.(iq.le.nqmx)) then
-            write(*,*)'Change tracer name ',trim(tnom(iq)),' to ?'
+          if ((iq.ge.1).and.(iq.le.nqtot)) then
+            write(*,*)'Change tracer name ',trim(tname(iq)),' to ?'
             read(*,*) txt
-            tnom(iq)=txt
+            tname(iq)=txt
             write(*,*)'Do you want to change another tracer name (y/n)?'
             read(*,'(a)') yes 
@@ -843 +868 @@
 ! inapropiate value of iq; quit this option
             yes='n'
-          endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
+          endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
          enddo ! of do while (yes.ne.'y')
 
@@ -851 +876 @@
 c          mise a 0 des q (traceurs)
           write(*,*) 'Tracers set to 0 (1.E-30 in fact)'
-           DO iq =1, nqmx
+           DO iq =1, nqtot
              DO l=1,llm
                DO j=1,jjp1
@@ -862 +887 @@
 
 c          set surface tracers to zero
-           DO iq =1, nqmx
+           DO iq =1, nqtot
              DO ig=1,ngridmx
                  qsurf(ig,iq)=0.
@@ -872 +897 @@
         else if (trim(modif).eq.'q=x') then
              write(*,*) 'Which tracer do you want to modify ?'
-             do iq=1,nqmx
-               write(*,*)iq,' : ',trim(tnom(iq))
+             do iq=1,nqtot
+               write(*,*)iq,' : ',trim(tname(iq))
              enddo
-             write(*,*) '(choose between 1 and ',nqmx,')'
+             write(*,*) '(choose between 1 and ',nqtot,')'
              read(*,*) iq 
-             write(*,*)'mixing ratio of tracer ',trim(tnom(iq)),
+             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
      &                 ' ? (kg/kg)'
              read(*,*) val
@@ -887 +912 @@
                ENDDO
              ENDDO
-             write(*,*) 'SURFACE value of tracer ',trim(tnom(iq)),
+             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
      &                   ' ? (kg/m2)'
              read(*,*) val
@@ -942 +967 @@
              write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
              write(*,*) 'Which tracer do you want to set?'
-             do iq=1,nqmx
-               write(*,*)iq,' : ',trim(tnom(iq))
+             do iq=1,nqtot
+               write(*,*)iq,' : ',trim(tname(iq))
              enddo
-             write(*,*) '(choose between 1 and ',nqmx,')'
+             write(*,*) '(choose between 1 and ',nqtot,')'
              read(*,*) iq 
-             if ((iq.lt.1).or.(iq.gt.nqmx)) then
+             if ((iq.lt.1).or.(iq.gt.nqtot)) then
                ! wrong value for iq, go back to menu
                write(*,*) "wrong input value:",iq
@@ -953 +978 @@
              endif
              ! look for input file 'profile_tracer'
-             txt="profile_"//trim(tnom(iq))
+             txt="profile_"//trim(tname(iq))
              open(41,file=trim(txt),status='old',form='formatted',
      &            iostat=ierr)
@@ -970 +995 @@
                    q(:,:,l,iq)=profile(l+1)
                  enddo
-                 write(*,*)'OK, tracer ',trim(tnom(iq)),' initialized ',
-     &                     'using values from file ',trim(txt)
+                 write(*,*)'OK, tracer ',trim(tname(iq)),' initialized '
+     &                    ,'using values from file ',trim(txt)
                else
                  write(*,*)'problem reading file ',trim(txt),' !'
-                 write(*,*)'No modifications to tracer ',trim(tnom(iq))
+                 write(*,*)'No modifications to tracer ',trim(tname(iq))
                endif
              else
                write(*,*)'Could not find file ',trim(txt),' !'
-               write(*,*)'No modifications to tracer ',trim(tnom(iq))
+               write(*,*)'No modifications to tracer ',trim(tname(iq))
              endif
              
@@ -1504 +1529 @@
 
           
-      CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
-      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps) 
+      CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqtot)
+      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqtot,masse,ps) 
 C
 C Ecriture etat initial physique
 C
 
-!      do ig=1,ngridmx
-!         print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
-!         print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
-!      enddo
-!      stop
-
-
-      call physdem1(ngridmx,"restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
-     .              dtphys,float(day_ini),
-     .              time,tsurf,tsoil,emis,q2,qsurf,
-     .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
-     .              cloudfrac,totalfrac,hice,tnom)
+      call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
+     &              nqtot,dtphys,real(day_ini),0.0,
+     &              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
+      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nq,
+     &                dtphys,real(day_ini),
+     &                tsurf,tsoil,emis,q2,qsurf,
+     &                cloudfrac,totalfrac,hice)
 
 c=======================================================================
@@ -1534 +1554 @@
      *rage est differente de la valeur parametree llm =',i4//)
 
+      write(*,*) "newstart: All is well that ends well."
+
       end
 

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

@@ -1 @@
-      SUBROUTINE phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq,
-     .           day_ini,time,
-     .           tsurf,tsoil,emis,q2,qsurf,cloudfrac,totcloudfrac,hice)
-
-      USE surfdat_h
-      USE comgeomfi_h
-      USE tracer_h
-
-      implicit none
-
-c======================================================================
-c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
-c  Adaptation à Mars : Yann Wanherdrick 
-c Objet: Lecture de l etat initial pour la physique
-c======================================================================
+subroutine phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq, &
+                     day_ini,time,tsurf,tsoil, &
+                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice)
+
+  USE infotrac, ONLY: tname
+  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
+  use iostart, only: nid_start, open_startphy, close_startphy, &
+                     get_field, get_var, inquire_field, &
+                     inquire_dimension, inquire_dimension_length
+
+  implicit none
+
+!======================================================================
+! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
+!  Adaptation à Mars : Yann Wanherdrick 
+! Objet: Lecture de l etat initial pour la physique
+!======================================================================
 #include "netcdf.inc"
 #include "dimensions.h"
@@ -20 +22 @@
 #include "comcstfi.h"
 
-      INTEGER ngrid
-c======================================================================
-      INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
-      PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
+!======================================================================
+!  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
+!  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
 !======================================================================
 !  Arguments:
 !  ---------
 !  inputs:
-      character*(*) fichnom ! "startfi.nc" file
-      integer tab0
-      integer Lmodif
-      integer nsoil ! # of soil layers
-      integer nq
-      integer day_ini
-      real time
+  integer,intent(in) :: ngrid
+  character*(*),intent(in) :: fichnom ! "startfi.nc" file
+  integer,intent(in) :: tab0
+  integer,intent(in) :: Lmodif
+  integer,intent(in) :: nsoil ! # of soil layers
+  integer,intent(in) :: nq
+  integer,intent(in) :: day_ini
+  real,intent(in) :: time
 
 !  outputs:
-      real tsurf(ngrid,nbsrf) ! surface temperature
-      real tsoil(ngrid,nsoil,nbsrf) ! soil temperature
-      real emis(ngrid) ! surface emissivity
-      real q2(ngrid, llm+1) ! 
-      real qsurf(ngrid,nq) ! tracers on surface
-!      real co2ice(ngrid) ! co2 ice cover
-      real cloudfrac(ngrid,nlayermx)
-      real hice(ngrid), totcloudfrac(ngrid)
-
+  real,intent(out) :: tsurf(ngrid) ! surface temperature
+  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
+  real,intent(out) :: emis(ngrid) ! surface emissivity
+  real,intent(out) :: q2(ngrid, llm+1) ! 
+  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
+! real co2ice(ngrid) ! co2 ice cover
+  real,intent(out) :: cloudfrac(ngrid,nlayermx)
+  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
 
 !======================================================================
@@ -55 +56 @@
 
       real xmin,xmax ! to display min and max of a field
-c
+!
       INTEGER ig,iq,lmax
       INTEGER nid, nvarid
@@ -65 +66 @@
       CHARACTER*2 str2
       CHARACTER*1 yes
-c
+!
       REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
       INTEGER nqold
 
 ! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
-      logical :: oldtracernames=.false.
+!      logical :: oldtracernames=.false.
       integer :: count
       character(len=30) :: txt ! to store some text
+      
+      INTEGER :: indextime=1 ! index of selected time, default value=1
+      logical :: found
 
 !! added variables by AS to allow to read only slices of startfi
-      real :: toto(ngrid)
-      integer :: sta   !! subscript in starti where we start to read
-      integer, dimension(2) :: sta2d
-      integer, dimension(2) :: siz2d
-
-c AS: get the cursor that is stored in dimphys.h
-c ... this allows to read only a part of startfi horiz grid
-      sta = cursor
-
-c
-c ALLOCATE ARRAYS IN surfdat_h
-c
-      IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
-      IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
-      IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
-      IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
-      IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
-      IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
-      IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
-
-c
-c Ouvrir le fichier contenant l etat initial:
-c
-
-      ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
-      IF (ierr.NE.NF_NOERR) THEN
-        write(6,*)' Pb d''ouverture du fichier '//fichnom
-        CALL ABORT
-      ENDIF
-
-! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
-!                    qsurf02, ...)
-      count=0
-      do iq=1,nq
-        txt= " "
-        write(txt,'(a5,i2.2)')'qsurf',iq
-        ierr=NF_INQ_VARID(nid,txt,nvarid)
-        if (ierr.ne.NF_NOERR) then
-          ! did not find old tracer name
-          exit ! might as well stop here
-        else
-          ! found old tracer name
-          count=count+1
-        endif
-      enddo
-      if (count.eq.nq) then
-        write(*,*) "phyetat0:tracers seem to follow old naming ",
-     &             "convention (qsurf01,qsurf02,...)"
-        write(*,*) "   => will work for now ... "
-        write(*,*) "      but you should run newstart to rename them"
-        oldtracernames=.true.
-      endif
-
-c modifications possibles des variables de tab_cntrl
-      write(*,*) 'TABFI de phyeta0',Lmodif,tab0
-      call tabfi (ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
-     .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
-
-c
-c Lecture des latitudes (coordonnees):
-c
-      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <latitude> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
-         CALL abort
-      ENDIF
-c
-c Lecture des longitudes (coordonnees):
-c
-      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <longitude> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
-         CALL abort
-      ENDIF
-c
-c Lecture des aires des mailles:
-c
-      ierr = NF_INQ_VARID (nid, "area", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <area> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <area>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(area)
-      xmax = MAXVAL(area)
-      PRINT*,'Aires des mailles <area>:', xmin, xmax
-c
-c Lecture du geopotentiel au sol:
-c
-      ierr = NF_INQ_VARID (nid, "phisfi", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <phisfi> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,phisfi)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,phisfi)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <phisfi>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(phisfi)
-      xmax = MAXVAL(phisfi)
-      PRINT*,'Geopotentiel au sol <phisfi>:', xmin, xmax
-c
-c Lecture de l''albedo du sol nu:
-c
-      ierr = NF_INQ_VARID (nid, "albedodat", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <albedodat> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,albedodat)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,albedodat)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <albedodat>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(albedodat)
-      xmax = MAXVAL(albedodat)
-      PRINT*,'Albedo du sol nu <albedodat>:', xmin, xmax
-c
-c ZMEA
-c
-      ierr = NF_INQ_VARID (nid, "ZMEA", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zmea)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zmea)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <ZMEA>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      DO i = 1, ngrid
-         xmin = MIN(zmea(i),xmin)
-         xmax = MAX(zmea(i),xmax)
-      ENDDO
-      PRINT*,'<zmea>:', xmin, xmax
-c
-c ZSTD
-c
-      ierr = NF_INQ_VARID (nid, "ZSTD", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zstd)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zstd)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <ZSTD>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      DO i = 1, ngrid
-         xmin = MIN(zstd(i),xmin)
-         xmax = MAX(zstd(i),xmax)
-      ENDDO
-      PRINT*,'<zstd>:', xmin, xmax
-c
-c ZSIG
-c
-      ierr = NF_INQ_VARID (nid, "ZSIG", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zsig)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zsig)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <ZSIG>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      DO i = 1, ngrid
-         xmin = MIN(zsig(i),xmin)
-         xmax = MAX(zsig(i),xmax)
-      ENDDO
-      PRINT*,'<zsig>:', xmin, xmax
-c
-c ZGAM
-c
-      ierr = NF_INQ_VARID (nid, "ZGAM", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zgam)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zgam)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <ZGAM>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      DO i = 1, ngrid
-         xmin = MIN(zgam(i),xmin)
-         xmax = MAX(zgam(i),xmax)
-      ENDDO
-      PRINT*,'<zgam>:', xmin, xmax
-c
-c ZTHE
-c
-      ierr = NF_INQ_VARID (nid, "ZTHE", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zthe)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zthe)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <ZTHE>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      DO i = 1, ngrid
-         xmin = MIN(zthe(i),xmin)
-         xmax = MAX(zthe(i),xmax)
-      ENDDO
-      PRINT*,'<zthe>:', xmin, xmax
-c
-c CO2 ice cover
-c
-! Ehouarn: from now on, there is no "co2ice" standalone field; it is supposed
-! to be with stored in qsurf(i_co2_ice)
-!      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
-!      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Le champ <co2ice> est absent'
+!      real :: toto(ngrid)
+!      integer :: sta   !! subscript in starti where we start to read
+!      integer, dimension(2) :: sta2d
+!      integer, dimension(2) :: siz2d
+
+! AS: get the cursor that is stored in dimphys.h
+! ... this allows to read only a part of startfi horiz grid
+!sta = cursor
+
+!
+! ALLOCATE ARRAYS IN surfdat_h
+!
+IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
+IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
+IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
+IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
+IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
+IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
+IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
+
+! open physics initial state file:
+call open_startphy(fichnom)
+
+
+! possibility to modify tab_cntrl in tabfi
+write(*,*)
+write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
+call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
+                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
+
+!c
+!c Lecture des latitudes (coordonnees):
+!c
+!      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
+!      IF (ierr.NE.NF_NOERR) THEN
+!         PRINT*, 'phyetat0: Le champ <latitude> est absent'
 !         CALL abort
 !      ENDIF
 !#ifdef NC_DOUBLE
-!      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2ice)
+!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
 !#else
-!      ierr = NF_GET_VAR_REAL(nid, nvarid, co2ice)
+!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
 !#endif
 !      IF (ierr.NE.NF_NOERR) THEN
-!         PRINT*, 'phyetat0: Lecture echouee pour <co2ice>'
+!         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
+!         CALL abort
+!      ENDIF
+!c
+!c Lecture des longitudes (coordonnees):
+!c
+!      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
+!      IF (ierr.NE.NF_NOERR) THEN
+!         PRINT*, 'phyetat0: Le champ <longitude> est absent'
+!         CALL abort
+!      ENDIF
+!#ifdef NC_DOUBLE
+!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
+!#else
+!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
+!#endif
+!      IF (ierr.NE.NF_NOERR) THEN
+!         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
+!         CALL abort
+!      ENDIF
+!c
+!c Lecture des aires des mailles:
+!c
+!      ierr = NF_INQ_VARID (nid, "area", nvarid)
+!      IF (ierr.NE.NF_NOERR) THEN
+!         PRINT*, 'phyetat0: Le champ <area> est absent'
+!         CALL abort
+!      ENDIF
+!#ifdef NC_DOUBLE
+!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
+!#else
+!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
+!#endif
+!      IF (ierr.NE.NF_NOERR) THEN
+!         PRINT*, 'phyetat0: Lecture echouee pour <area>'
 !         CALL abort
 !      ENDIF
 !      xmin = 1.0E+20
 !      xmax = -1.0E+20
-!      xmin = MINVAL(co2ice)
-!      xmax = MAXVAL(co2ice)
-!      PRINT*,'CO2 ice cover <co2ice>:', xmin, xmax
-c
-c Lecture des temperatures du sol:
-c
-      ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <tsurf> est absent'
-         PRINT*, '          Mais je vais essayer de lire TS**'
-         IF (nbsrf.GT.99) THEN
-            PRINT*, "Trop de sous-mailles"
-            CALL abort
-         ENDIF
-         DO nsrf = 1, nbsrf
-           WRITE(str2,'(i2.2)') nsrf
-           ierr = NF_INQ_VARID (nid, "TS"//str2, nvarid)
-           IF (ierr.NE.NF_NOERR) THEN
-           PRINT*, "phyetat0: Le champ <TS"//str2//"> est absent"
-              CALL abort
-           ENDIF
-#ifdef NC_DOUBLE
-           ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsurf(1,nsrf))
-#else
-           ierr = NF_GET_VAR_REAL(nid, nvarid, tsurf(1,nsrf))
-#endif
-           IF (ierr.NE.NF_NOERR) THEN
-             PRINT*, "phyetat0: Lecture echouee pour <TS"//str2//">"
-             CALL abort
-           ENDIF
-           xmin = 1.0E+20
-           xmax = -1.0E+20
-           xmin = MINVAL(tsurf)
-           xmax = MAXVAL(tsurf)
-           PRINT*,'Temperature du sol TS**:', nsrf, xmin, xmax
-         ENDDO
-      ELSE
-         PRINT*, 'phyetat0: Le champ <tsurf> est present'
-         PRINT*, '          J ignore donc les autres temperatures TS**'
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, tsurf)
-#else
-         ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, tsurf)
-#endif
-         IF (ierr.NE.NF_NOERR) THEN
-            PRINT*, "phyetat0: Lecture echouee pour <TSURF>"
-            CALL abort
-         ENDIF
-         xmin = 1.0E+20
-         xmax = -1.0E+20
-         xmin = MINVAL(tsurf)
-         xmax = MAXVAL(tsurf)
-         PRINT*,'Temperature du sol <tsurf>', xmin, xmax
-         IF (nbsrf >= 2) THEN
-            DO nsrf = 2, nbsrf
-               DO i = 1, ngrid
-                  tsurf(i,nsrf) = tsurf(i,1)
-               ENDDO
-            ENDDO
-         ENDIF
-      ENDIF
-c
-c Lecture des temperatures du sol profond:
-c
-!      IF (nsoil.GT.99 .OR. nbsrf.GT.99) THEN
-!         PRINT*, "Trop de couches ou sous-mailles"
-!         CALL abort
-!      ENDIF
-!      DO nsrf = 1, nbsrf
-!         DO isoil=1, nsoil
-!            WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
-!            ierr = NF_INQ_VARID (nid, 'tsoil', nvarid)
-!            IF (ierr.NE.NF_NOERR) THEN
-!               PRINT*, "phyetat0: Le champ <tsoil> est absent"
-!               PRINT*, "          Il prend donc la valeur de surface"
-!               DO i=1, ngrid
-!                  tsoil(i,isoil,nsrf)=tsurf(i,nsrf)
-!               ENDDO
-!            ELSE
-!#ifdef NC_DOUBLE
-!              ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,1,nsrf))
-!#else
-!              ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,1,nsrf))
-!#endif
-!               IF (ierr.NE.NF_NOERR) THEN
-!                  PRINT*, "Lecture echouee pour <tsoil>"
-!                  CALL abort
-!               ENDIF
-!            ENDIF
-!         ENDDO
-!      ENDDO
-!      xmin = 1.0E+20
-!      xmax = -1.0E+20
-!      xmin = MINVAL(tsoil)
-!      xmax = MAXVAL(tsoil)
-!      PRINT*,'Temperatures du sol profond <tsoil>', xmin, xmax
-c
-c Surface emissivity
-c
-      ierr = NF_INQ_VARID (nid, "emis", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <emis> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta,ngrid, emis)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, emis)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <emis>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(emis)
-      xmax = MAXVAL(emis)
-      PRINT*,'Surface emissivity <emis>:', xmin, xmax
-
-c
-c Cloud fraction (added by BC 2010)
-c
-      ierr = NF_INQ_VARID (nid, "cloudfrac", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <cloudfrac> est absent'
-      cloudfrac(:,:)=0.5
-!         CALL abort
-      ENDIF
-      sta2d = (/sta,1/)
-      siz2d = (/ngrid, nlayermx/)
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,cloudfrac)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,cloudfrac)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <cloudfrac>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(cloudfrac)
-      xmax = MAXVAL(cloudfrac)
-      PRINT*,'Cloud fraction <cloudfrac>:', xmin, xmax
-
-
-c
-c Total cloud fraction (added by BC 2010)
-c
-      ierr = NF_INQ_VARID (nid, "totcloudfrac", nvarid)
-!      ierr = NF_INQ_VARID (nid, "totalfrac", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <totcloudfrac> est absent'
-      totcloudfrac(:)=0.5
-!         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, totcloudfrac)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, totcloudfrac)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <totcloudfrac>'
-         !CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(totcloudfrac)
-      xmax = MAXVAL(totcloudfrac)
-      PRINT*,'Cloud fraction <totcloudfrac>:', xmin, xmax
-
-
-
-c
-c Height of oceanic ice (added by BC 2010)
-c
-      ierr = NF_INQ_VARID (nid, "hice", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <hice> est absent'
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, hice)
-#else
-      ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, hice)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <hice>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(hice)
-      xmax = MAXVAL(hice)
-      PRINT*,'Height of oceanic ice <hice>:', xmin, xmax
-
-
-c
-c pbl wind variance
-c
-      ierr = NF_INQ_VARID (nid, "q2", nvarid)
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Le champ <q2> est absent'
-         CALL abort
-      ENDIF
-      sta2d = (/sta,1/)
-      siz2d = (/ngrid, llm+1/)
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,q2)
-#else
-      ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,q2)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'phyetat0: Lecture echouee pour <q2>'
-         CALL abort
-      ENDIF
-      xmin = 1.0E+20
-      xmax = -1.0E+20
-      xmin = MINVAL(q2)
-      xmax = MAXVAL(q2)
-      PRINT*,'pbl wind variance <q2>:', xmin, xmax
-c
-c tracer on surface
-c
-
-      IF(nq.GE.1) THEN
-         nqold=nq
-         DO iq=1,nq
-!            str7(1:5)='qsurf'
-!            WRITE(str7(6:7),'(i2.2)') iq
-!            ierr = NF_INQ_VARID (nid,str7,nvarid)
-           IF (oldtracernames) THEN
-             txt=" "
-             write(txt,'(a5,i2.2)')'qsurf',iq
-           ELSE
-             txt=noms(iq)
-!             if (txt.eq."h2o_vap") then
-               ! There is no surface tracer for h2o_vap;
-               ! "h2o_ice" should be loaded instead
-!               txt="h2o_ice"
-!               write(*,*) 'phyetat0: loading surface tracer',
-!     &                     ' h2o_ice instead of h2o_vap'
-!             endif
-           ENDIF ! of IF (oldtracernames) THEN
-           ierr=NF_INQ_VARID(nid,txt,nvarid)
-           IF (ierr.NE.NF_NOERR) THEN
-             write(*,*) 'PHYETAT0: WARNING : surface tracer',trim(txt),
-     &                  ' not found in file'
-             write(*,*) trim(txt), ' set to 0'
-             do ig=1,ngrid
-               qsurf(ig,iq)=0.
-             end do
-             nqold=min(iq-1,nqold)
-           ELSE
-           toto(:) = 0.
-#ifdef NC_DOUBLE
-           ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,toto)
-#else
-           ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,toto)
-#endif
-             IF (ierr.NE.NF_NOERR) THEN
-               PRINT*, 'phyetat0: Lecture echouee pour <',trim(txt),'>'
-               CALL abort
-             ENDIF
-             qsurf(1:ngrid,iq) = toto(1:ngrid)
-           ENDIF
-           xmin = 1.0E+20
-           xmax = -1.0E+20
-           xmin = MINVAL(qsurf(1:ngrid,iq))
-           xmax = MAXVAL(qsurf(1:ngrid,iq))
-           PRINT*,'tracer on surface <',trim(txt),'>:',xmin,xmax
-         ENDDO
-         if ((nqold.lt.nq).and.(nqold.ge.1)) then
-c        case when new tracer are added in addition to old ones
-             write(*,*)'qsurf 1 to ', nqold,'were already present'
-             write(*,*)'qsurf ', nqold+1,' to ', nq,'are new'
-             write(*,*)' and initialized to zero'
-             qsurf(:,nqold+1:nq)=0.0
-!            yes=' '
-!            do while ((yes.ne.'y').and.(yes.ne.'n'))
-!             write(*,*) 'Would you like to reindex qsurf # 1 ->',nqold
-!             write(*,*) 'to #',nq-nqold+1,'->', nq,'   (y or n) ?'
-!             read(*,fmt='(a)') yes
-!            end do
-!            if (yes.eq.'y') then
-!              write(*,*) 'OK, let s reindex qsurf'
-!                 do ig=1,ngrid
-!                    do iq=nq,nq-nqold+1,-1
-!                       qsurf(ig,iq)=qsurf(ig,iq-nq+nqold)
-!                    end do
-!                    do iq=nq-nqold,1,-1
-!                       qsurf(ig,iq)= 0.
-!                    end do
-!                 end do
-!            end if
-         end if ! of if ((nqold.lt.nq).and.(nqold.ge.1))
-      ENDIF ! of IF(nq.GE.1)
+!      xmin = MINVAL(area)
+!      xmax = MAXVAL(area)
+!      PRINT*,'Aires des mailles <area>:', xmin, xmax
+
+! Load surface geopotential:
+call get_field("phisfi",phisfi,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <phisfi>"
+  call abort
+else
+  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
+             minval(phisfi), maxval(phisfi)
+endif
+
+! Load bare ground albedo:
+call get_field("albedodat",albedodat,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <albedodat>"
+  call abort
+else
+  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
+             minval(albedodat), maxval(albedodat)
+endif
+
+! ZMEA
+call get_field("ZMEA",zmea,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZMEA>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZMEA> range:", &
+             minval(zmea), maxval(zmea)
+endif
+
+! ZSTD
+call get_field("ZSTD",zstd,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZSTD>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZSTD> range:", &
+             minval(zstd), maxval(zstd)
+endif
+
+! ZSIG
+call get_field("ZSIG",zsig,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZSIG>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZSIG> range:", &
+             minval(zsig), maxval(zsig)
+endif
+
+! ZGAM
+call get_field("ZGAM",zgam,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZGAM>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZGAM> range:", &
+             minval(zgam), maxval(zgam)
+endif
+
+! ZTHE
+call get_field("ZTHE",zthe,found)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <ZTHE>"
+  call abort
+else
+  write(*,*) "phyetat0: <ZTHE> range:", &
+             minval(zthe), maxval(zthe)
+endif
+
+! Surface temperature :
+call get_field("tsurf",tsurf,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <tsurf>"
+  call abort
+else
+  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
+             minval(tsurf), maxval(tsurf)
+endif
+
+! Surface emissivity
+call get_field("emis",emis,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <emis>"
+  call abort
+else
+  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
+             minval(emis), maxval(emis)
+endif
+
+! Cloud fraction (added by BC 2010)
+call get_field("cloudfrac",cloudfrac,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <cloudfrac>"
+  call abort
+else
+  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
+             minval(cloudfrac), maxval(cloudfrac)
+endif
+
+! Total cloud fraction (added by BC 2010)
+call get_field("totcloudfrac",totcloudfrac,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <totcloudfrac>"
+  call abort
+else
+  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
+             minval(totcloudfrac), maxval(totcloudfrac)
+endif
+
+! Height of oceanic ice (added by BC 2010)
+call get_field("hice",hice,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <hice>"
+  call abort
+else
+  write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
+             minval(hice), maxval(hice)
+endif
+
+! pbl wind variance
+call get_field("q2",q2,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: Failed loading <q2>"
+  call abort
+else
+  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
+             minval(q2), maxval(q2)
+endif
+
+! tracer on surface
+if (nq.ge.1) then
+  do iq=1,nq
+    txt=tname(iq)
+    if (txt.eq."h2o_vap") then
+      ! There is no surface tracer for h2o_vap;
+      ! "h2o_ice" should be loaded instead
+      txt="h2o_ice"
+      write(*,*) 'phyetat0: loading surface tracer', &
+                           ' h2o_ice instead of h2o_vap'
+    endif
+    call get_field(txt,qsurf(:,iq),found,indextime)
+    if (.not.found) then
+      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
+      write(*,*) "         ",trim(txt)," is set to zero"
+    else
+      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
+                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
+    endif
+  enddo
+endif ! of if (nq.ge.1)
 
 ! Call to soil_settings, in order to read soil temperatures,
 ! as well as thermal inertia and volumetric heat capacity
 
-      PRINT*,'SOIL INIT'
-      call soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
-c
-c Fermer le fichier:
-c
-      ierr = NF_CLOSE(nid)
-c
-
-      RETURN
-      END
+call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
+
+!
+! close file:
+!
+call close_startphy
+
+END SUBROUTINE phyetat0

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

@@ -1 @@
-      subroutine physdem1(ngrid,filename,lonfi,latfi,nsoil,nq,
-     .                   phystep,day_ini,
-     .                   time,tsurf,tsoil,emis,q2,qsurf,
-     .                   airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe,
-     .                   cloudfrac,totcloudfrac,hice,nametrac)
-
-      use surfdat_h
-      use comsoil_h
-      USE comgeomfi_h
- 
-      implicit none
-
-!==================================================================
-!     
-!     Purpose
-!     -------
-!     Create physics (re-)start data file "restartfi.nc"
-!     
-!     Called by
-!     ---------
-!     rcm1d.F90
-!     physiq.F90
-!     newstart.F
-!     
-!     Calls
-!     -----
-!     none
-!
-!==================================================================
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comvert.h"
-#include "comgeom2.h"
-#include "control.h"
-#include "comdissnew.h"
-#include "logic.h"
-#include "ener.h"
-#include "netcdf.inc"
-#include "dimphys.h"
-#include"callkeys.h"
-
-      INTEGER :: ngrid
-
-      INTEGER nid,iq
-      INTEGER, parameter :: ivap=1
-      REAL, parameter :: qsolmax= 150.0
-      character (len=*) :: filename
-      character (len=7) :: str7
-
-      REAL day_ini
-      INTEGER nsoil,nq
-      integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid
-
-      REAL phystep,time
-      REAL latfi(ngrid), lonfi(ngrid)
-!      REAL champhys(ngrid)
-      REAL tsurf(ngrid)
-      INTEGER length
-      PARAMETER (length=100)
-      REAL tab_cntrl(length)
-
-!     added by BC
-      REAL hice(ngrid),cloudfrac(ngrid,nlayermx)
-      REAL totcloudfrac(ngrid)
-
-!     AS: name of tracers added as an argument to avoid using nqmx in commons (i.e. adv trac)
-!       nota: physdem1 is used both in physiq and newstart.
-      character*20 nametrac(nq)   ! name of the tracer
-
-
-!      EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
-!      EXTERNAL exner_hyb , SSUM
-
-#include "serre.h"
-#include "fxyprim.h"
+module phyredem
+
+implicit none
+
+contains
+
+subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
+                         phystep,day_ini,time,airefi, &
+                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
+! create physics restart file and write time-independent variables
+  use comsoil_h, only: inertiedat, volcapa, mlayer
+  use comgeomfi_h, only: area
+  use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, &
+                       albedice, emisice, emissiv, &
+                       iceradius, dtemisice, phisfi
+  use iostart, only : open_restartphy, close_restartphy, & 
+                      put_var, put_field, length
+  use mod_grid_phy_lmdz, only : klon_glo
+
+  implicit none
 #include "planete.h"
 #include "comcstfi.h"
-
-      real tsoil(ngrid,nsoil),emis(ngrid)
-      real q2(ngrid, llm+1),qsurf(ngrid,nq)
-      real airefi(ngrid)
-      real alb(ngrid),ith(ngrid,nsoil)
-      real pzmea(ngrid),pzstd(ngrid)
-      real pzsig(ngrid),pzgam(ngrid),pzthe(ngrid)
-      integer ig
-
-! flag which identifies if we are using old tracer names (qsurf01,...)
-      logical :: oldtracernames=.false.
-      integer :: count
-      character(len=30) :: txt ! to store some text
-! indexes of water vapour & water ice tracers (if any):
-      integer :: i_h2o_vap=0
-      integer :: i_h2o_ice=0
-c-----------------------------------------------------------------------
-
-      ! copy airefi(:) to area(:)
-      CALL SCOPY(ngrid,airefi,1,area,1)
-      ! note: area() is defined in comgeomfi.h
-
-      DO ig=1,ngrid
-         albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat.h
-         zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat.h
-         zstd(ig)=pzstd(ig) ! note: zstd() is defined in surfdat.h
-         zsig(ig)=pzsig(ig) ! note: zsig() is defined in surfdat.h
-         zgam(ig)=pzgam(ig) ! note: zgam() is defined in surfdat.h
-         zthe(ig)=pzthe(ig) ! note: zthe() is defined in surfdat.h
-      ENDDO
-
-      inertiedat(:,:)=ith(:,:) ! note inertiedat() is defined in comsoil.h
-c
-c  things to store in the physics start file:
-c
-      ierr = NF_CREATE(adjustl(filename),NF_CLOBBER, nid)
-      IF (ierr.NE.NF_NOERR) THEN
-        WRITE(6,*)'physdem1: Problem creating file ',adjustl(filename)
-        write(6,*) NF_STRERROR(ierr)
-        CALL ABORT
-      ENDIF
-c
-      print*,'we got this far'
-
-      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 18,
-     .                       "Physics start file")
-c
-      ierr = NF_DEF_DIM (nid,"index",length,idim1)
-      if (ierr.ne.NF_NOERR) then
-        WRITE(6,*)'physdem1: Problem defining index dimension'
-        write(6,*) NF_STRERROR(ierr)
-        call abort
-      endif
-c
-      ierr = NF_DEF_DIM (nid,"physical_points",ngrid,idim2)
-      if (ierr.ne.NF_NOERR) then
-        WRITE(6,*)'physdem1: Problem defining physical_points dimension'
-        write(6,*) NF_STRERROR(ierr)
-        call abort
-      endif
-c
-      ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoil,idim3)
-      if (ierr.ne.NF_NOERR) then
-      WRITE(6,*)'physdem1: Problem defining subsurface_layers dimension'
-        write(6,*) NF_STRERROR(ierr)
-        call abort
-      endif
-c
-      ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
-      if (ierr.ne.NF_NOERR) then
-        WRITE(6,*)'physdem1: Problem defining nlayer+1 dimension'
-        write(6,*) NF_STRERROR(ierr)
-        call abort
-      endif
-c
-      ierr = NF_DEF_DIM (nid,"number_of_advected_fields",nq,idim5)
-      if (ierr.ne.NF_NOERR) then
-        WRITE(6,*)'physdem1: Problem defining advected fields dimension'
-        WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
-        write(6,*) NF_STRERROR(ierr)
-      endif
-
-      ierr = NF_DEF_DIM (nid,"nlayer",llm,idim6)
-      if (ierr.ne.NF_NOERR) then
-        WRITE(6,*)'physdem1: Problem defining nlayer dimension'
-        write(6,*) NF_STRERROR(ierr)
-        call abort
-      endif
-
-      ierr = NF_ENDDEF(nid) ! exit NetCDF define mode
-
-c clear tab_cntrl(:) array
-      DO ierr = 1, length
-         tab_cntrl(ierr) = 0.0
-      ENDDO
-
-      write(*,*) "physdem1: ngrid: ",ngrid
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c Fill control array tab_cntrl(:) with paramleters for this run
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-c Informations on the physics grid
-      tab_cntrl(1) = float(ngrid)  ! number of nodes on physics grid
-      tab_cntrl(2) = float(nlayermx) ! number of atmospheric layers
-      tab_cntrl(3) = day_ini + int(time)         ! initial day 
-      tab_cntrl(4) = time -int(time)            ! initiale time of day
-
-c Informations about Mars, used by dynamics and physics
-      tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
-      tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
-      tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
-      tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
-      tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
-      tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
-
-      tab_cntrl(11) = phystep  ! time step in the physics
-      tab_cntrl(12) = 0.
-      tab_cntrl(13) = 0.
-
-c Informations about Mars, only for physics
-      tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
-      tab_cntrl(15) = periastr  ! min. star-planet distance (AU)
-      tab_cntrl(16) = apoastr   ! max. star-planet distance (AU) 
-      tab_cntrl(17) = peri_day  ! date of periastron (sols since N. spring)
-      tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98
-
-c Boundary layer and turbulence
-      tab_cntrl(19) = z0        ! surface roughness (m) ~0.01
-      tab_cntrl(20) = lmixmin   ! mixing length ~100
-      tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
-
-c Optical properties of polar caps and ground emissivity
-      tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
-      tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
-      tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
-      tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
-      tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
-      tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
-      tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
-      tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
-      tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
-
-      tab_cntrl(28) = 0. 
-      tab_cntrl(29) = 0.
-      tab_cntrl(30) = 0.
+  character(len=*), intent(in) :: filename
+  real,intent(in) :: lonfi(ngrid)
+  real,intent(in) :: latfi(ngrid)
+  integer,intent(in) :: nsoil
+  integer,intent(in) :: ngrid
+  integer,intent(in) :: nlay
+  integer,intent(in) :: nq
+  real,intent(in) :: phystep
+  real,intent(in) :: day_ini
+  real,intent(in) :: time
+  real,intent(in) :: airefi(ngrid)
+  real,intent(in) :: alb(ngrid)
+  real,intent(in) :: ith(ngrid,nsoil)
+  real,intent(in) :: pzmea(ngrid)
+  real,intent(in) :: pzstd(ngrid)
+  real,intent(in) :: pzsig(ngrid)
+  real,intent(in) :: pzgam(ngrid)
+  real,intent(in) :: pzthe(ngrid)
+ 
+  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
+  
+  ! Create physics start file
+  call open_restartphy(filename)
+
+! tab_cntrl() contains run parameters
+  tab_cntrl(:)=0 ! initialization
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! Fill control array tab_cntrl(:) with paramleters for this run
+!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+! Informations on the physics grid
+  tab_cntrl(1) = float(klon_glo)  ! number of nodes on physics grid
+  tab_cntrl(2) = float(nlay) ! number of atmospheric layers
+  tab_cntrl(3) = day_ini + int(time)         ! final day 
+  tab_cntrl(4) = time -int(time)            ! final time of day
+
+! Informations about Mars, used by dynamics and physics
+  tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
+  tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
+  tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
+  tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
+  tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
+  tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
+
+  tab_cntrl(11) = phystep  ! time step in the physics
+  tab_cntrl(12) = 0.
+  tab_cntrl(13) = 0.
+
+! Informations about Mars, only for physics
+  tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
+  tab_cntrl(15) = periastr  ! min. star-planet distance (AU)
+  tab_cntrl(16) = apoastr   ! max. star-planet distance (AU) 
+  tab_cntrl(17) = peri_day  ! date of periastron (sols since N. spring)
+  tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98
+
+! Boundary layer and turbulence
+  tab_cntrl(19) = z0        ! surface roughness (m) ~0.01
+  tab_cntrl(20) = lmixmin   ! mixing length ~100
+  tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
+
+! Optical properties of polar caps and ground emissivity
+  tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
+  tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
+  tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
+  tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
+  tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
+  tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
+  tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
+  tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
+  tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
+
+  tab_cntrl(28) = 0. 
+  tab_cntrl(29) = 0.
+  tab_cntrl(30) = 0.
 
 ! Soil properties:
-      tab_cntrl(35) = volcapa ! soil volumetric heat capacity
-      
-
-!      write(*,*) "physdem1: tab_cntrl():",tab_cntrl
-      
-      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
-         CALL abort
-      ENDIF
-      ! define variable
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'physdem1: Failed to define controle'
-         CALL abort
-      ENDIF
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
-     .                        "Control parameters")
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'physdem1: Failed to define controle title attribute'
-         CALL abort
-      ENDIF
-      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'physdem1: Failed to swich out of NetCDF define mode'
-         CALL abort
-      ENDIF
-      ! write variable
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-         PRINT*, 'physdem1: Failed to write controle data'
-         CALL abort
-      ENDIF
-
-! write mid-layer depths mlayer() !known from comsoil.h
-
-      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
-      ! define variable
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid,"soildepth",NF_DOUBLE,1,idim3,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid,"soildepth",NF_FLOAT,1,idim3,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
-     .                        "Soil mid-layer depth")
-      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
-      ! write variable
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,mlayer)
-#endif
-
-c
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
-     .                        "Longitudes of physics grid")
-      ierr = NF_ENDDEF(nid)
-
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
-#endif
-
-c
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
-     .                        "Latitudes of physics grid")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
-#endif
-
-c
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "area", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "area", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
-     .                        "Mesh area")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
-#endif
-
-c
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "phisfi", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "phisfi", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
-     .                        "Geopotential at the surface")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phisfi)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,phisfi)
-#endif
-
-c
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "albedodat", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "albedodat", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
-     .                        "Albedo of bare ground")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedodat)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,albedodat)
-#endif
-
-c
-c   some data for Francois Lott's programs
-c
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
-     .                        "Relief: mean relief")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
-#endif
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
-     .                        "Relief: standard deviation")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
-#endif
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
-     .                        "Relief: sigma parameter")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
-#endif
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
-     .                        "Relief: gamma parameter")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
-#endif
-c
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
-     .                        "Relief: theta parameter")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
-#endif
-
-c Write the physical fields
-
-
-! Soil Thermal inertia
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid,"inertiedat",NF_DOUBLE,
-     &                   2,(/idim2,idim3/),nvarid)
-#else
-      ierr = NF_DEF_VAR (nid,"inertiedat",NF_FLOAT,
-     &                   2,(/idim2,idim3/),nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 20,
-     .                        "Soil thermal inertia")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,inertiedat)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,inertiedat)
-#endif
-
-
-
-
-!  Surface temperature
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "tsurf", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "tsurf", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
-     .                        "Surface temperature")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsurf)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsurf)
-#endif
-
-! Soil temperature
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid,"tsoil",NF_DOUBLE,2,(/idim2,idim3/),nvarid)
-#else
-!      ierr = NF_DEF_VAR (nid, "tsoil", NF_FLOAT, 2, idim2,nvarid)
-      ierr = NF_DEF_VAR (nid,"tsoil",NF_FLOAT,2,(/idim2,idim3/),nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 16,
-     .                        "Soil temperature")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil)
-#endif
-
-c emissivity
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "emis", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "emis", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 18,
-     .                        "Surface emissivity")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,emis)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,emis)
-#endif
-
-c planetary boundary layer
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid,"q2", NF_DOUBLE, 2, (/idim2,idim4/),nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "q2", NF_FLOAT, 2,(/idim2,idim4/),nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
-     .                        "pbl wind variance")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q2)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,q2)
-#endif
-      IF (ierr.NE.NF_NOERR) THEN
-        PRINT*, 'physdem1: Failed to write q2'
-        CALL abort
-      ENDIF
-
-
-
-
-
-
-
-c cloud fraction (added by BC 2010)
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 2,
-     & (/idim2,idim6/),nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 2,
-     & (/idim2,idim6/),nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
-     .                        "Cloud fraction")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cloudfrac)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,cloudfrac)
-#endif
-
-c total cloud fraction (added by BC 2010)
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid,"totcloudfrac", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "totcloudfrac", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
-     .                        "Total fraction")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,totcloudfrac)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,totcloudfrac)
-#endif
-
-c oceanic ice (added by BC 2010)
-
-      ierr = NF_REDEF (nid)
-#ifdef NC_DOUBLE
-      ierr = NF_DEF_VAR (nid, "hice", NF_DOUBLE, 1, idim2,nvarid)
-#else
-      ierr = NF_DEF_VAR (nid, "hice", NF_FLOAT, 1, idim2,nvarid)
-#endif
-      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
-     .                        "Height of oceanic ice")
-      ierr = NF_ENDDEF(nid)
-#ifdef NC_DOUBLE
-      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,hice)
-#else
-      ierr = NF_PUT_VAR_REAL (nid,nvarid,hice)
-#endif
-
-
-
-
-
-
-
-
-
-
-
-c tracers
-
-      if (nq > 0) then
-
-! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
-!                    qsurf02, ...)
-      count=0
-      do iq=1,nq
-        txt= " "
-        write(txt,'(a1,i2.2)')'q',iq
-        if (txt.ne.nametrac(iq)) then ! use tracer names stored in dynamics
-          ! did not find old tracer name
-          exit ! might as well stop here
-        else
-          ! found old tracer name
-          count=count+1
-        endif
-      enddo
-      if (count.eq.nq) then
-        write(*,*) "physdem1:tracers seem to follow old naming ",
-     &             "convention (qsurf01,qsurf02,...)"
-
-        call abort
-        !write(*,*) "   => will work for now ... "
-        !write(*,*) "      but you should run newstart to rename them"
-        !oldtracernames=.true.
-        ! Moreover, if computing water cycle with ice, move surface ice
-        ! back to qsurf(nq)
-        !IF (iceparty) THEN
-        !  write(*,*)'physdem1: moving surface water ice to index ',nq
-        !  qsurf(1:ngrid,nq)=qsurf(1:ngrid,nq-1)
-        !  qsurf(1:ngrid,nq-1)=0
-        !ENDIF
-      endif ! of if (count.eq.nq)
-
-      IF(nq.GE.1) THEN
-! preliminary stuff: look for water vapour & water ice tracers (if any)
-        if (.not.oldtracernames) then
-         do iq=1,nq
-           if (nametrac(iq).eq."h2o_vap") then
-             i_h2o_vap=iq
-           endif
-           if (nametrac(iq).eq."h2o_ice") then
-             i_h2o_ice=iq
-           endif
-         enddo ! of iq=1,nq
-         ! handle special case of only water vapour tracer (no ice)
-         if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then
-          ! then the index of (surface) ice is i_h2o_vap
-            print*,'water vapour but no water ice, WTF?'
-            call abort
-            i_h2o_ice=i_h2o_vap
-         endif
-        endif ! of if (.not.oldtracernames) 
-
-         DO iq=1,nq
-           IF (oldtracernames) THEN
-             txt=" "
-             write(txt,'(a5,i2.2)')'qsurf',iq
-           ELSE
-             txt=nametrac(iq)
-
-
-!     in new version, h2o_vap acts as liquid surface tracer
-!     so the section below is now redundant
-!     ------------------------------------------------------------------
-!             ! Exception: there is no water vapour surface tracer
-!             if (txt.eq."h2o_vap") then
-!               write(*,*)"physdem1: skipping water vapour tracer"
-!               if (i_h2o_ice.eq.i_h2o_vap) then
-!               ! then there is no "water ice" tracer; but still
-!               ! there is some water ice on the surface
-!                 write(*,*)"          writting water ice instead"
-!                 txt="h2o_ice"
-!               else
-!               ! there is a "water ice" tracer which has been / will be
-!               ! delt with in due time
-!                 cycle
-!               endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
-!             endif ! of if (txt.eq."h2o_vap")
-!     ------------------------------------------------------------------
-
-
-
-           ENDIF ! of IF (oldtracernames)
-
-           ierr=NF_REDEF(nid)
-           IF (ierr.NE.NF_NOERR) THEN
-             PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
-             CALL abort
-           ENDIF
-#ifdef NC_DOUBLE
-           ierr=NF_DEF_VAR(nid,txt,NF_DOUBLE,1,idim2,nvarid)
-#else
-           ierr=NF_DEF_VAR(nid,txt,NF_FLOAT,1,idim2,nvarid)
-#endif
-           IF (ierr.NE.NF_NOERR) THEN
-             PRINT*, 'physdem1: Failed to define ',trim(txt)
-             CALL abort
-           ENDIF
-           ierr=NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
-     &                        "tracer on surface")
-           IF (ierr.NE.NF_NOERR) THEN
-             PRINT*, 'physdem1: Failed to define ',trim(txt),
-     &               ' title attribute'
-             CALL abort
-           ENDIF
-           ierr=NF_ENDDEF(nid)
-           IF (ierr.NE.NF_NOERR) THEN
-             PRINT*, 'physdem1: Failed to swich out of define mode'
-             CALL abort
-           ENDIF
-           
-#ifdef NC_DOUBLE
-            ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,qsurf(1,iq))
-#else
-            ierr=NF_PUT_VAR_REAL (nid,nvarid,qsurf(1,iq))
-#endif
-           IF (ierr.NE.NF_NOERR) THEN
-             PRINT*, 'physdem1: Failed to write ',trim(txt)
-             CALL abort
-           ENDIF
-         ENDDO ! of DO iq=1,nq
-      ENDIF ! of IF(nq.GE.1)
-
-      endif ! of if tracer
-
-c close file
-      ierr = NF_CLOSE(nid)
-
-      RETURN
-
-      END
+  tab_cntrl(35) = volcapa ! soil volumetric heat capacity
+
+  call put_var("controle","Control parameters",tab_cntrl)
+  
+  ! Write the mid-layer depths
+  call put_var("soildepth","Soil mid-layer depth",mlayer)
+  
+  ! Write longitudes
+  call put_field("longitude","Longitudes of physics grid",lonfi)
+  
+  ! Write latitudes
+  call put_field("latitude","Latitudes of physics grid",latfi)
+  
+  ! Write mesh areas
+  call put_field("area","Mesh area",area)
+  
+  ! Write surface geopotential
+  call put_field("phisfi","Geopotential at the surface",phisfi)
+  
+  ! Write surface albedo
+  call put_field("albedodat","Albedo of bare ground",albedodat)
+  
+  ! Subgrid topogaphy variables
+  call put_field("ZMEA","Relief: mean relief",zmea)
+  call put_field("ZSTD","Relief: standard deviation",zstd)
+  call put_field("ZSIG","Relief: sigma parameter",zsig)
+  call put_field("ZGAM","Relief: gamma parameter",zgam)
+  call put_field("ZTHE","Relief: theta parameter",zthe)
+  
+  ! Soil thermal inertia
+  call put_field("inertiedat","Soil thermal inertia",inertiedat)
+  
+  ! Close file
+  call close_restartphy
+  
+end subroutine physdem0
+
+subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
+                    phystep,time,tsurf,tsoil,emis,q2,qsurf, &
+                    cloudfrac,totcloudfrac,hice)
+  ! write time-dependent variable to restart file
+  use iostart, only : open_restartphy, close_restartphy, & 
+                      put_var, put_field
+  use infotrac, only: tname
+  implicit none
+!======================================================================
+!#include "temps.h"
+#include "comcstfi.h"
+#include "planete.h"
+!======================================================================
+  character(len=*),intent(in) :: filename
+  integer,intent(in) :: nsoil
+  integer,intent(in) :: ngrid
+  integer,intent(in) :: nlay
+  integer,intent(in) :: nq
+  real,intent(in) :: phystep
+  real,intent(in) :: time
+  real,intent(in) :: tsurf(ngrid)
+  real,intent(in) :: tsoil(ngrid,nsoil)
+  real,intent(in) :: emis(ngrid)
+  real,intent(in) :: q2(ngrid,nlay+1)
+  real,intent(in) :: qsurf(ngrid,nq)
+  real,intent(in) :: cloudfrac(ngrid,nlay)
+  real,intent(in) :: totcloudfrac(ngrid)
+  real,intent(in) :: hice(ngrid)
+
+  integer :: iq
+  
+  ! Open file
+  call open_restartphy(filename)
+
+  ! First variable to write must be "Time", in order to correctly
+  ! set time counter in file
+  !call put_var("Time","Temps de simulation",time)
+  
+  ! Surface temperature
+  call put_field("tsurf","Surface temperature",tsurf)
+  
+  ! Soil temperature
+  call put_field("tsoil","Soil temperature",tsoil)
+  
+  ! Emissivity of the surface
+  call put_field("emis","Surface emissivity",emis)
+  
+  ! Planetary Boundary Layer
+  call put_field("q2","pbl wind variance",q2)
+  
+! cloud fraction and sea ice (NB: these should be optional... to be improved)
+  call put_field("cloudfrac","Cloud fraction",cloudfrac)
+  call put_field("totcloudfrac","Total fraction",totcloudfrac)
+  call put_field("hice","Height of oceanic ice",hice)
+
+! tracers
+  if (nq>0) then
+    do iq=1,nq
+      call put_field(tname(iq),"tracer on surface",qsurf(:,iq))
+    enddo
+  endif ! of if (nq>0)
+
+! close file
+      CALL close_restartphy
+!$OMP BARRIER
+
+end subroutine physdem1
+
+end module phyredem

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

@@ -10 +10 @@
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
       use watercommon_h, only : RLVTT, Psat_water,epsi
-      use gases_h
+      use gases_h, only: gnom, gfrac
       use radcommon_h, only: sigma, eclipse, glat, grav
       use radii_mod, only: h2o_reffrad, co2_reffrad, h2o_cloudrad
-      use aerosol_mod
-      use surfdat_h
-      use comdiurn_h
-      use comsaison_h
-      use comsoil_h
-      USE comgeomfi_h
-      USE tracer_h
+      use aerosol_mod, only: iaero_co2, iaero_h2o
+      use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe, &
+                           dryness, watercaptag
+      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
+      use comsaison_h, only: mu0, fract, dist_star, declin
+      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
+      USE comgeomfi_h, only: long, lati, area, totarea
+      USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
+                          alpha_lift, alpha_devil, qextrhor, &
+                          igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, &
+                          igcm_co2_ice
+      use control_mod, only: ecritphy, iphysiq, nday
+      use phyredem, only: physdem0, physdem1
 
       implicit none
@@ -80 +86 @@
 !    pt(ngrid,nlayer)      Temperature (K)
 !    pq(ngrid,nlayer,nq)   Advected fields
-!    pudyn(ngrid,nlayer)    \ 
+!    pudyn(ngrid,nlayer)    \
 !    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
 !    ptdyn(ngrid,nlayer)     / corresponding variables
@@ -123 +129 @@
 #include "comcstfi.h"
 #include "planete.h"
-#include "control.h"
+!#include "control.h"
 #include "netcdf.inc"
 
@@ -621 +627 @@
                             ! need epsi for the wvp definitions in callcorrk.F
 
+         if (ngrid.ne.1) then ! no need to create a restart file in 1d
+           call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
+                         ptimestep,pday+nday,time_phys,area, &
+                         albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+         endif
+         
       endif        !  (end of "if firstcall")
 
@@ -626 +638 @@
 ! 1.2   Initializations done at every physical timestep:
 ! ---------------------------------------------------
-
 !     Initialize various variables
 !     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -655 +666 @@
 !    Compute variations of g with latitude (oblate case)
 !
-        if (oblate .eq. .false.) then
+        if (oblate .eqv. .false.) then
            glat(:) = g
         else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
@@ -734 +745 @@
                call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
                ztim1,ztim2,ztim3,mu0,fract, flatten)
-           else if(diurnal .eq. .false.) then
+           else if(diurnal .eqv. .false.) then
 
                call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad,flatten)
@@ -745 +756 @@
            if(rings_shadow) then 
                 write(*,*) 'Rings shadow activated'
-                if(diurnal .eq. .false.) then ! we need to compute the daily average insolation 
+                if(diurnal .eqv. .false.) then ! we need to compute the daily average insolation 
                     pas = 1./nb_hours
                     ptime_day = 0.
@@ -784 +795 @@
             endif
 
-
            if (corrk) then
 
@@ -1531 +1541 @@
 
          if(meanOLR)then
-            if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
+            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
                ! to record global radiative balance
                open(92,file="rad_bal.out",form='formatted',position='append')
@@ -1607 +1617 @@
 
          if(meanOLR)then
-            if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
+            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
                ! to record global water balance
                open(98,file="h2o_bal.out",form='formatted',position='append')
@@ -1681 +1691 @@
             endif
 
-            write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
-#ifdef CPP_PARA
-! for now in parallel we use a different routine to write restartfi.nc
-            call phyredem(ngrid,"restartfi.nc",           &
-                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
-                    cloudfrac,totcloudfrac,hice)
-#else
-            call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
-                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
-                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
-                    cloudfrac,totcloudfrac,hice,noms)
-#endif
-         endif
+            if (ngrid.ne.1) then
+              write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
+!#ifdef CPP_PARA
+!! for now in parallel we use a different routine to write restartfi.nc
+!            call phyredem(ngrid,"restartfi.nc",           &
+!                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
+!                    cloudfrac,totcloudfrac,hice)
+!#else
+!            call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
+!                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
+!                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
+!                    cloudfrac,totcloudfrac,hice,noms)
+!#endif
+              call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
+                      ptimestep,ztime_fin, &
+                      tsurf,tsoil,emis,q2,qsurf_hist, &
+                      cloudfrac,totcloudfrac,hice)
+            endif
+         
+         endif ! of if (lastcall)
 
 

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

@@ -2 +2 @@
 
 ! to use  'getin'
-      use ioipsl_getincom 
-      use surfdat_h
-      use comdiurn_h
-      use comsaison_h
-      use comsoil_h
-      USE comgeomfi_h
-
+      use ioipsl_getincom, only: getin
+      use infotrac, only: nqtot, tname
+      use surfdat_h, only: albedodat, phisfi, dryness, watercaptag,
+     &                     zmea, zstd, zsig, zgam, zthe,
+     &                     emissiv, emisice, albedice, iceradius,
+     &                     dtemisice
+      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
+!      use comsaison_h
+      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
+      USE comgeomfi_h, only: lati, long, area
+      use control_mod, only: day_step, ecritphy
+      use phyredem, only: physdem0,physdem1
+      use comgeomphy, only: initcomgeomphy
       implicit none
 
@@ -41 +47 @@
 #include "comcstfi.h"
 #include "planete.h"
-#include "control.h"
+!#include "control.h"
 #include "comvert.h"
 #include "netcdf.inc"
@@ -61 +67 @@
       REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
       REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
-      REAL psurf,tsurf      
+      REAL psurf,tsurf(1)      
       REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
       REAL gru,grv   ! prescribed "geostrophic" background wind
       REAL temp(nlayermx)   ! temperature at the middle of the layers
-      REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
-      REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
+      REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
+      REAL,ALLOCATABLE :: qsurf(:)      ! tracer surface budget (e.g. kg.m-2)
       REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
 !      REAL co2ice           ! co2ice layer (kg.m-2) !not used anymore
@@ -72 +78 @@
       integer :: i_h2o_ice=0  ! tracer index of h2o ice
       integer :: i_h2o_vap=0  ! tracer index of h2o vapor
-      REAL emis             ! surface layer
+      REAL emis(1)             ! surface layer
       REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
       REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
@@ -80 +86 @@
       REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
       REAL dpsurf    
-      REAL dq(nlayermx,nqmx)
-      REAL dqdyn(nlayermx,nqmx)
+      REAL,ALLOCATABLE :: dq(:,:)
+      REAL,ALLOCATABLE :: dqdyn(:,:)
 
 c   Various intermediate variables
@@ -88 +94 @@
       REAL phi(nlayermx),h(nlayermx),s(nlayermx)
       REAL pks, ptif, w(nlayermx)
-      REAL qtotinit, mqtot(nqmx),qtot
       INTEGER ierr, aslun
       REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
@@ -112 +117 @@
 
 !     added by AS to avoid the use of adv trac common
-      character*20 nametrac(nqmx)   ! name of the tracer (no need for adv trac common)
-
+      character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
+
+      real :: latitude, longitude
 
 c=======================================================================
 c INITIALISATION
 c=======================================================================
+! initialize "serial/parallel" related stuff
+      CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
+      call initcomgeomphy
 
       !! those are defined in surfdat_h.F90
@@ -227 +236 @@
           ! read number of tracers:
           read(90,*,iostat=ierr) nq
+          nqtot=nq ! set value of nqtot (in infotrac module) as nq
           if (ierr.ne.0) then
             write(*,*) "rcm1d: error reading number of tracers"
             write(*,*) "   (first line of traceur.def) "
             stop
+          endif
+          if (nq>0) then
+            allocate(tname(nq))
+            allocate(q(nlayermx,nq))
+            allocate(qsurf(nq))
+            allocate(dq(nlayermx,nq))
+            allocate(dqdyn(nlayermx,nq))
           else
-            ! check that the number of tracers is indeed nqmx
-            if (nq.ne.nqmx) then
-              write(*,*) "rcm1d: error, wrong number of tracers:"
-              write(*,*) "nq=",nq," whereas nqmx=",nqmx
-              stop
-            endif
+            allocate(tname(1)) ! tname(1) is used below, even if nq=0
           endif
         
           do iq=1,nq
           ! minimal version, just read in the tracer names, 1 per line
-            read(90,*,iostat=ierr) nametrac(iq)
+            read(90,*,iostat=ierr) tname(iq)
             if (ierr.ne.0) then
               write(*,*) 'rcm1d: error reading tracer names...'
@@ -253 +265 @@
          i_h2o_vap=0
          do iq=1,nq
-           if (nametrac(iq)=="co2_ice") then
+           if (tname(iq)=="co2_ice") then
              i_co2_ice=iq
-           elseif (nametrac(iq)=="h2o_ice") then
+           elseif (tname(iq)=="h2o_ice") then
              i_h2o_ice=iq
-           elseif (nametrac(iq)=="h2o_vap") then
+           elseif (tname(iq)=="h2o_vap") then
              i_h2o_vap=iq
            endif
@@ -275 +287 @@
 
       else
-       nq=nqmx
+       nq=nqtot
+       if (nq>0) then
+            allocate(tname(nq))
+            allocate(q(nlayermx,nq))
+            allocate(qsurf(nq))
+            allocate(dq(nlayermx,nq))
+            allocate(dqdyn(nlayermx,nq))
+        else
+            allocate(tname(1)) ! tname(1) is used below, even if nq=0
+        endif
        ! Check that tracer boolean is compliant with number of tracers
        ! -- otherwise there is an error (and more generally we have to be consistent)
@@ -290 +311 @@
         do iq=1,nq
           write(str7,'(a1,i2.2)')'q',iq
-          nametrac(iq)=str7
+          tname(iq)=str7
         enddo
         ! actually, we'll need at least one "co2_ice" tracer
         ! (for surface CO2 ice)
-        nametrac(1)="co2_ice"
         i_co2_ice=1
+        tname(i_co2_ice)="co2_ice"
       endif ! of if (tracer)
 
@@ -315 +336 @@
       phisfi(1)=0.E+0
      !!! LATITUDE. can be set. 
-      lati(1)=0 ! default value for lati(1)
+      latitude=0 ! default value for latitude
       PRINT *,'latitude (in degrees) ?'
-      call getin("latitude",lati(1))
-      write(*,*) " latitude = ",lati(1)
-      lati(1)=lati(1)*pi/180.E+0
+      call getin("latitude",latitude)
+      write(*,*) " latitude = ",latitude
+      latitude=latitude*pi/180.E+0
      !!! LONGITUDE. useless in 1D.
-      long(1)=0.E+0
-      long(1)=long(1)*pi/180.E+0
+      longitude=0.E+0
+      longitude=longitude*pi/180.E+0
 
 !!! CALL INIFIS
@@ -332 +353 @@
 !!! - physical frequency: nevermind, in inifis this is a simple print
       cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
-      CALL inifis(1,llm,day0,daysec,dtphys,
-     .            lati,long,area,rad,g,r,cpp)
+      CALL inifis(1,llm,nq,day0,daysec,dtphys,
+     .            latitude,longitude,area,rad,g,r,cpp)
 !!! We check everything is OK.
       PRINT *,"CHECK"
@@ -558 +579 @@
 c  emissivity / surface co2 ice ( + h2o ice??)
 c  -------------------------------------------
-      emis=emissiv ! default value for emissivity
+      emis(1)=emissiv ! default value for emissivity
       PRINT *,'Emissivity of bare ground ?'
-      call getin("emis",emis)
-      write(*,*) " emis = ",emis
-      emissiv=emis ! we do this so that condense_co2 sets things to the right
+      call getin("emis",emis(1))
+      write(*,*) " emis = ",emis(1)
+      emissiv=emis(1) ! we do this so that condense_co2 sets things to the right
                    ! value if there is no snow
 
@@ -573 +594 @@
             ! if we have some CO2 ice on the surface, change emissivity
             if (lati(1).ge.0) then ! northern hemisphere
-              emis=emisice(1)
+              emis(1)=emisice(1)
             else ! southern hemisphere
-              emis=emisice(2)
+              emis(1)=emisice(2)
             endif
          ENDIF
@@ -673 +694 @@
       call profile(nlayer+1,tmp1,tmp2)
 
-      tsurf=tmp2(0)
+      tsurf(1)=tmp2(0)
       DO ilayer=1,nlayer
         temp(ilayer)=tmp2(ilayer)
       ENDDO
       print*,"check"
-      PRINT*,"INPUT SURFACE TEMPERATURE",tsurf
+      PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
       PRINT*,"INPUT TEMPERATURE PROFILE",temp
 
@@ -698 +719 @@
       DO isoil=1,nsoil
          inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
-         tsoil(isoil)=tsurf  ! soil temperature
+         tsoil(isoil)=tsurf(1)  ! soil temperature
       ENDDO
 
@@ -711 +732 @@
 
 
-c  Ecriture de "startfi"
+c  Write a "startfi" file
 c  --------------------
-c  (Ce fichier sera aussitot relu au premier
-c   appel de "physiq", mais il est necessaire pour passer
-c   les variables purement physiques a "physiq"...
-
-      call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
-     &              dtphys,float(day0),
-     &              time,tsurf,tsoil,emis,q2,qsurf,
-     &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
-     &              cloudfrac,totcloudfrac,hice,nametrac)
+c  This file will be read during the first call to "physiq".
+c  It is needed to transfert physics variables to "physiq"...
+
+      call physdem0("startfi.nc",long,lati,nsoilmx,1,llm,nq,
+     &              dtphys,real(day0),time,area,
+     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
+      call physdem1("startfi.nc",nsoilmx,1,llm,nq,
+     &                dtphys,time,
+     &                tsurf,tsoil,emis,q2,qsurf,
+     &                cloudfrac,totcloudfrac,hice)
+
+!      call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
+!     &              dtphys,float(day0),
+!     &              time,tsurf,tsoil,emis,q2,qsurf,
+!     &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
+!     &              cloudfrac,totcloudfrac,hice,nametrac)
 
 c=======================================================================
@@ -775 +803 @@
 
       CALL physiq (1,llm,nq,
-     .     nametrac,
+     .     tname,
      ,     firstcall,lastcall,
      ,     day,time,dtphys,
@@ -858 +886 @@
       ENDDO   ! fin de la boucle temporelle
 
+      write(*,*) "rcm1d: Everything is cool."
+
 c    ========================================================
       end                       !rcm1d
@@ -874 +904 @@
 
 #include "../dyn3d/disvert.F"
+#include "../dyn3d/abort_gcm.F"

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

@@ -4 +4 @@
      &          capcal,fluxgrd)
 
-      use comsoil_h
+      use comsoil_h, only: layer, mlayer, volcapa
 
       implicit none
@@ -17 +17 @@
 !-----------------------------------------------------------------------
 
-#include "dimensions.h"
-#include "dimphys.h"
+      include "dimensions.h"
+      include "dimphys.h"
 
 
@@ -25 +25 @@
 !  ---------
 !  inputs:
-      integer ngrid     ! number of (horizontal) grid-points 
-      integer nsoil     ! number of soil layers 
-      logical firstcall ! identifier for initialization call 
-      logical lastcall
-      real therm_i(ngrid,nsoil) ! thermal inertia
-      real timestep         ! time step
-      real tsurf(ngrid)   ! surface temperature
+      integer,intent(in) :: ngrid       ! number of (horizontal) grid-points 
+      integer,intent(in) :: nsoil       ! number of soil layers 
+      logical,intent(in) :: firstcall ! identifier for initialization call 
+      logical,intent(in) :: lastcall
+      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
+      real,intent(in) :: timestep           ! time step
+      real,intent(in) :: tsurf(ngrid)   ! surface temperature
 ! outputs:
-      real tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
-      real capcal(ngrid) ! surface specific heat
-      real fluxgrd(ngrid) ! surface diffusive heat flux
+      real,intent(out) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
+      real,intent(out) :: capcal(ngrid) ! surface specific heat
+      real,intent(out) :: fluxgrd(ngrid) ! surface diffusive heat flux
 
 ! local saved variables:
@@ -50 +50 @@
 
 
-!      print*,'tsoil=',tsoil
-!      print*,'tsurf=',tsurf
-
 ! 0. Initialisations and preprocessing step
       if (firstcall) then
@@ -58 +55 @@
       !       and not changed by soil.F 
 
-      ALLOCATE(mthermdiff(ngrid,0:nsoilmx-1)) ! mid-layer thermal diffusivity
-      ALLOCATE(thermdiff(ngrid,nsoilmx-1))    ! inter-layer thermal diffusivity
-      ALLOCATE(coefq(0:nsoilmx-1))              ! q_{k+1/2} coefficients
-      ALLOCATE(coefd(ngrid,nsoilmx-1))        ! d_k coefficients
-      ALLOCATE(alph(ngrid,nsoilmx-1))         ! alpha_k coefficients
-      ALLOCATE(beta(ngrid,nsoilmx-1))         ! beta_k coefficients
+      ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity
+      ALLOCATE(thermdiff(ngrid,nsoil-1))    ! inter-layer thermal diffusivity
+      ALLOCATE(coefq(0:nsoil-1))              ! q_{k+1/2} coefficients
+      ALLOCATE(coefd(ngrid,nsoil-1))        ! d_k coefficients
+      ALLOCATE(alph(ngrid,nsoil-1))         ! alpha_k coefficients
+      ALLOCATE(beta(ngrid,nsoil-1))         ! beta_k coefficients
 
 ! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
@@ -187 +184 @@
       enddo
 
-      if (lastcall) then
-        IF ( ALLOCATED(mthermdiff)) DEALLOCATE(mthermdiff)
-        IF ( ALLOCATED(thermdiff)) DEALLOCATE(thermdiff)
-        IF ( ALLOCATED(coefq)) DEALLOCATE(coefq)
-        IF ( ALLOCATED(coefd)) DEALLOCATE(coefd)
-        IF ( ALLOCATED(alph)) DEALLOCATE(alph)
-        IF ( ALLOCATED(beta)) DEALLOCATE(beta)
-      endif
-
       end
 

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

@@ -1 @@
-      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
-
-      use comsoil_h
-
+      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime)
+
+!      use netcdf
+      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
+      use iostart, only: inquire_field_ndims, get_var, get_field,
+     &                   inquire_field, inquire_dimension_length
       implicit none
 
@@ -9 +11 @@
 !
 !  Purpose: Read and/or initialise soil depths and properties
+!
+! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using
+!                      r4 or r8 restarts independently of having compiled
+!                      the GCM in r4 or r8)
+!                June 2013 TN : Possibility to read files with a time axis
+!
 !
 !  This subroutine reads from a NetCDF file (opened by the caller)
@@ -26 +34 @@
 !======================================================================
 
-#include "dimensions.h"
-#include "dimphys.h"
-
-#include "netcdf.inc"
 !======================================================================
 !  arguments
 !  ---------
 !  inputs:
-      integer nid       ! Input Netcdf file ID 
-      integer ngrid     ! # of horizontal grid points
-      integer nsoil     ! # of soil layers
-      integer sta       ! # at which reading starts
-      real tsurf(ngrid) ! surface temperature
+      integer,intent(in) :: nid ! Input Netcdf file ID 
+      integer,intent(in) :: ngrid       ! # of horizontal grid points
+      integer,intent(in) :: nsoil       ! # of soil layers
+      real,intent(in) :: tsurf(ngrid)   ! surface temperature
+      integer,intent(in) :: indextime   ! position on time axis
 !  output:
-      real tsoil(ngrid,nsoilmx) ! soil temperature
+      real,intent(out) :: tsoil(ngrid,nsoil)    ! soil temperature
 
 !======================================================================
@@ -50 +54 @@
       integer ndims     ! # of dimensions of read <inertiedat> data
       integer ig,iloop  ! loop counters
+      
+      integer edges(3),corner(3) ! to read a specific time
 
       logical :: olddepthdef=.false. ! flag
@@ -70 +76 @@
       real,parameter :: default_volcapa=1.e6
 
-      integer, dimension(2) :: sta2d
-      integer, dimension(2) :: siz2d
-
-!======================================================================
-
-      !! this is defined in comsoil_h
-      IF (.not.ALLOCATED(layer)) 
-     . ALLOCATE(layer(nsoilmx))
-      IF (.not.ALLOCATED(mlayer)) 
-     . ALLOCATE(mlayer(0:nsoilmx-1))
-      IF (.not.ALLOCATED(inertiedat)) 
-     . ALLOCATE(inertiedat(ngrid,nsoilmx))
-
-      !! this is defined in dimphys.h
-      sta = cursor
+      logical :: found,ok
+      
+!======================================================================
 
 ! 1. Depth coordinate
@@ -91 +85 @@
 ! 1.1 Start by reading how many layers of soil there are
 
-        ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
-        if (ierr.ne.NF_NOERR) then
-         write(*,*)'soil_settings: Problem reading <subsurface_layers>'
-         call abort
-        endif
-
-        ierr=NF_INQ_DIMLEN(nid,dimid,dimlen)
-        if (ierr.ne.NF_NOERR) then
-         write(*,*)'soil_settings: Problem getting <subsurface_layers>',
-     &             'length'
-         call abort
-        endif
+        dimlen=inquire_dimension_length("subsurface_layers")
 
         if (dimlen.ne.nsoil) then
@@ -112 +95 @@
         endif
 
-      sta2d = (/sta,1/)
-      siz2d = (/ngrid,dimlen/)
-
-
 ! 1.2 Find out the # of dimensions <inertiedat> was defined as using
-!     (in ye old days, thermal inertia was only given at the "surface")
-      ! Look for thermal inertia data
-      ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)
-      if (ierr.NE.NF_NOERR) then
-         write(*,*)'soil_settings: Field <inertiedat> not found!'
-         call abort
-      endif
-
-      ! Read the # of dimensions <inertidat> was defined as using
-      ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims)
-      ! if (ndims.eq.1) then we have the "old 2D-surface" format
+
+      ndims=inquire_field_ndims("inertiedat")
 
 ! 1.3 Read depths values or set olddepthdef flag and values
@@ -147 +117 @@
         enddo
       else ! Look for depth 
-        ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
-        if (ierr.ne.NF_NOERR) then
-          write(*,*)'soil_settings: Field <soildepth> not found!'
-          call abort
-        endif
         ! read <depth> coordinate
         if (interpol) then !put values in oldmlayer
-          if (.not.allocated(oldmlayer)) then
-             allocate(oldmlayer(dimlen),stat=ierr)
-          endif
-#ifdef NC_DOUBLE
-           ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldmlayer)
-#else
-           ierr = NF_GET_VAR_REAL(nid,nvarid,oldmlayer)
-#endif
-          if (ierr.ne.NF_NOERR) then
-           write(*,*)'soil_settings: Problem while reading <soildepth>'
-           call abort
+          call get_var("soildepth",oldmlayer,found)
+          if (.not.found) then
+            write(*,*)'soil_settings: Problem while reading <soildepth>'
           endif
         else ! put values in mlayer
-#ifdef NC_DOUBLE
-           ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
-#else
-           ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
-#endif
-          if (ierr.ne.NF_NOERR) then
-           write(*,*)'soil_settings: Problem while reading <soildepth>'
-           call abort
+          call get_var("soildepth",mlayer,found)
+          if (.not.found) then
+            write(*,*)'soil_settings: Problem while reading <soildepth>'
           endif
         endif !of if (interpol)
@@ -183 +135 @@
       ! default mlayer distribution, following a power law:
       !  mlayer(k)=lay1*alpha**(k-1/2)
-        lay1=2.e-4  !mars case (nsoilmx=18)
-!        lay1=3.e-2  !earth case (nsoilmx=13)
+        lay1=2.e-4
         alpha=2
         do iloop=0,nsoil-1
@@ -200 +151 @@
       enddo
 
-! 2. Volumetric heat capacity (note: it is declared in comsoil.h)
+! 2. Volumetric heat capacity (note: it is declared in comsoil_h)
 ! ---------------------------
 ! "volcapa" is (so far) 0D and written in "controle" table of startfi file 
@@ -214 +165 @@
         volcapa=default_volcapa
       endif
-! Look for it
-!      ierr=NF_INQ_VARID(nid,"volcapa",nvarid)
-!      if (ierr.NE.NF_NOERR) then
-!         write(*,*)'soil_settings: Field <volcapa> not found!'
-!         write(*,*)'Initializing Volumetric heat capacity to ',
-!     &             default_volcapa
-!         volcapa=default_volcapa
-!      else
-!#ifdef NC_DOUBLE
-!       ierr = NF_GET_VAR_DOUBLE(nid,nvarid,volcapa)
-!#else
-!       ierr = NF_GET_VAR_REAL(nid,nvarid,volcapa)
-!#endif
-!        if (ierr.ne.NF_NOERR) then
-!         write(*,*)'soil_settings: Problem while reading <volcapa>'
-!         call abort
-!       endif
-!      endif
-
-! 3. Thermal inertia (note: it is declared in comsoil.h)
+
+! 3. Thermal inertia (note: it is declared in comsoil_h)
 ! ------------------
 
 ! 3.1 Look (again) for thermal inertia data (to reset nvarid)
-      ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)
-      if (ierr.NE.NF_NOERR) then
-         write(*,*)'soil_settings: Field <inertiedat> not found!'
-         call abort
-      endif
 
 ! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
@@ -251 +179 @@
        ! Read Surface thermal inertia
        allocate(surfinertia(ngrid))
-#ifdef NC_DOUBLE
-       ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, surfinertia)
-#else
-       ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, surfinertia)
-#endif
-        if (ierr.NE.NF_NOERR) then
-         write(*,*)'soil_settings: Problem while reading <inertiedat>'
+       call get_field("inertiedat",surfinertia,found)
+       if (.not.found) then
+         write(*,*) "soil_settings: Failed loading <inertiedat>"
          call abort
-        endif
+       endif
        
        write(*,*)' => Building soil thermal inertia (using reference sur
@@ -277 +201 @@
             stop
            endif
-         endif
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldinertiedat)
-#else
-        ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldinertiedat)
-#endif
-        if (ierr.NE.NF_NOERR) then
-         write(*,*)'soil_settings: Problem while reading <inertiedat>'
-         call abort
+         endif ! of if (.not.allocated(oldinertiedat))
+        call get_field("inertiedat",oldinertiedat,found)
+        if (.not.found) then
+          write(*,*) "soil_settings: Failed loading <inertiedat>"
+          call abort
         endif
        else ! put values in therm_i
-#ifdef NC_DOUBLE
-        ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,inertiedat)
-#else
-        ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,inertiedat)
-#endif
-        if (ierr.NE.NF_NOERR) then
-         write(*,*)'soil_settings: Problem while reading <inertiedat>'
-         call abort
-        endif
+         call get_field("inertiedat",inertiedat,found)
+         if (.not.found) then
+           write(*,*) "soil_settings: Failed loading <inertiedat>"
+           call abort
+         endif
+!        endif
        endif ! of if (interpol)
       endif ! of if (ndims.eq.1)
@@ -303 +220 @@
 ! -------------------------
 
-      ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
-      if (ierr.ne.NF_NOERR) then
+!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
+      ok=inquire_field("tsoil")
+!      if (ierr.ne.nf90_noerr) then
+      if (.not.ok) then
         write(*,*)'soil_settings: Field <tsoil> not found!'
         write(*,*)' => Building <tsoil> from surface values <tsurf>'
@@ -320 +239 @@
            endif
          endif
-#ifdef NC_DOUBLE
-        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldtsoil)
-#else
-        ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldtsoil)
-#endif
-        if (ierr.ne.NF_NOERR) then
-         write(*,*)'soil_settings: Problem while reading <tsoil>'
-         call abort
-        endif
+         call get_field("tsoil",oldtsoil,found)
+         if (.not.found) then
+           write(*,*) "soil_settings: Failed loading <tsoil>"
+           call abort
+         endif
        else ! put values in tsoil
-#ifdef NC_DOUBLE
-        ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,tsoil)
-#else
-        ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,tsoil)
-#endif
-        if (ierr.ne.NF_NOERR) then
-         write(*,*)'soil_settings: Problem while reading <tsoil>'
-         call abort
-        endif
+         call get_field("tsoil",tsoil,found,timeindex=indextime)
+         if (.not.found) then
+           write(*,*) "soil_settings: Failed loading <tsoil>"
+           call abort
+         endif
        endif ! of if (interpol) 
-      endif
+      endif! of if (.not.ok)
 
 ! 5. If necessary, interpolate soil temperatures and thermal inertias

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

@@ -19 +19 @@
 c=======================================================================
 
+      use infotrac, only: iniadvtrac, nqtot, tname
       USE comsoil_h
-      USE comgeomfi_h
-
+      USE comgeomfi_h, ONLY: lati, long, area
+!      use control_mod
+      use comgeomphy, only: initcomgeomphy
       implicit none
 
@@ -32 +34 @@
 #include "logic.h"
 #include "temps.h"
-#include "control.h"
+!#include "control.h"
 #include "ener.h"
 
 #include "dimphys.h"
 #include "planete.h"
-#include"advtrac.h"
+!#include"advtrac.h"
 #include "netcdf.inc"
 
@@ -48 +50 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL teta(ip1jmp1,llm)                    ! temperature potentielle 
-      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
+      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
       REAL pks(ip1jmp1)                      ! exner (f pour filtre)
       REAL pk(ip1jmp1,llm)
@@ -63 +65 @@
       REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
       REAL co2ice(ngridmx)      ! CO2 ice layer
-      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
+      REAL q2(ngridmx,nlayermx+1)
+      REAL,ALLOCATABLE :: qsurf(:,:)
       REAL emis(ngridmx)
       INTEGER start,length
@@ -83 +86 @@
       REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
       REAL co2iceS(ip1jmp1)
-      REAL q2S(ip1jmp1,llm+1),qsurfS(ip1jmp1,nqmx)
+      REAL q2S(ip1jmp1,llm+1)
+      REAL,ALLOCATABLE :: qsurfS(:,:)
       REAL emisS(ip1jmp1)
 
@@ -122 +126 @@
 c-----------------------------------------------------------------------
 
+      CALL defrun_new(99, .TRUE. )
       grireg   = .TRUE.
       cursor = 1 ! added by AS in dimphys.
+
+! initialize "serial/parallel" related stuff
+      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
+      call initcomgeomphy
 
       ! ALLOCATE ARRAYS IN comgeomfi_h (usually done in inifis)
@@ -134 +143 @@
 c Lecture des donnees
 c=======================================================================
-! Load tracer names:
-      call iniadvtrac(nq,numvanle)
-! Initialize global tracer indexes (stored in tracer.h)
-! ... this has to be done before phyetat0
-      call initracer(ngridmx,nqmx,tnom)
+! Load tracer number and names:
+      call iniadvtrac(nqtot,numvanle)
+
+! allocate arrays:
+      allocate(q(ip1jmp1,llm,nqtot))
+      allocate(qsurf(ngridmx,nqtot))
+      allocate(qsurfS(ip1jmp1,nqtot))
 
       fichnom = 'start.nc'
-      CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
+      CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
      .       ps,phis,timedyn)
 
@@ -173 +184 @@
       Lmodif=0
 
-      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,
+      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqtot,day_ini_fi,
      .      timefi,
      .      tsurf,tsoil,emis,q2,qsurf,
@@ -286 +297 @@
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
       call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
-      call gr_fi_dyn(nqmx,ngridmx,iip1,jjp1,qsurf,qsurfS)
+      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
       call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
@@ -390 +401 @@
 
 c-----------------------------------------------------------------------
-c Ecriture du champs  q  ( q[1,nqmx] )
-c-----------------------------------------------------------------------
-      do iq=1,nqmx
-        call write_archive(nid,ntime,tnom(iq),'tracer','kg/kg',
+c Ecriture du champs  q  ( q[1,nqtot] )
+c-----------------------------------------------------------------------
+      do iq=1,nqtot
+        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
      &         3,q(1,1,iq))
       end do
 c-----------------------------------------------------------------------
-c Ecriture du champs  qsurf  ( qsurf[1,nqmx] )
-c-----------------------------------------------------------------------
-      do iq=1,nqmx
-        txt=trim(tnom(iq))//"_surf"
+c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
+c-----------------------------------------------------------------------
+      do iq=1,nqtot
+        txt=trim(tname(iq))//"_surf"
         call write_archive(nid,ntime,txt,'Tracer on surface',
      &  'kg.m-2',2,qsurfS(1,iq))
@@ -448 +459 @@
       ierr=NF_CLOSE(nid)
 
+      write(*,*) "start2archive: All is well that ends well."
+
       end 

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

@@ -65 +65 @@
 C             ---------
 C
-      INTEGER kgrid
-      REAL ptim1,ptim2,ptim3, pflat
-      REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)
-      REAL psilat(kgrid), pcolat(kgrid)
+      INTEGER,INTENT(IN) :: kgrid
+      REAL,INTENT(IN) :: ptim1,ptim2,ptim3, pflat
+      REAL,INTENT(IN) :: psilon(kgrid),pcolon(kgrid)
+      REAL,INTENT(IN) :: psilat(kgrid), pcolat(kgrid)
+      REAL,INTENT(OUT) :: pmu0(kgrid),pfract(kgrid)
 C
       INTEGER jl
@@ -100 +101 @@
         ztim3=pcolat(jl)*ptim3
         pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)
-        pmu0(jl)=pmu0(jl)/SQRT(pcolat(jl)**2 + (rap**2) * (psilat(jl)**2))
+        pmu0(jl)=pmu0(jl)/SQRT(pcolat(jl)**2+(rap**2)*(psilat(jl)**2))
 
       ENDDO

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

@@ -1 +1 @@
       SUBROUTINE surfini(ngrid,nq,qsurf,psolaralb)
 
-      USE surfdat_h
-      USE tracer_h
+      USE surfdat_h, only: albedodat, albedice
+      USE tracer_h, only: igcm_co2_ice
+      use comgeomfi_h, only: lati
+      use planetwide_mod, only: planetwide_maxval, planetwide_minval
 
       IMPLICIT NONE
@@ -17 +19 @@
 #include "callkeys.h"
 c
-      INTEGER ngrid,ig,icap
-      INTEGER nq
-      REAL  piceco2(ngrid),psolaralb(ngrid)
-      REAL qsurf(ngrid,nq) !tracer on surface (kg/m2)
+      INTEGER,INTENT(IN) :: ngrid
+      INTEGER,INTENT(IN) :: nq
+      REAL,INTENT(OUT) :: psolaralb(ngrid)
+      REAL,INTENT(IN) :: qsurf(ngrid,nq) !tracer on surface (kg/m2)
 
-      EXTERNAL ISMIN,ISMAX
-      INTEGER ISMIN,ISMAX
+      INTEGER :: ig,icap
+      REAL :: min_albedo,max_albedo
 c
 c=======================================================================
 
-c
-c     calcul de piceco2 (kg/m2) a l'etat initial
-c     ------------------------------------------
 
       DO ig=1,ngrid
          psolaralb(ig)=albedodat(ig)
-!         psolaralb(ig,2)=albedodat(ig)
       ENDDO
 
-      PRINT*,'surfini: minimum des donnees albedo',
-     s     albedodat(ISMIN(ngrid,albedodat,1))
-      PRINT*,'surfini: maximum des donnees albedo',
-     s     albedodat(ISMAX(ngrid,albedodat,1))
-
-c      calcul de psolaralb 
-c      -------------------
-!      DO ig=1,ngrid
-c        IF (water) THEN 
-c          if (qsurf(ig,nq).gt.0.005) then
-c             psolaralb(ig,1) = 0.4
-c             psolaralb(ig,2) = 0.4 
-c           endif 
-c         ENDIF
-!      ENDDO ! of DO ig=1,ngrid     
-c IF there is more than 5 pr. um of h2o ice but no C02 ice, surface albedo is set to 0.4. 
+      call planetwide_minval(albedodat,min_albedo)
+      call planetwide_maxval(albedodat,max_albedo)
+      write(*,*) 'surfini: minimum corrected albedo',min_albedo
+      write(*,*) 'surfini: maximum corrected albedo',max_albedo
 
       if (igcm_co2_ice.ne.0) then
@@ -57 +43 @@
         DO ig=1,ngrid
           IF (qsurf(ig,igcm_co2_ice) .GT. 0.) THEN
-             !!!! no good in parallel
-             IF(ig.GT.ngrid/2+1) THEN
-                icap=2
+             IF(lati(ig).LT.0.) THEN
+                icap=2 ! Southern hemisphere
              ELSE
-                icap=1
+                icap=1 ! Northern hemisphere
              ENDIF
              psolaralb(ig) = albedice(icap)
-!             psolaralb(ig,2) =  albedice(icap)
           END IF
         ENDDO ! of DO ig=1,ngrid     
@@ -72 +56 @@
       endif
 
-      PRINT*,'surfini: minimum des donnees albedo',
-     s     psolaralb(ISMIN(ngrid,psolaralb,1))
-      PRINT*,'surfini: maximum des donnees albedo',
-     s     psolaralb(ISMAX(ngrid,psolaralb,1))
+      call planetwide_minval(psolaralb,min_albedo)
+      call planetwide_maxval(psolaralb,max_albedo)
+      write(*,*) 'surfini: minimum corrected albedo',min_albedo
+      write(*,*) 'surfini: maximum corrected albedo',max_albedo
 
-      RETURN
       END

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

@@ -44 +44 @@
 c=======================================================================
 ! to use  'getin'
-      use ioipsl_getincom 
-
-      USE surfdat_h
-      use comsoil_h
-      USE comgeomfi_h
+      use ioipsl_getincom , only: getin
+
+      use surfdat_h, only: albedice, emisice, iceradius, dtemisice,
+     &                     emissiv
+      use comsoil_h, only: volcapa
+      use iostart, only: get_var
+      use mod_phys_lmdz_para, only: is_parallel
 
       implicit none
  
-#include "dimensions.h"
-#include "dimphys.h"
 #include "comcstfi.h"
 #include "planete.h"
@@ -63 +63 @@
 c-----------------------------------------------------------------------
 
-      INTEGER ngrid
-
 c Arguments
 c ---------
-      INTEGER nid,nvarid,tab0
-      INTEGER*4 day_ini
-      INTEGER Lmodif
-      INTEGER lmax
-      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
+      INTEGER,INTENT(IN) :: ngrid,nid,tab0
+      INTEGER*4,INTENT(OUT) :: day_ini
+      INTEGER,INTENT(IN) :: Lmodif
+      INTEGER,INTENT(OUT) :: lmax
+      REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
 
 c Variables
 c ---------
-      INTEGER length
-      parameter (length = 100)
+      INTEGER,PARAMETER :: length=100
       REAL tab_cntrl(length) ! array in which are stored the run's parameters
-      INTEGER  ierr
+      INTEGER  ierr,nvarid
       INTEGER size
       CHARACTER modif*20
-
+      LOGICAL :: found
+      
+      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
+
+      IF (nid.eq.0) then
+c-----------------------------------------------------------------------
+c  Initialization of various physical constants to defaut values (nid = 0 case)
+c-----------------------------------------------------------------------
+      ELSE
 c-----------------------------------------------------------------------
 c  Initialization of physical constants by reading array tab_cntrl(:)
@@ -88 +93 @@
 c Read 'controle' array
 c
-      ierr = NF_INQ_VARID (nid, "controle", nvarid)
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Tabfi: Could not find <controle> data"
-         CALL abort
-      ENDIF
-#ifdef NC_DOUBLE
-      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
-#else
-      ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
-#endif
-      IF (ierr .NE. NF_NOERR) THEN
-         PRINT*, "Tabfi: Failed reading <controle> array"
-         CALL abort
-      ENDIF
-
-       print*,'tabfi: tab_cntrl',tab_cntrl
+
+       call get_var("controle",tab_cntrl,found)
+       if (.not.found) then
+         write(*,*)"tabfi: Failed reading <controle> array"
+         call abort
+       else
+         write(*,*)'tabfi: tab_cntrl',tab_cntrl
+       endif
 c
 c  Initialization of some physical constants
 c informations on physics grid
-      if(ngrid.ne.tab_cntrl(tab0+1)) then
-         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
-         print*,tab_cntrl(tab0+1),ngrid
-      endif
+!      if(ngrid.ne.tab_cntrl(tab0+1)) then
+!         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
+!         print*,tab_cntrl(tab0+1),ngrid
+!      endif
       lmax = nint(tab_cntrl(tab0+2))
       day_ini = tab_cntrl(tab0+3)
@@ -146 +143 @@
 c soil properties
       volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
-
-
-
-
 c-----------------------------------------------------------------------
 c       Save some constants for later use (as routine arguments)
@@ -160 +153 @@
       p_rad=rad
 
+      ENDIF    ! end of (nid = 0) 
 
 c-----------------------------------------------------------------------
@@ -209 +203 @@
 c-----------------------------------------------------------------------
 c        Modifications...
+! NB: Modifying controls should only be done by newstart, and in seq mode
+      if ((Lmodif.eq.1).and.is_parallel) then
+        write(*,*) "tabfi: Error modifying tab_control should",
+     &             " only happen in serial mode (eg: by newstart)"
+        stop
+      endif
 c-----------------------------------------------------------------------
 

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

@@ -32 +32 @@
 c=======================================================================
 
+      use comsoil_h, only: nsoilmx
       implicit none
 
 #include "dimensions.h"
-#include "dimphys.h"
 #include "paramet.h"
-#include "control.h"
 #include "comvert.h"
 #include "comgeom.h"

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

@@ -39 +39 @@
 !
 !=================================================================
-
-      USE surfdat_h
-#ifdef CPP_PARA
+      use surfdat_h, only: phisfi
+      use control_mod, only: ecritphy, day_step, iphysiq
       USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
      &                               is_master, gather
       USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
-#endif
       implicit none
 
 ! Commons
-#include "dimensions.h"
-#include "dimphys.h"
-#include "paramet.h"
-#include "control.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "netcdf.inc"
-#include "temps.h"
+      include "dimensions.h"
+      include "paramet.h"
+      include "comvert.h"
+      include "comgeom.h"
+      include "netcdf.inc"
+      include "temps.h"
 
 ! Arguments on input:
@@ -73 +69 @@
       real*4,save :: date
 
-      REAL phis(ip1jmp1) ! surface geopotential on extended lonxlat grid
+      REAL phis(ip1jmp1)
 
       integer irythme
@@ -96 +92 @@
       integer,save :: n_nom_def
       integer :: n
-      integer,parameter :: n_nom_def_max=99
-      character(len=20),save :: nom_def(n_nom_def_max)
+      integer,parameter :: n_nom_def_max=199
+      character(len=120),save :: nom_def(n_nom_def_max)
       logical,save :: firstcall=.true.
-
-      integer dimvert
+      
+#ifndef MESOSCALE
 
 #ifdef CPP_PARA
@@ -113 +109 @@
       real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
 #else
-      logical,parameter :: is_parallel=.false.
-      logical,parameter :: is_master=.true.
-      logical,parameter :: is_mpi_root=.true.
       real phisfi_glo(ngrid) ! surface geopotential on global physics grid
 #endif
+
 !***************************************************************
 !Sortie des variables au rythme voulu
 
-      irythme = int(ecritphy) ! sortie au rythme de ecritphy
+      irythme = ecritphy ! sortie au rythme de ecritphy
 !     irythme = iconser  ! sortie au rythme des variables de controle
 !     irythme = iphysiq  ! sortie a tous les pas physique
@@ -129 +123 @@
 !***************************************************************
 
-
 ! At very first call, check if there is a "diagfi.def" to use and read it
 ! -----------------------------------------------------------------------
@@ -135 +128 @@
          firstcall=.false.
 
-!      Open diagfi.def definition file if there is one:
+  !      Open diagfi.def definition file if there is one:
          open(99,file="diagfi.def",status='old',form='formatted',
      s   iostat=ierr2)
@@ -170 +163 @@
       end if
 
-
 ! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
 ! ------------------------------------------------------------------------
@@ -185 +177 @@
            stop
          endif
-
+         
          zitau = -1 ! initialize zitau
 
@@ -240 +232 @@
         ! are undefined; so set them to 1
         iphysiq=1 
-        !JL11 irythme=1
+        irythme=1
         ! NB: 
       endif
@@ -546 +538 @@
       endif
 
+#endif
       end
-

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

@@ -10 +10 @@
 ! (yielding the sought time series of the variable)
 
-#ifdef CPP_PARA
+! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
+
+use comsoil_h, only: nsoilmx
+use control_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
-#endif
 
 implicit none
 
 #include"dimensions.h"
-#include"dimphys.h"
 #include"paramet.h"
-#include"control.h"
 #include"netcdf.inc"
 
@@ -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 common parameter set in 'dimphys.h'
+! Note: nsoilmx is a parameter set in 'comsoil_h'
 
 ! Local variables:
-real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data on extended lonxlat grid
-! Note iip1,jjp1 known from paramet.h; nsoilmx known from dimphys.h
-real,dimension(iip1,jjp1) :: data2 ! to store 2D data on extended lonxlat grid
-real :: data0 ! to store 0D data
+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 :: data0 ! to store 0D data
 integer :: i,j,l ! for loops
 integer :: ig0
 
-real,save :: date ! time counter (in elapsed days)
+real*4,save :: date ! time counter (in elapsed days)
 integer,save :: isample ! sample rate at which data is to be written to output
 integer,save :: ntime=0 ! counter to internally store time steps
@@ -63 +63 @@
 real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
 real px2(ngrid)
-#else
-logical,parameter :: is_master=.true.
-logical,parameter :: is_mpi_root=.true.
 #endif
 
@@ -83 +80 @@
   
   ! Set output sample rate
-  isample=int(ecritphy) ! same as for diagfi outputs
+  isample=ecritphy ! same as for diagfi outputs
   ! Note ecritphy is known from control.h
   
@@ -128 +125 @@
      ierr=NF_INQ_VARID(nid,"time",varid)
      ! Add the current value of date to the "time" array
-#ifdef NC_DOUBLE
-     ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
-#else
+!#ifdef NC_DOUBLE
+!     ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
+!#else
      ierr=NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
-#endif
+!#endif
      if (ierr.ne.NF_NOERR) then
       write(*,*)"writediagsoil: Failed writing date to time variable"
@@ -203 +200 @@
   
   ! B.3. Write the slab of data
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
-#else
+!#ifdef NC_DOUBLE
+!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
+!#else
   ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3)
-#endif
+!#endif
   if (ierr.ne.NF_NOERR) then
     write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
@@ -273 +270 @@
   
   ! B.3. Write the slab of data
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
-#else
+!#ifdef NC_DOUBLE
+!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
+!#else
   ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2)
-#endif
+!#endif
   if (ierr.ne.NF_NOERR) then
     write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
@@ -312 +309 @@
 
   ! B.3. Write the data
-#ifdef NC_DOUBLE
-  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0)
-#else
+!#ifdef NC_DOUBLE
+!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0)
+!#else
   ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data0)
-#endif
+!#endif
   if (ierr.ne.NF_NOERR) then
     write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&

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

@@ -48 +48 @@
       use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
 #endif
+      use control_mod, only: ecritphy, iphysiq, day_step
 
       implicit none
@@ -55 +56 @@
 #include "dimphys.h"
 #include "paramet.h"
-#include "control.h"
+!#include "control.h"
 #include "comvert.h"
 #include "comgeom.h"
@@ -114 +115 @@
 !Sortie des variables au rythme voulu
 
-      irythme = int(ecritphy)*iradia ! sortie au rythme de ecritphy*iradia
+      irythme = ecritphy*iradia ! sortie au rythme de ecritphy*iradia
 !EM+JL if the spetra need to be output more frequently, need to define a ecritSpec...
 !     irythme = iphysiq  ! sortie a tous les pas physique

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

@@ -48 +48 @@
       use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
 #endif
+      use control_mod, only: ecritphy, iphysiq, day_step
 
       implicit none
@@ -55 +56 @@
 #include "dimphys.h"
 #include "paramet.h"
-#include "control.h"
+!#include "control.h"
 #include "comvert.h"
 #include "comgeom.h"
@@ -114 +115 @@
 !Sortie des variables au rythme voulu
 
-      irythme = int(ecritphy)*iradia ! sortie au rythme de ecritphy
+      irythme = ecritphy*iradia ! sortie au rythme de ecritphy
 !EM+JL if the spetra need to be output more frequently, need to define a ecritSpec...
 !     irythme = iphysiq  ! sortie a tous les pas physique

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

@@ -1 +1 @@
 subroutine wstats(ngrid,nom,titre,unite,dim,px)
 
-#ifdef CPP_PARA
 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
-#endif
 
 implicit none
@@ -47 +45 @@
 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
-logical,parameter :: is_master=.true.
-logical,parameter :: is_mpi_root=.true.
-integer,parameter :: klon_mpi_begin=1
 #endif
 
@@ -290 +285 @@
 endif ! of if (is_master)
 
-end
+end subroutine wstats
 
 !======================================================
@@ -367 +362 @@
 endif
 
-end
+end subroutine inivar
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
+
+! This subroutine defines variable 'name' in a (pre-existing and opened)
+! NetCDF file (known from its NetCDF ID 'nid').
+! The number of dimensions 'nbdim' of the variable, as well as the IDs of
+! corresponding dimensions must be set (in array 'dimids').
+! Upon successfull definition of the variable, 'nvarid' contains the
+! NetCDF ID of the variable.
+! The variables' attributes 'title' (Note that 'long_name' would be more
+! appropriate) and 'units' are also set. 
+
+implicit none
+
+#include "netcdf.inc"
+
+integer,intent(in) :: nid ! NetCDF file ID
+character(len=*),intent(in) :: name ! the variable's name
+character(len=*),intent(in) :: title ! 'title' attribute of variable
+character(len=*),intent(in) :: units ! 'units' attribute of variable
+integer,intent(in) :: nbdim ! number of dimensions of the variable
+integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
+                                              ! the variable is defined along
+integer,intent(out) :: nvarid ! NetCDF ID of the variable
+integer,intent(out) :: ierr ! returned NetCDF staus code
+
+! 1. Switch to NetCDF define mode 
+ierr=NF_REDEF(nid)
+
+! 2. Define the variable
+#ifdef NC_DOUBLE
+ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
+#else
+ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
+#endif
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+! 3. Write attributes
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
+                     len_trim(adjustl(title)),adjustl(title))
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
+                     len_trim(adjustl(units)),adjustl(units))
+if(ierr/=NF_NOERR) then
+   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
+   write(*,*) NF_STRERROR(ierr)
+   stop ""
+endif
+
+! 4. Switch out of NetCDF define mode
+ierr = NF_ENDDEF(nid)
+
+end subroutine def_var_stats
+

trunk/LMDZ.GENERIC/makegcm_g95

@@ -15 +15 @@
 set bands="32x36"
 set scatterers="1"
+set full=""
 ########################################################################
 # path a changer contenant les sources et les objets du modele
@@ -273 +274 @@
              latitudes and vertical layers respectively.
 
--t ntrac   Selects the number of tracers present in the model
-
-             Options -d and -t overwrite file 
-             $LMDGCM/libf/grid/dimensions.h
-             which contains the 3 dimensions of the
-             horizontal grid 
-             im, jm, lm plus the number of tracers passively advected
-             by the dynamics ntrac,
-             in 4 PARAMETER FORTRAN format 
-             with a new file:
-             $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
-             If the file does not exist already
-             it is created by the script
-             $LMDGCM/libf/grid/dimension/makdim
-
 -s nscat   Number of radiatively active scatterers
 
@@ -331 +317 @@
            there is no need to specify  -Ldirn.
 
+-full      Full (re)compilation (from scratch)
+
 eod
 exit
@@ -367 +355 @@
      case -olddyn
         set dyntype="olddyn" ; shift; goto top
+
+     case -full
+        set full="full" ; shift ; goto top
 
      case -filtre
@@ -527 +518 @@
 # Build the appropriate 'dimensions.h' file
 cd dimension
-./makdim $ntrac $dim
+./makdim $dim
 # echo contents of dimensions.h to standard output
 cat $libf/grid/dimensions.h
@@ -559 +550 @@
 echo dimension $dimension dim $dim
 if ( $dimension == 1 ) then
-  echo pas de dynamique
-  set dyn="L_DYN= DYN= L_FILTRE= "
+  echo "No dynamics"
+##  set dyn="L_DYN= DYN= L_FILTRE= "
+## NB: we still need to have L_DYN=libdyn3d to reach routines and module
+## objects which are located in dyn3d
+  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
 endif
 endif
@@ -588 +582 @@
 
 echo "dimc $dimc"
+
+#cleanup for a full recompilation, if requested
+if ("$full" == "full") then
+# remove makefile and $libo
+  cd $model
+  \rm -f makefile
+  \rm -rf $libo/*
+endif
 
 ########################################################################

trunk/LMDZ.GENERIC/makegcm_gfortran

@@ -15 +15 @@
 set bands="32x36"
 set scatterers="1"
+set full=""
 ########################################################################
 # path a changer contenant les sources et les objets du modele
@@ -184 +185 @@
    set optim90="-O3 -funroll-loops "
    set optimtru90="-O3 -funroll-loops "
-   set opt_link=" -L$NCDFLIB -lnetcdf -lnetcdff "
+   set opt_link=" -L$NCDFLIB -lnetcdf "
 
    #NB: on gnome -O3 ==> NaNs ...
@@ -274 +275 @@
              latitudes and vertical layers respectively.
 
--t ntrac   Selects the number of tracers present in the model
-
-             Options -d and -t overwrite file 
-             $LMDGCM/libf/grid/dimensions.h
-             which contains the 3 dimensions of the
-             horizontal grid 
-             im, jm, lm plus the number of tracers passively advected
-             by the dynamics ntrac,
-             in 4 PARAMETER FORTRAN format 
-             with a new file:
-             $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
-             If the file does not exist already
-             it is created by the script
-             $LMDGCM/libf/grid/dimension/makdim
-
 -s nscat   Number of radiatively active scatterers
 
@@ -332 +318 @@
            there is no need to specify  -Ldirn.
 
+-full      Full (re)compilation (from scratch)
+
 eod
 exit
@@ -368 +356 @@
      case -olddyn
         set dyntype="olddyn" ; shift; goto top
+
+     case -full
+        set full="full" ; shift ; goto top
 
      case -filtre
@@ -532 +523 @@
 # Build the appropriate 'dimensions.h' file
 cd dimension
-./makdim $ntrac $dim
+./makdim $dim
 # echo contents of dimensions.h to standard output
 cat $libf/grid/dimensions.h
@@ -564 +555 @@
 echo dimension $dimension dim $dim
 if ( $dimension == 1 ) then
-  echo pas de dynamique
-  set dyn="L_DYN= DYN= L_FILTRE= "
+  echo "No dynamics"
+##  set dyn="L_DYN= DYN= L_FILTRE= "
+## NB: we still need to have L_DYN=libdyn3d to reach routines and module
+## objects which are located in dyn3d
+  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
 endif
 endif
@@ -593 +587 @@
 
 echo "dimc $dimc"
+
+#cleanup for a full recompilation, if requested
+if ("$full" == "full") then
+# remove makefile and $libo
+  cd $model
+  \rm -f makefile
+  \rm -rf $libo/*
+endif
 
 ########################################################################
@@ -657 +659 @@
    set f77="gfortran -ffree-line-length-264"
    set f90="gfortran -ffree-line-length-264"
-   set opt_link=" -L$LIBOGCM -L$NCDFLIB -lnetcdff -lnetcdf "
+#   set opt_link=" -L$LIBOGCM -L$NCDFLIB -lnetcdff -lnetcdf "
+   set opt_link=" -L$LIBOGCM -L$NCDFLIB -lnetcdf "
 else if $SUN then
    set f77=f90

trunk/LMDZ.GENERIC/makegcm_ifort

@@ -15 +15 @@
 set bands="32x36"
 set scatterers="1"
+set full=""
 ########################################################################
 # path a changer contenant les sources et les objets du modele
@@ -279 +280 @@
              latitudes and vertical layers respectively.
 
--t ntrac   Selects the number of tracers present in the model
-
-             Options -d and -t overwrite file 
-             $LMDGCM/libf/grid/dimensions.h
-             which contains the 3 dimensions of the
-             horizontal grid 
-             im, jm, lm plus the number of tracers passively advected
-             by the dynamics ntrac,
-             in 4 PARAMETER FORTRAN format 
-             with a new file:
-             $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
-             If the file does not exist already
-             it is created by the script
-             $LMDGCM/libf/grid/dimension/makdim
-
 -s nscat   Number of radiatively active scatterers
 
@@ -337 +323 @@
            there is no need to specify  -Ldirn.
 
+-full      Full (re)compilation (from scratch)
+
 eod
 exit
@@ -373 +361 @@
      case -olddyn
         set dyntype="olddyn" ; shift; goto top
+
+     case -full
+        set full="full" ; shift ; goto top
 
      case -filtre
@@ -534 +525 @@
 # Build the appropriate 'dimensions.h' file
 cd dimension
-./makdim $ntrac $dim
+./makdim $dim
 # echo contents of dimensions.h to standard output
 cat $libf/grid/dimensions.h
@@ -566 +557 @@
 echo dimension $dimension dim $dim
 if ( $dimension == 1 ) then
-  echo pas de dynamique
-  set dyn="L_DYN= DYN= L_FILTRE= "
+  echo "No dynamics"
+##  set dyn="L_DYN= DYN= L_FILTRE= "
+## NB: we still need to have L_DYN=libdyn3d to reach routines and module
+## objects which are located in dyn3d
+  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
 endif
 endif
@@ -595 +589 @@
 
 echo "dimc $dimc"
+
+#cleanup for a full recompilation, if requested
+if ("$full" == "full") then
+# remove makefile and $libo
+  cd $model
+  \rm -f makefile
+  \rm -rf $libo/*
+endif
 
 ########################################################################

trunk/LMDZ.GENERIC/makegcm_pgf90

@@ -15 +15 @@
 set bands="32x36"
 set scatterers="1"
+set full=""
 ########################################################################
 # path a changer contenant les sources et les objets du modele
@@ -271 +272 @@
              latitudes and vertical layers respectively.
 
--t ntrac   Selects the number of tracers present in the model
-
-             Options -d and -t overwrite file 
-             $LMDGCM/libf/grid/dimensions.h
-             which contains the 3 dimensions of the
-             horizontal grid 
-             im, jm, lm plus the number of tracers passively advected
-             by the dynamics ntrac,
-             in 4 PARAMETER FORTRAN format 
-             with a new file:
-             $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
-             If the file does not exist already
-             it is created by the script
-             $LMDGCM/libf/grid/dimension/makdim
-
 -s nscat   Number of radiatively active scatterers
 
@@ -329 +315 @@
            there is no need to specify  -Ldirn.
 
+-full      Full (re)compilation (from scratch)
+
 eod
 exit
@@ -365 +353 @@
      case -olddyn
         set dyntype="olddyn" ; shift; goto top
+
+     case -full
+        set full="full" ; shift ; goto top
 
      case -filtre
@@ -529 +520 @@
 # Build the appropriate 'dimensions.h' file
 cd dimension
-./makdim $ntrac $dim
+./makdim $dim
 # echo contents of dimensions.h to standard output
 cat $libf/grid/dimensions.h
@@ -561 +552 @@
 echo dimension $dimension dim $dim
 if ( $dimension == 1 ) then
-  echo pas de dynamique
-  set dyn="L_DYN= DYN= L_FILTRE= "
+  echo "No dynamics"
+##  set dyn="L_DYN= DYN= L_FILTRE= "
+## NB: we still need to have L_DYN=libdyn3d to reach routines and module
+## objects which are located in dyn3d
+  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
 endif
 endif
@@ -590 +584 @@
 
 echo "dimc $dimc"
+
+#cleanup for a full recompilation, if requested
+if ("$full" == "full") then
+# remove makefile and $libo
+  cd $model
+  \rm -f makefile
+  \rm -rf $libo/*
+endif
 
 ########################################################################