Changeset 1216 for trunk


Ignore:
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:
14 added
42 edited
4 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

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

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

    r787 r1216  
    66c    Auteur :  P. Le Van, F. Hourdin
    77c   .........
    8 
     8      USE infotrac, ONLY: tname, nqtot
    99      IMPLICIT NONE
    1010c=======================================================================
    … …  
    7171#include "comvert.h"
    7272#include "comgeom2.h"
    73 #include "control.h"
    74 
    75 #include "advtrac.h"
     73!#include "control.h"
     74
     75!#include "advtrac.h"
    7676!! this is to get tnom (tracers name)
    7777
    … …  
    8585      REAL pteta(iip1,jjp1,llm)
    8686      REAL pmasse(iip1,jjp1,llm)
    87       REAL pq(iip1,jjp1,llm,nqmx)
     87      REAL pq(iip1,jjp1,llm,nqtot)
    8888      REAL pphis(iip1,jjp1)
    8989      REAL pphi(iip1,jjp1,llm)
    … …  
    9292      REAL pducov(iip1,jjp1,llm)
    9393      REAL pdteta(iip1,jjp1,llm)
    94       REAL pdq(iip1,jjp1,llm,nqmx)
     94      REAL pdq(iip1,jjp1,llm,nqtot)
    9595c
    9696      REAL pw(iip1,jjp1,llm)
    … …  
    103103      REAL pdufi(iip1,jjp1,llm)
    104104      REAL pdhfi(iip1,jjp1,llm)
    105       REAL pdqfi(iip1,jjp1,llm,nqmx)
     105      REAL pdqfi(iip1,jjp1,llm,nqtot)
    106106      REAL pdpsfi(iip1,jjp1)
    107107      logical tracer
    … …  
    116116c
    117117      REAL zufi(ngridmx,llm), zvfi(ngridmx,llm)
    118       REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqmx)
     118      REAL ztfi(ngridmx,llm),zqfi(ngridmx,llm,nqtot)
    119119c
    120120      REAL zvervel(ngridmx,llm)
    121121c
    122122      REAL zdufi(ngridmx,llm),zdvfi(ngridmx,llm)
    123       REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqmx)
     123      REAL zdtfi(ngridmx,llm),zdqfi(ngridmx,llm,nqtot)
    124124      REAL zdpsrf(ngridmx)
    125125c
    … …  
    170170
    171171c
    172       IF (firstcal) THEN
    173          latfi(1)=rlatu(1)
    174          lonfi(1)=0.
    175          DO j=2,jjm
    176             DO i=1,iim
    177                latfi((j-2)*iim+1+i)= rlatu(j)
    178                lonfi((j-2)*iim+1+i)= rlonv(i)
    179             ENDDO
    180          ENDDO
    181          latfi(ngridmx)= rlatu(jjp1)
    182          lonfi(ngridmx)= 0.
    183 
    184          ! build airefi(), mesh area on physics grid
    185          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
    186          ! Poles are single points on physics grid
    187          airefi(1)=airefi(1)*iim
    188          airefi(ngridmx)=airefi(ngridmx)*iim
    189 
    190          CALL inifis(ngridmx,llm,day_ini,daysec,dtphys,
    191      .                latfi,lonfi,airefi,rad,g,r,cpp)
    192       ENDIF
     172!      IF (firstcal) THEN
     173!         latfi(1)=rlatu(1)
     174!         lonfi(1)=0.
     175!         DO j=2,jjm
     176!            DO i=1,iim
     177!               latfi((j-2)*iim+1+i)= rlatu(j)
     178!               lonfi((j-2)*iim+1+i)= rlonv(i)
     179!            ENDDO
     180!         ENDDO
     181!         latfi(ngridmx)= rlatu(jjp1)
     182!         lonfi(ngridmx)= 0.
     183!
     184!         ! build airefi(), mesh area on physics grid
     185!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
     186!         ! Poles are single points on physics grid
     187!         airefi(1)=airefi(1)*iim
     188!         airefi(ngridmx)=airefi(ngridmx)*iim
     189!
     190!         CALL inifis(ngridmx,llm,day_ini,daysec,dtphys,
     191!     .                latfi,lonfi,airefi,rad,g,r,cpp)
     192!      ENDIF
    193193
    194194c
    … …  
    278278c   43.bis Taceurs (en kg/kg)
    279279c   --------------------------
    280       DO iq=1,nqmx
     280      DO iq=1,nqtot
    281281         DO l=1,llm
    282282            zqfi(1,l,iq) = pq(1,1,l,iq)
    … …  
    425425
    426426      CALL physiq (ngridmx,llm,nq,
    427      .     tnom,
     427     .     tname,
    428428     ,     debut,lafin,
    429429     ,     rday_ecri,heure,dtphys,
    … …  
    467467
    468468
    469 c   62. humidite specifique
     469c   62. traceurs
    470470c   ---------------------
    471471
    472       DO iq=1,nqmx
     472      DO iq=1,nqtot
    473473         DO l=1,llm
    474474            DO i=1,iip1

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

    r1214 r1216  
    1 !-----------------------------------------------------------------------
    2 ! INCLUDE 'control.h'
    3 ! For Fortran 77/Fortran 90 compliance always use line continuation
    4 ! symbols '&' in columns 73 and 6
    5 !
    61
    7       COMMON/control/nday,day_step,                                     &
    8      &              iperiod,iconser,idissip,iphysiq ,                   &
    9      &              periodav,ecritphy,anneeref
     2module control_mod
    103
    11       INTEGER   nday,day_step,iperiod,iconser,                          &
    12      &          idissip,iphysiq,anneeref
    13       REAL periodav, ecritphy
     4  implicit none
    145
    15 !-----------------------------------------------------------------------
     6  integer,save :: nday ! # of days to run
     7  integer,save :: day_step ! # of dynamical time steps per day
     8  integer,save :: iperiod  ! make a Matsuno step before avery iperiod-1 LF steps
     9  integer,save :: iconser !
     10  integer,save :: idissip ! apply dissipation every idissip dynamical step
     11  integer,save :: iphysiq ! call physics every iphysiq dynamical steps
     12  integer,save :: anneeref ! reference year # ! not used
     13  real,save :: periodav
     14  integer,save :: ecritphy ! output data in "diagfi.nc" every ecritphy dynamical steps
     15
     16end module control_mod

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

    r1006 r1216  
    3838      USE ioipsl_getincom
    3939      use sponge_mod,only: callsponge,nsponge,mode_sponge,tetasponge
     40      use control_mod,only: nday, day_step, iperiod, anneeref,
     41     &                      iconser, idissip, iphysiq, ecritphy
    4042      IMPLICIT NONE
    4143
    4244#include "dimensions.h"
    4345#include "paramet.h"
    44 #include "control.h"
     46!#include "control.h"
    4547#include "logic.h"
    4648#include "serre.h"

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

    r993 r1216  
    11      SUBROUTINE dynetat0(fichnom,nq,vcov,ucov,
    22     .                    teta,q,masse,ps,phis,time)
     3      use infotrac, only: tname, nqtot
    34      IMPLICIT NONE
    45
    … …  
    3031#include "serre.h"
    3132#include "logic.h"
    32 #include"advtrac.h"
     33!#include"advtrac.h"
    3334
    3435c   Arguments:
    … …  
    328329!           WRITE(str3(2:3),'(i2.2)') iq
    329330!           ierr =  NF_INQ_VARID (nid, str3, nvarid)
    330 ! NB: tracers are now read in using their name ('tnom' from advtrac.h)
    331 !           write(*,*) "  loading tracer:",trim(tnom(iq))
    332            ierr=NF_INQ_VARID(nid,tnom(iq),nvarid)
     331! NB: tracers are now read in using their name ('tname' from infotrac)
     332!           write(*,*) "  loading tracer:",trim(tname(iq))
     333           ierr=NF_INQ_VARID(nid,tname(iq),nvarid)
    333334           IF (ierr .NE. NF_NOERR) THEN
    334335!              PRINT*, "dynetat0: Le champ <"//str3//"> est absent"
    335               PRINT*, "dynetat0: Le champ <"//trim(tnom(iq))//
     336              PRINT*, "dynetat0: Le champ <"//trim(tname(iq))//
    336337     &                "> est absent"
    337338              PRINT*, "          Il est donc initialise a zero"
    … …  
    346347             IF (ierr .NE. NF_NOERR) THEN
    347348!                 PRINT*, "dynetat0: Lecture echouee pour "//str3
    348                PRINT*, "dynetat0: Lecture echouee pour "//trim(tnom(iq))
     349               PRINT*,"dynetat0: Lecture echouee pour "//trim(tname(iq))
    349350               CALL abort
    350351             ENDIF
    … …  
    354355c        case when new tracer are added in addition to old ones
    355356             write(*,*)'tracers 1 to ', nqold,'were already present'
    356              write(*,*)'tracers ', nqold+1,' to ', nqmx,'are new'
     357             write(*,*)'tracers ', nqold+1,' to ', nqtot,'are new'
    357358             write(*,*)' and initialized to zero'
    358              q(:,:,:,nqold+1:nqmx)=0.0
     359             q(:,:,:,nqold+1:nqtot)=0.0
    359360!             yes=' '
    360361!            do while ((yes.ne.'y').and.(yes.ne.'n'))

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

    r993 r1216  
    11      SUBROUTINE dynredem0(fichnom,idayref,anneeref,phis,nq)
     2      use infotrac, only: tname
    23      IMPLICIT NONE
    34c=======================================================================
    … …  
    1617#include "netcdf.inc"
    1718#include "serre.h"
    18 #include "advtrac.h"
     19!#include "advtrac.h"
    1920c   Arguments:
    2021c   ----------
    … …  
    902903!               ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 12,
    903904!     .                          "Traceurs "//str3)
    904              txt="Traceur "//trim(tnom(iq))
    905 #ifdef NC_DOUBLE
    906                ierr=NF_DEF_VAR(nid,tnom(iq),NF_DOUBLE,4,dims4,nvarid)
    907 #else
    908                ierr=NF_DEF_VAR(nid,tnom(iq),NF_FLOAT,4,dims4,nvarid)
     905             txt="Traceur "//trim(tname(iq))
     906#ifdef NC_DOUBLE
     907               ierr=NF_DEF_VAR(nid,tname(iq),NF_DOUBLE,4,dims4,nvarid)
     908#else
     909               ierr=NF_DEF_VAR(nid,tname(iq),NF_FLOAT,4,dims4,nvarid)
    909910#endif
    910911               ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",
    … …  
    954955      SUBROUTINE dynredem1(fichnom,time,
    955956     .                     vcov,ucov,teta,q,nq,masse,ps)
     957      use infotrac, only: nqtot, tname
    956958      IMPLICIT NONE
    957959c=================================================================
    … …  
    963965#include "comvert.h"
    964966#include "comgeom.h"
    965 #include"advtrac.h"
     967!#include"advtrac.h"
    966968
    967969      INTEGER nq, l
    … …  
    969971      REAL teta(ip1jmp1,llm)                   
    970972      REAL ps(ip1jmp1),masse(ip1jmp1,llm)                   
    971       REAL q(iip1,jjp1,llm,nqmx)
     973      REAL q(iip1,jjp1,llm,nqtot)
    972974      REAL q3d(iip1,jjp1,llm) !temporary variable
    973975      CHARACTER*(*) fichnom
    … …  
    10521054!            WRITE(str3(2:3),'(i2.2)') iq
    10531055!            ierr = NF_INQ_VARID(nid, str3, nvarid)
    1054             ierr=NF_INQ_VARID(nid,tnom(iq),nvarid)
     1056            ierr=NF_INQ_VARID(nid,tname(iq),nvarid)
    10551057            IF (ierr .NE. NF_NOERR) THEN
    10561058!               PRINT*, "Variable "//str3//" n est pas definie"
    1057               PRINT*, "Variable "//trim(tnom(iq))//" n est pas definie"
     1059              PRINT*,"Variable "//trim(tname(iq))//" n est pas definie"
    10581060              CALL abort
    10591061            ENDIF

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

    r1006 r1216  
    11      PROGRAM gcm
    22
     3      use infotrac, only: iniadvtrac, nqtot, iadv
    34      use sponge_mod,only: callsponge,mode_sponge,sponge
     5      use control_mod, only: nday, day_step, iperiod, iphysiq,
     6     &                       iconser, ecritphy, idissip
     7      use comgeomphy, only: initcomgeomphy
    48      IMPLICIT NONE
    59
    … …  
    4246#include "logic.h"
    4347#include "temps.h"
    44 #include "control.h"
     48!#include "control.h"
    4549#include "ener.h"
    4650#include "netcdf.inc"
    4751#include "serre.h"
    4852#include "tracstoke.h"
    49 #include"advtrac.h"
     53!#include"advtrac.h"
    5054
    5155      INTEGER*4  iday ! jour julien
    … …  
    5660      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    5761      real, dimension(ip1jmp1,llm) :: teta   ! temperature potentielle
    58       REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
     62      REAL,allocatable :: q(:,:,:)           ! champs advectes
    5963      REAL ps(ip1jmp1)                       ! pression  au sol
    6064      REAL pext(ip1jmp1)                     ! pression  extensive
    … …  
    7983c   tendances dynamiques
    8084      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
    81       REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
     85      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
     86      REAL,ALLOCATABLE :: dq(:,:,:)
    8287
    8388c   tendances de la dissipation
    … …  
    8792c   tendances physiques
    8893      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
    89       REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
     94      REAL dhfi(ip1jmp1,llm),dpfi(ip1jmp1)
     95      REAL,ALLOCATABLE :: dqfi(:,:,:)
    9096
    9197c   variables pour le fichier histoire
    … …  
    123129      LOGICAL tracer
    124130          data tracer/.true./
    125       INTEGER nq
     131!      INTEGER nq
    126132
    127133C Calendrier
    … …  
    150156      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
    151157
     158c-----------------------------------------------------------------------
     159c    variables pour l'initialisation de la physique :
     160c    ------------------------------------------------
     161      INTEGER ngridmx
     162      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
     163      REAL zcufi(ngridmx),zcvfi(ngridmx)
     164      REAL latfi(ngridmx),lonfi(ngridmx)
     165      REAL airefi(ngridmx)
     166      SAVE latfi, lonfi, airefi
     167      INTEGER i,j
    152168
    153169c-----------------------------------------------------------------------
    … …  
    159175
    160176c-----------------------------------------------------------------------
    161 c  Initialize tracers using iniadvtrac (Ehouarn, oct 2008)
    162       CALL iniadvtrac(nq,numvanle)
    163 
    164 
    165       CALL dynetat0("start.nc",nqmx,vcov,ucov,
     177      CALL defrun_new( 99, .TRUE. )
     178
     179!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     180! FH 2008/05/02
     181! A nettoyer. On ne veut qu'une ou deux routines d'interface
     182! dynamique -> physique pour l'initialisation
     183!#ifdef CPP_PHYS
     184      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     185      call initcomgeomphy
     186!#endif
     187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     188
     189! Initialize tracers
     190      CALL iniadvtrac(nqtot,numvanle)
     191! Allocation de la tableau q : champs advectes   
     192      allocate(q(ip1jmp1,llm,nqtot))
     193      allocate(dq(ip1jmp1,llm,nqtot))
     194      allocate(dqfi(ip1jmp1,llm,nqtot))
     195
     196      CALL dynetat0("start.nc",nqtot,vcov,ucov,
    166197     .              teta,q,masse,ps,phis, time_0)
    167 
    168       CALL defrun_new( 99, .TRUE. )
    169198
    170199c  on recalcule eventuellement le pas de temps
    … …  
    196225     *                tetagdiv, tetagrot , tetatemp              )
    197226c
     227
     228c-----------------------------------------------------------------------
     229c   Initialisation de la physique :
     230c   -------------------------------
     231
     232!      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
     233         latfi(1)=rlatu(1)
     234         lonfi(1)=0.
     235         zcufi(1) = cu(1)
     236         zcvfi(1) = cv(1)
     237         DO j=2,jjm
     238            DO i=1,iim
     239               latfi((j-2)*iim+1+i)= rlatu(j)
     240               lonfi((j-2)*iim+1+i)= rlonv(i)
     241               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
     242               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
     243            ENDDO
     244         ENDDO
     245         latfi(ngridmx)= rlatu(jjp1)
     246         lonfi(ngridmx)= 0.
     247         zcufi(ngridmx) = cu(ip1jm+1)
     248         zcvfi(ngridmx) = cv(ip1jm-iim)
     249
     250         ! build airefi(), mesh area on physics grid
     251         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
     252         ! Poles are single points on physics grid
     253         airefi(1)=airefi(1)*iim
     254         airefi(ngridmx)=airefi(ngridmx)*iim
     255
     256! Initialisation de la physique: pose probleme quand on tourne
     257! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
     258! Il faut une cle CPP_PHYS
     259!#ifdef CPP_PHYS
     260!         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
     261         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys,
     262     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
     263     &                1)
     264!     &                iflag_phys)
     265!#endif
     266!         call_iniphys=.false.
     267!      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
    198268
    199269      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    … …  
    229299     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
    230300
    231       CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
     301      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqtot)
    232302
    233303      ecripar = .TRUE.
    … …  
    237307
    238308c   Quelques initialisations pour les traceurs
    239       call initial0(ijp1llm*nqmx,dq)
     309      call initial0(ijp1llm*nqtot,dq)
    240310c     istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
    241311c     istphy=istdyn/iphysiq
    … …  
    328398       IF( forward. OR . leapf )  THEN
    329399
    330         DO iq = 1, nqmx
     400        DO iq = 1, nqtot
    331401c
    332402         IF ( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
    333403            CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
    334404
    335          ELSE IF( iq.EQ. nqmx )   THEN
     405         ELSE IF( iq.EQ. nqtot )   THEN
    336406c
    337407            iapp_tracvl = 5
    … …  
    341411c
    342412
    343             CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
     413            CALL vanleer(numvanle,iapp_tracvl,nqtot,q,pbaru,pbarv,
    344414     *                      p, masse, dq,  iadv(1), teta, pk     )
    345415
    … …  
    413483
    414484
    415         CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
     485        CALL calfis( nqtot, lafin ,rdayvrai,rday_ecri,time  ,
    416486     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
    417487     $     du,dv,dteta,dq,w, dufi,dvfi,dhfi,dqfi,dpfi,tracer)
    … …  
    421491c      ------------------------------
    422492!        if(1.eq.2)then
    423           CALL addfi( nqmx, dtphys, leapf, forward   ,
     493          CALL addfi( nqtot, dtphys, leapf, forward   ,
    424494     $                  ucov, vcov, teta , q   ,ps , masse,
    425495     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
    … …  
    556626c                 iav=0
    557627c              ENDIF
    558 c              CALL writedynav(histaveid, nqmx, itau,vcov ,
     628c              CALL writedynav(histaveid, nqtot, itau,vcov ,
    559629c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
    560630c           ENDIF
    … …  
    569639       CALL test_period ( ucov,vcov,teta,q,p,phis )
    570640       CALL dynredem1("restart.nc",0.0,
    571      .                     vcov,ucov,teta,q,nqmx,masse,ps)
     641     .                     vcov,ucov,teta,q,nqtot,masse,ps)
    572642
    573643              CLOSE(99)
    … …  
    636706                  iav=0
    637707               ENDIF
    638 c              CALL writedynav(histaveid, nqmx, itau,vcov ,
     708c              CALL writedynav(histaveid, nqtot, itau,vcov ,
    639709c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
    640710
    … …  
    644714                 IF(itau.EQ.itaufin)
    645715     . CALL dynredem1("restart.nc",0.0,
    646      .                     vcov,ucov,teta,q,nqmx,masse,ps)
     716     .                     vcov,ucov,teta,q,nqtot,masse,ps)
    647717
    648718                 forward = .TRUE.

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

    r1214 r1216  
     1MODULE infotrac
     2
     3IMPLICIT NONE
     4! nqtot : total number of tracers and higher order of moment, water vapor and liquid included
     5  INTEGER, SAVE :: nqtot
     6  INTEGER,allocatable :: iadv(:)   ! tracer advection scheme number
     7  CHARACTER(len=20),allocatable ::  tname(:) ! tracer name
     8
     9CONTAINS
     10
    111      subroutine iniadvtrac(nq,numvanle)
    212!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    … …  
    818      IMPLICIT NONE
    919
    10 #include "dimensions.h"
    11 #include "advtrac.h"
    12 #include "control.h"
     20!#include "dimensions.h"
     21!#include "advtrac.h"
     22!#include "control.h"
    1323
    1424! routine arguments:
    … …  
    2030      INTEGER :: iq
    2131      INTEGER :: ierr
    22 
    23 
    24       if (nqmx > 0) then
     32      CHARACTER(len=3) :: qname
    2533
    2634! Look for file traceur.def
    27       OPEN(90,file='traceur.def',form='formatted',status='old',
    28      &        iostat=ierr)
     35      OPEN(90,file='traceur.def',form='formatted',status='old', &
     36              iostat=ierr)
    2937      IF (ierr.eq.0) THEN
    3038        write(*,*) "iniadvtrac: Reading file traceur.def"
    … …  
    3543          write(*,*) "   (first line of traceur.def) "
    3644          stop
    37         else
    38           ! check that the number of tracers is indeed nqmx
    39           if (nq.ne.nqmx) then
    40             write(*,*) "iniadvtrac: error, wrong number of tracers:"
    41             write(*,*) "nq=",nq," whereas nqmx=",nqmx
    42             stop
    43           endif
    4445        endif
     46       
     47        ! allocate arrays:
     48        allocate(iadv(nq))
     49        allocate(tname(nq))
    4550       
    4651        ! initialize advection schemes to Van-Leer for all tracers
    … …  
    4954        enddo
    5055       
    51 
    52 
    53 !     MODIFICATION TO TEST OTHER SCHEMES BY RDW
    54 !        do iq=1,nq
    55 !           iadv(iq)=1
    56 !        enddo
    57 !        print*,'IADV SET TO 1 IN iniadvtrac!!!!'
    58 
    5956        do iq=1,nq
    6057        ! minimal version, just read in the tracer names, 1 per line
    61           read(90,*,iostat=ierr) tnom(iq)
     58          read(90,*,iostat=ierr) tname(iq)
    6259          if (ierr.ne.0) then
    6360            write(*,*) 'iniadvtrac: error reading tracer names...'
    … …  
    6562          endif
    6663        enddo !of do iq=1,nq
    67       ELSE
    68         write(*,*) "iniadvtrac: can't find file traceur.def..."
    69         stop
     64        close(90) ! done reading tracer names, close file
    7065      ENDIF ! of IF (ierr.eq.0)
    7166
    72 c  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
    73 c  ...................................................................
    74 c
    75 c     iadv = 1    shema  transport type "humidite specifique LMD" 
    76 c     iadv = 2    shema   amont
    77 c     iadv = 3    shema  Van-leer
    78 c     iadv = 4    schema  Van-leer + humidite specifique
    79 c                        Modif F.Codron
    80 c
    81 c
    82       DO  iq = 1, nqmx
    83        IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
    84      * ,' pour le traceur no ', iq
    85        IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'
    86 
    87      * ,' traceur no ', iq
    88        IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
    89      * ,'le traceur no ', iq
     67!  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
     68!  ...................................................................
     69!
     70!     iadv = 1    shema  transport type "humidite specifique LMD" 
     71!     iadv = 2    shema   amont
     72!     iadv = 3    shema  Van-leer
     73!     iadv = 4    schema  Van-leer + humidite specifique
     74!                        Modif F.Codron
     75!
     76!
     77      DO  iq = 1, nq-1
     78       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'&
     79       ,' pour le traceur no ', iq
     80       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'  &
     81       ,' traceur no ', iq
     82       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour' &
     83       ,'le traceur no ', iq
    9084
    9185       IF( iadv(iq).EQ.4 )  THEN
    92          PRINT *,' Le shema  Van-Leer + humidite specifique ',
    93      * ' est  uniquement pour la vapeur d eau .'
     86         PRINT *,' Le shema  Van-Leer + humidite specifique ',          &
     87       ' est  uniquement pour la vapeur d eau .'
    9488         PRINT *,' Corriger iadv( ',iq, ')  et repasser ! '
    9589         CALL ABORT
    … …  
    9791
    9892       IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 )   THEN
    99         PRINT *,' Erreur dans le choix de iadv (nqmx).Corriger et '
    100      * ,' repasser car  iadv(iq) = ', iadv(iq)
     93        PRINT *,' Erreur dans le choix de iadv (nqtot).Corriger et '    &
     94       ,' repasser car  iadv(iq) = ', iadv(iq)
    10195         CALL ABORT
    10296       ENDIF
    10397      ENDDO
    10498
    105 !!!! AS: compiler complains about iadv(nqmx) when there is nqmx=0
    106 !!!! AS: so I commented those lines and changed nqmx-1 for nqmx above
    107 !       IF( iadv(nqmx).EQ.1 ) PRINT *,' Choix du shema humidite '
    108 !     * ,'specifique pour la vapeur d''eau'
    109 !       IF( iadv(nqmx).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'
    110 !     * ,' vapeur d''eau '
    111 !       IF( iadv(nqmx).EQ.3 ) PRINT *,' Choix du shema  Van-Leer '
    112 !     * ,' pour la vapeur d''eau'
    113 !       IF( iadv(nqmx).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + '
    114 !     * ,' humidite specifique pour la vapeur d''eau'
     99       IF( iadv(nq).EQ.1 ) PRINT *,' Choix du shema humidite '          &
     100       ,'specifique pour la vapeur d''eau'
     101       IF( iadv(nq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'  &
     102       ,' vapeur d''eau '
     103       IF( iadv(nq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer '         &
     104       ,' pour la vapeur d''eau'
     105       IF( iadv(nq).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + '       &
     106       ,' humidite specifique pour la vapeur d''eau'
    115107!
    116 !c
    117 !!       IF( (iadv(nqmx).LE.0).OR.(iadv(nqmx).GT.4) )   THEN
    118 !!     MODIFICATION TO TEST WITHOUT TRACER ADVECTION BY RDW
    119 !       IF( (iadv(nqmx).LT.0).OR.(iadv(nqmx).GT.4) )   THEN
    120 !        PRINT *,' Erreur dans le choix de iadv (nqmx).Corriger et '
    121 !     * ,' repasser car  iadv(nqmx) = ', iadv(nqmx)
    122 !         CALL ABORT
    123 !       ENDIF
     108       IF( (iadv(nq).LE.0).OR.(iadv(nq).GT.4) )   THEN
     109        PRINT *,' Erreur dans le choix de iadv (nqtot).Corriger et '    &
     110       ,' repasser car  iadv(nqtot) = ', iadv(nqtot)
     111         CALL ABORT
     112       ENDIF
    124113
    125114      first = .TRUE.
    126       numvanle = nqmx + 1
    127       DO  iq = 1, nqmx
     115      numvanle = nq + 1
     116      DO  iq = 1, nq
    128117        IF(((iadv(iq).EQ.3).OR.(iadv(iq).EQ.4)).AND.first ) THEN
    129118          numvanle = iq
    … …  
    131120        ENDIF
    132121      ENDDO
    133 c
    134       DO  iq = 1, nqmx
     122!
     123      DO  iq = 1, nq
    135124
    136125      IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
    137           PRINT *,' Il y a discontinuite dans le choix du shema de ',
    138      *    'Van-leer pour les traceurs . Corriger et repasser . '
     126          PRINT *,' Il y a discontinuite dans le choix du shema de ',   &
     127          'Van-leer pour les traceurs . Corriger et repasser . '
    139128           CALL ABORT
    140129      ENDIF
    141130
    142131      ENDDO
    143 c
    144       endif ! of if nqmx > 0
     132!
     133      end subroutine iniadvtrac
    145134
    146       end
     135END MODULE infotrac

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

    r135 r1216  
    11      SUBROUTINE iniconst
    22
     3      use control_mod, only: iphysiq, idissip
    34      IMPLICIT NONE
    45c
    … …  
    1314#include "comconst.h"
    1415#include "temps.h"
    15 #include "control.h"
     16!#include "control.h"
    1617#include "comvert.h"
    1718

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

    r253 r1216  
    88c   -------------
    99
     10      use control_mod, only: idissip, iperiod
    1011      IMPLICIT NONE
    1112#include "dimensions.h"
    … …  
    1415#include "comconst.h"
    1516#include "comvert.h"
    16 #include "control.h"
     17!#include "control.h"
    1718
    1819      LOGICAL lstardis

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

    r253 r1216  
    33
    44      COMMON/logic/ purmats,physic,forward,leapf,apphys,grireg,
    5      *  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,hybrid,autozlevs
     5     *  statcl,conser,apdiss,apdelq,saison,ecripar,fxyhypb,ysinus,
     6     &  hybrid,autozlevs
    67
    78      LOGICAL purmats,physic,forward,leapf,apphys,grireg,statcl,conser,

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

    r135 r1216  
    11      SUBROUTINE test_period ( ucov, vcov, teta, q, p, phis )
     2     
     3      USE infotrac, ONLY: nqtot
     4      IMPLICIT NONE
    25c
    36c     Auteur : P. Le Van 
    … …  
    1417c
    1518      REAL ucov(ip1jmp1,llm), vcov(ip1jm,llm), teta(ip1jmp1,llm) ,
    16      ,      q(ip1jmp1,llm,nqmx), p(ip1jmp1,llmp1), phis(ip1jmp1)
     19     ,      q(ip1jmp1,llm,nqtot), p(ip1jmp1,llmp1), phis(ip1jmp1)
    1720c
    1821c   .....  Variables  locales  .....
    … …  
    5154     
    5255c
    53       DO nq =1, nqmx
     56      DO nq =1, nqtot
    5457        DO l =1, llm
    5558          DO ij = 1, ip1jmp1, iip1

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

    r135 r1216  
    11#!/bin/bash
    22
    3 nqmx=$1
    4 shift
     3#nqmx=$1
     4#shift
    55for i in $* ; do
    66   list=$list.$i
    77done
    8 
    9 fichdim=dimensions${list}.t${nqmx}
     8fichdim=dimensions${list}
     9#fichdim=dimensions${list}.t${nqmx}
    1010
    1111echo $fichdim
    … …  
    1818   jm=$2
    1919   lm=$3
    20    n2=$1
     20#   n2=$1
    2121   ndm=1
    2222
    … …  
    5151!   INCLUDE 'dimensions.h'
    5252!
    53 !   dimensions.h contient les dimensions du modele
    54 !   ndm est tel que iim=2**ndm
    55 !   nqmx est la dimension de la variable traceur q
     53!   dimensions.h contains the model dimensions
     54!   
    5655!-----------------------------------------------------------------------
    5756
    58       INTEGER, parameter :: iim= $im
     57      INTEGER, parameter :: iim=$im
    5958      INTEGER, parameter :: jjm=$jm
    6059      INTEGER, parameter :: llm=$lm
    6160      INTEGER, parameter :: ndm=$ndm
    62 
    63       integer, parameter :: nqmx=$nqmx
    6461
    6562!-----------------------------------------------------------------------

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

    r1194 r1216  
    33! symbols '&' in columns 73 and 6
    44!
    5       COMMON/callkeys/callrad,corrk,calldifv,UseTurbDiff,calladj        &
     5! Group commons according to their type for minimal performance impact
     6
     7      COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj      &
    68     &   , co2cond,callsoil                                             &
    7      &   , season,diurnal,tlocked,rings_shadow,iradia,lwrite            &
    8      &   , iaervar,iddist,topdustref,callstats,calleofdump              &
     9     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
     10     &   , callstats,calleofdump                                        &
    911     &   , enertest                                                     &
    1012     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
    11      &   , Nmix_co2, radfixed, dusttau                                  &
     13     &   , radfixed                                                     &
    1214     &   , meanOLR, specOLR                                             &
    13      &   , kastprof, Tstrat                                             &
    14      &   , nosurf, intheat, flatten, oblate                             &     
    15      &   , newtonian, tau_relax, testradtimes                           &
     15     &   , kastprof                                                     &
     16     &   , nosurf, oblate                                               &     
     17     &   , newtonian, testradtimes                                      &
    1618     &   , check_cpp_match, force_cpp                                   &
    1719     &   , rayleigh                                                     &
    18      &   , stelbbody, stelTbb                                           &
    19      &   , tplanet                                                      &
    20      &   , obs_tau_col_tropo                                            &
    21      &   , obs_tau_col_strato                                           &
    22      &   , pres_bottom_tropo                                            &
    23      &   , pres_top_tropo                                               &
    24      &   , pres_bottom_strato                                           &
    25      &   , pres_top_strato                                              &
    26      &   , size_tropo                                                   &
    27      &   , size_strato                                                  &
    28      &   , startype, Fat1AU                                             &
     20     &   , stelbbody                                                    &
    2921     &   , nearco2cond                                                  &
    30      &   , tracer, mass_redistrib, varactive, varfixed, satval          &
     22     &   , tracer, mass_redistrib, varactive, varfixed                  &
    3123     &   , sedimentation,water,watercond,waterrain                      &
    3224     &   , aeroco2,aeroh2o,aeroh2so4,aeroback2lay                       &
    3325     &   , aerofixco2,aerofixh2o                                        &
    34      &   , hydrology, sourceevol, icetstep, albedosnow                  &
    35      &   , maxicethick, Tsaldiff                                        &
    36      &   , CLFfixval, CLFvarying                                        &
    37      &   , n2mixratio                                                   &
    38      &   , co2supsat                                                    &
    39      &   , cloudlvl                                                     &
    40      &   , pceil                                                        &
    41      &   , Rmean, J2, MassPlanet                                        &
     26     &   , hydrology, sourceevol                                        &
     27     &   , CLFvarying                                                   &
    4228     &   , strictboundcorrk                                                                                     
    4329
     30      COMMON/callkeys_i/iaervar,iddist,iradia,startype
     31     
     32      COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb,     &
     33     &                  Tstrat,tplanet,obs_tau_col_tropo,               &
     34     &                  obs_tau_col_strato,pres_bottom_tropo,           &
     35     &                  pres_top_tropo,pres_bottom_strato,              &
     36     &                  pres_top_strato,size_tropo,size_strato,satval,  &
     37     &                  CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
     38     &                  maxicethick,Tsaldiff,tau_relax,cloudlvl,        &
     39     &                  icetstep,intheat,flatten,Rmean,J2,MassPlanet
     40     
    4441      logical callrad,corrk,calldifv,UseTurbDiff                        &
    4542     &   , calladj,co2cond,callsoil                                     &

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

    r787 r1216  
     1module comsoil_h
    12
    2        module comsoil_h
     3implicit none
     4! nsoilmx : number of subterranean layers
     5!integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
     6!integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
     7  integer, parameter :: nsoilmx = 18
    38
    4        implicit none
     9  real,save,allocatable,dimension(:) :: layer      ! soil layer depths
     10  real,save,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
     11  real,save,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia
     12  real,save :: volcapa    ! soil volumetric heat capacity
     13       ! NB: volcapa is read fromn control(35) from physicq start file
     14       !     in physdem (or set via tabfi, or initialized in
     15       !                 soil_settings.F)
    516
    6        real,allocatable,dimension(:) :: layer      ! soil layer depths
    7        real,allocatable,dimension(:) :: mlayer     ! soil mid-layer depths
    8        real,allocatable,dimension(:,:) :: inertiedat ! soil thermal inertia
    9        real volcapa    ! soil volumetric heat capacitya
     17contains
    1018
    11        end module comsoil_h
     19  subroutine ini_comsoil_h(ngrid)
     20 
     21  implicit none
     22  integer,intent(in) :: ngrid ! number of atmospheric columns
     23 
     24    if (.not.allocated(layer)) allocate(layer(nsoilmx)) !soil layer depths
     25    if (.not.allocated(mlayer)) allocate(mlayer(0:nsoilmx-1)) ! soil mid-layer depths
     26    if (.not.allocated(inertiedat)) allocate(inertiedat(ngrid,nsoilmx)) ! soil thermal inertia
     27 
     28  end subroutine ini_comsoil_h
    1229
     30end module comsoil_h
     31

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

    r869 r1216  
    77
    88      use radinc_h, only : naerkind
    9       use gases_h
     9      use gases_h, only: gfrac, igas_co2
    1010      use radii_mod, only : co2_reffrad
    1111      use aerosol_mod, only : iaero_co2
    12       USE surfdat_h
    13       USE comgeomfi_h
    14       USE tracer_h
     12      USE surfdat_h, only: albedodat, albedice, emisice, emissiv
     13      USE comgeomfi_h, only: lati
     14      USE tracer_h, only: noms, rho_co2
    1515
    1616
    … …  
    453453!     --------------------------------------------------------------------
    454454      DO ig=1,ngrid
    455          IF(ig.GT.ngrid/2+1) THEN
    456             icap=2
     455         IF(lati(ig).LT.0.) THEN
     456            icap=2 ! Southern Hemisphere
    457457         ELSE
    458             icap=1
     458            icap=1 ! Nortnern hemisphere
    459459         ENDIF
    460460

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

    r863 r1216  
    77! nlayermx : number of atmospheric layers
    88      integer, parameter :: nlayermx = llm
    9 ! nsoilmx : number of subterranean layers
     9! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
    1010      !integer, parameter :: nsoilmx = 4 ! for a test
    11       integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
     11      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
    1212      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
    1313!-----------------------------------------------------------------------

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

    r993 r1216  
    3636
    3737      USE comsoil_h
    38  
     38!      use control_mod
    3939      implicit none
    4040
    … …  
    4949#include "logic.h"
    5050#include "serre.h"
    51 #include "control.h"
     51!#include "control.h"
    5252
    5353#include "netcdf.inc"

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

    r1194 r1216  
    1       SUBROUTINE inifis(ngrid,nlayer,
     1      SUBROUTINE inifis(ngrid,nlayer,nq,
    22     $           day_ini,pdaysec,ptimestep,
    33     $           plat,plon,parea,
    … …  
    66      use radinc_h, only : naerkind
    77      use datafile_mod, only: datadir
    8       use comdiurn_h
    9 
    10       !! to be conservative, but why this include here?
    11       use surfdat_h
    12       use comsaison_h
    13 
    14       USE comgeomfi_h
     8      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
     9      use comgeomfi_h, only: long, lati, area, totarea
     10      use comsoil_h, only: ini_comsoil_h
     11      use control_mod, only: ecritphy
    1512
    1613!=======================================================================
    … …  
    6158
    6259
    63       REAL prad,pg,pr,pcpp,pdaysec,ptimestep
     60      REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
    6461 
    65       INTEGER ngrid,nlayer
    66       REAL plat(ngrid),plon(ngrid),parea(ngrid)
     62      INTEGER,INTENT(IN) :: ngrid,nlayer,nq
     63      REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
    6764      integer day_ini
    6865      INTEGER ig,ierr
    … …  
    10299         STOP
    103100      ENDIF
     101
     102      ! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
     103      ! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
     104      call getin("ecritphy",ecritphy)
    104105
    105106! --------------------------------------------------------------
    … …  
    691692         coslon(ig)=cos(plon(ig))
    692693      ENDDO
    693          
     694
    694695      pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
    695696
    696       RETURN
     697      ! allocate "comsoil_h" arrays
     698      call ini_comsoil_h(ngrid)
     699     
    697700      END

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

    r965 r1216  
    11      subroutine inistats(ierr)
    22
    3 #ifdef CPP_PARA
    43      use mod_phys_lmdz_para, only : is_master
    5 #endif
    64      implicit none
    75
    … …  
    1412#include "netcdf.inc"
    1513
    16 #ifndef CPP_PARA
    17       logical,parameter :: is_master=.true.
    18 #endif
    1914      integer,intent(out) :: ierr
    2015      integer :: nid
    … …  
    3227      nsteppd=nint(daysec/dtphys)
    3328      write (*,*) 'nsteppd=',nsteppd
    34       write (*,*) 'istime=',istime
    35 
    36 
    37 
    3829      if (abs(float(nsteppd)-daysec/dtphys).gt.1.e-8*daysec)
    39      ,   stop'Dans Instat:  1jour .ne. n pas physiques'
    40 
     30     &   stop'Dans Instat:  1jour .ne. n pas physiques'
    4131
    4232      if(mod(nsteppd,istime).ne.0)
    43      ,   stop'Dans Instat:  1jour .ne. n*istime pas physiques'
     33     &   stop'Dans Instat:  1jour .ne. n*istime pas physiques'
    4434
    4535      istats=nsteppd/istime
    … …  
    5747      ! only the master needs do this
    5848
    59       ierr = NF_CREATE("stats.nc",NF_CLOBBER,nid)
     49      ierr = NF_CREATE("stats.nc",IOR(NF_CLOBBER,NF_64BIT_OFFSET),nid)
    6050      if (ierr.ne.NF_NOERR) then
    6151         write (*,*) NF_STRERROR(ierr)
    … …  
    7060
    7161      ierr = NF_ENDDEF(nid)
    72       call def_var(nid,"Time","Time",
    73      .            "days since 0000-00-0 00:00:00",1,
    74      .            idim_time,nvarid,ierr)
     62      call def_var_stats(nid,"Time","Time",
     63     &            "days since 0000-00-0 00:00:00",1,
     64     &            idim_time,nvarid,ierr)
    7565! Time is initialised later by mkstats subroutine
    7666
    77       call def_var(nid,"latitude","latitude","degrees_north",1,
    78      .            idim_lat,nvarid,ierr)
     67      call def_var_stats(nid,"latitude","latitude",
     68     &            "degrees_north",1,idim_lat,nvarid,ierr)
    7969#ifdef NC_DOUBLE
    8070      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlatu/pi*180)
    … …  
    8272      ierr = NF_PUT_VAR_REAL (nid,nvarid,rlatu/pi*180)
    8373#endif
    84       call def_var(nid,"longitude","East longitude","degrees_east",1,
    85      .            idim_lon,nvarid,ierr)
     74      call def_var_stats(nid,"longitude","East longitude",
     75     &            "degrees_east",1,idim_lon,nvarid,ierr)
    8676#ifdef NC_DOUBLE
    8777      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,rlonv/pi*180)
    … …  
    10696      ierr = NF_PUT_VAR_REAL (nid,nvarid,pseudoalt)
    10797#endif
    108       call def_var(nid,"aps","hybrid pressure at midlayers"," ",
    109      .            1,idim_llm,nvarid,ierr)
     98      call def_var_stats(nid,"aps","hybrid pressure at midlayers"
     99     & ," ",1,idim_llm,nvarid,ierr)
    110100#ifdef NC_DOUBLE
    111101      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,aps)
    … …  
    114104#endif
    115105
    116       call def_var(nid,"bps","hybrid sigma at midlayers"," ",
    117      .            1,idim_llm,nvarid,ierr)
     106      call def_var_stats(nid,"bps","hybrid sigma at midlayers"
     107     & ," ",1,idim_llm,nvarid,ierr)
    118108#ifdef NC_DOUBLE
    119109      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,bps)
    … …  
    125115
    126116      endif ! of if (is_master)
    127       end
     117      end subroutine inistats
     118

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

    r993 r1216  
    11      SUBROUTINE iniwrite(nid,idayref,phis)
    22
    3       USE comsoil_h
    4 
     3      use comsoil_h, only: mlayer, nsoilmx
    54      IMPLICIT NONE
    65
    … …  
    2322#include "dimensions.h"
    2423#include "paramet.h"
    25 !include "comconst.h"
    2624#include "comcstfi.h"
    2725#include "comvert.h"
    2826#include "comgeom.h"
    29 #include "temps.h"
    3027#include "ener.h"
    3128#include "logic.h"
    3229#include "netcdf.inc"
    3330#include "serre.h"
    34 #include"dimphys.h"
    3531
    3632c   Arguments:
    3733c   ----------
    3834
    39       integer nid        ! NetCDF file ID
    40       INTEGER*4 idayref  ! date (initial date for this run)
    41       REAL phis(ip1jmp1) ! surface geopotential
     35      integer,intent(in) :: nid        ! NetCDF file ID
     36      INTEGER*4,intent(in) :: idayref  ! date (initial date for this run)
     37      real,intent(in) :: phis(ip1jmp1) ! surface geopotential
    4238
    4339c   Local:
    … …  
    5753         tab_cntrl(l)=0.
    5854      ENDDO
    59       tab_cntrl(1)  = FLOAT(iim)
    60       tab_cntrl(2)  = FLOAT(jjm)
    61       tab_cntrl(3)  = FLOAT(llm)
    62       tab_cntrl(4)  = FLOAT(idayref)
     55      tab_cntrl(1)  = real(iim)
     56      tab_cntrl(2)  = real(jjm)
     57      tab_cntrl(3)  = real(llm)
     58      tab_cntrl(4)  = real(idayref)
    6359      tab_cntrl(5)  = rad
    6460      tab_cntrl(6)  = omeg
    … …  
    252248
    253249!-------------------------------
    254 ! (soil) depth variable mlayer() (known from comsoil.h)
     250! (soil) depth variable mlayer() (known from comsoil_h)
    255251!-------------------------------
    256252      ierr=NF_REDEF (nid) ! Enter NetCDF (re-)define mode

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

    r965 r1216  
    11subroutine iniwritesoil(nid,ngrid)
    2 
    3 use comsoil_h, only : inertiedat, mlayer
    4 #ifdef CPP_PARA
    5 use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
    6 use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    7 #endif
    82
    93! initialization routine for 'writediagoil'. Here we create/define
    … …  
    115! (time-independent) parameters.
    126
     7use comsoil_h, only: mlayer, inertiedat, nsoilmx
     8
    139implicit none
    1410
    1511#include"dimensions.h"
    16 #include"dimphys.h"
    1712#include"paramet.h"
    1813#include"comcstfi.h"
    … …  
    3833integer :: i,j,l,ig0
    3934
    40 #ifdef CPP_PARA
    41 ! Added to work in parallel mode
    42 real dx3_glop(klon_glo,nsoilmx)
    43 real dx3_glo(iim,jjp1,nsoilmx) ! to store a global 3D data set
    44 #else
    45 logical,parameter :: is_master=.true.
    46 logical,parameter :: is_mpi_root=.true.
    47 #endif
    48 
    4935! 1. Define the dimensions
    50 if (is_master) then
    5136! Switch to NetCDF define mode
    5237ierr=NF_REDEF(nid)
    … …  
    170155ierr=NF_PUT_VAR_REAL(nid,varid,mlayer)
    171156#endif
    172 ! Note mlayer(0:nsoilmx-1) known from comsoil.h
     157! Note mlayer(0:nsoilmx-1) known from comsoil_h
    173158if (ierr.ne.NF_NOERR) then
    174159  write(*,*)"iniwritesoil: Error, could not write depth variable"
    … …  
    197182! Note no need to write time variable here; it is done in writediagsoil.
    198183
    199 endif ! of if (is_master)
    200 
    201184! 3. Other variables to be included
    202185
    203186! 3.1 mesh area surrounding each horizontal point
    204 if (is_master) then
    205187ierr=NF_REDEF(nid) ! switch to NetCDF define mode
    206188
    … …  
    236218endif
    237219
    238 endif ! of if (is_master)
    239 
    240 
    241220! 3.2 Thermal inertia
    242 if (is_master) then
    243221ierr=NF_REDEF(nid) ! switch to NetCDF define mode
    244222
    … …  
    262240ierr=NF_PUT_ATT_TEXT(nid,varid,"units",len_trim(text),text)
    263241
    264 endif !of if (is_master)
    265 
    266242! Recast data along 'dynamics' grid
    267 #ifdef CPP_PARA
    268   ! gather field on a "global" (without redundant longitude) array
    269   call Gather(inertiedat,dx3_glop)
    270 !$OMP MASTER
    271   if (is_mpi_root) then
    272     call Grid1Dto2D_glo(dx3_glop,dx3_glo)
    273     ! copy dx3_glo() to dx3(:) and add redundant longitude
    274     data3(1:iim,:,:)=dx3_glo(1:iim,:,:)
    275     data3(iip1,:,:)=data3(1,:,:)
    276   endif
    277 !$OMP END MASTER
    278 !$OMP BARRIER
    279 #else
    280 ! Note: inertiedat is known from comsoil.h
     243! Note: inertiedat is known from comsoil_h
     244
    281245do l=1,nsoilmx
    282246  ! handle the poles
    … …  
    294258  enddo
    295259enddo ! of do l=1,nsoilmx
    296 #endif
    297 
    298 ! Write data3 to file
    299 if (is_master) then
    300  ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
    301  ! Write
    302 #ifdef NC_DOUBLE
    303  ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
    304 #else
    305  ierr=NF_PUT_VAR_REAL(nid,varid,data3)
    306 #endif
    307  if (ierr.ne.NF_NOERR) then
     260
     261! Write data2 to file
     262ierr=NF_ENDDEF(nid) ! switch out of NetCDF define mode
     263! Write
     264#ifdef NC_DOUBLE
     265ierr=NF_PUT_VAR_DOUBLE(nid,varid,data3)
     266#else
     267ierr=NF_PUT_VAR_REAL(nid,varid,data3)
     268#endif
     269if (ierr.ne.NF_NOERR) then
    308270  write(*,*)"iniwritesoil: Error, could not write th_inertia variable"
    309  endif
    310 endif ! of if (is_master)
     271endif
    311272
    312273end subroutine iniwritesoil

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

    r787 r1216  
    11program kcm1d
    22
     3  use infotrac, only: nqtot
    34  use radinc_h,      only: NAERKIND
    45  use radcommon_h,   only: L_NSPECTI, L_NSPECTV, sigma
    56  use watercommon_h, only: mH2O
    6   use ioipsl_getincom
    7   use comsaison_h
     7  use ioipsl_getincom, only: getin
     8  use comsaison_h, only: mu0, fract, dist_star
     9!  use control_mod
    810
    911  implicit none
    … …  
    3234#include "comcstfi.h"
    3335#include "planete.h"
    34 #include "control.h"
     36!#include "control.h"
    3537
    3638  ! --------------------------------------------------------------
    … …  
    4446  real plev(nlayermx+1)   ! Intermediate pressure levels [Pa]
    4547  real temp(nlayermx)     ! temperature at the middle of the layers [K]
    46   real q(nlayermx,nqmx)   ! tracer mixing ratio [kg/kg]
    47   real vmr(nlayermx,nqmx) ! tracer mixing ratio [mol/mol]
    48   real qsurf(nqmx)        ! tracer surface budget [kg/kg] ????
     48  real,allocatable :: q(:,:)   ! tracer mixing ratio [kg/kg]
     49  real,allocatable :: vmr(:,:) ! tracer mixing ratio [mol/mol]
     50  real,allocatable :: qsurf(:)        ! tracer surface budget [kg/kg] ????
    4951  real psurf,psurf_n,tsurf     
    5052
    … …  
    7779  logical firstcall,lastcall,global1d
    7880
    79   character*20 nametrac(nqmx)   ! name of the tracer (no need for adv trac common)
     81  character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
    8082
    8183  ! --------------
    … …  
    219221           write(*,*) "   (first line of traceur.def) "
    220222           stop
    221         else
    222            ! check that the number of tracers is indeed nqmx
    223            if (nq.ne.nqmx) then
    224               write(*,*) "kcm1d: error, wrong number of tracers:"
    225               write(*,*) "nq=",nq," whereas nqmx=",nqmx
    226               stop
    227            endif
    228223        endif
     224        nqtot=nq
     225        ! allocate arrays which depend on number of tracers
     226        allocate(nametrac(nq))
     227        allocate(q(nlayermx,nq))
     228        allocate(vmr(nlayermx,nq))
     229        allocate(qsurf(nq))
    229230
    230231        do iq=1,nq

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

    r993 r1216  
    33     &     q,qsurf,surfith,nid)
    44
    5       USE surfdat_h
    6       USE comsoil_h
    7       USE tracer_h
     5!      USE surfdat_h
     6      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, volcapa, inertiedat
     7      USE tracer_h, ONLY: igcm_co2_ice
     8      USE infotrac, ONLY: tname, nqtot
     9!      USE control_mod
    810
    911c=======================================================================
    … …  
    2628#include "dimensions.h"
    2729#include "dimphys.h"
    28 #include "planete.h"
     30!#include "planete.h"
    2931#include "paramet.h"
    3032#include "comconst.h"
    3133#include "comvert.h"
    3234#include "comgeom2.h"
    33 #include "control.h"
    34 #include "logic.h"
     35!#include "control.h"
     36!#include "logic.h"
    3537#include "ener.h"
    3638#include "temps.h"
    3739#include "netcdf.inc"
    38 #include"advtrac.h"
     40!#include"advtrac.h"
    3941c=======================================================================
    4042c   Declarations
    … …  
    7981      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
    8082      REAL h(iip1,jjp1,llm),ps(iip1,jjp1)
    81       REAL q(iip1,jjp1,llm,nqmx),qtot(iip1,jjp1,llm)
     83      REAL q(iip1,jjp1,llm,nqtot),qtot(iip1,jjp1,llm)
    8284
    8385c autre variables dynamique nouvelle grille
    … …  
    9799      REAL co2ice(ngridmx) ! CO2 ice layer
    98100      REAL emis(ngridmx)
    99       REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
     101      REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqtot)
    100102c     REAL phisfi(ngridmx)
    101103
    … …  
    119121      real co2iceS(iip1,jjp1)
    120122      real emisS(iip1,jjp1)
    121       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqmx)
     123      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
    122124
    123125      real ptotal, co2icetotal
    … …  
    283285      allocate(psold(imold+1,jmold+1))
    284286      allocate(phisold(imold+1,jmold+1))
    285       allocate(qold(imold+1,jmold+1,lmold,nqmx))
     287      allocate(qold(imold+1,jmold+1,lmold,nqtot))
    286288      allocate(co2iceold(imold+1,jmold+1))
    287289      allocate(tsurfold(imold+1,jmold+1))
    … …  
    296298      allocate(surfithold(imold+1,jmold+1))
    297299      allocate(mlayerold(nsoilold))
    298       allocate(qsurfold(imold+1,jmold+1,nqmx))
     300      allocate(qsurfold(imold+1,jmold+1,nqtot))
    299301
    300302      allocate(var (imold+1,jmold+1,llm))
    … …  
    703705
    704706! Surface tracers:     
    705       do iq=1,nqmx
     707      do iq=1,nqtot
    706708        ! initialize all surface tracers to zero
    707709        call initial0((jmold+1)*(imold+1), qsurfold(1,1,iq))
    … …  
    709711
    710712
    711 !      print*,'tnom=',tnom
     713!      print*,'tname=',tname
    712714!      print*,'nid',nid
    713715!      print*,'nvarid',nvarid
    714716!      stop
    715717
    716       DO iq=1,nqmx
    717           txt=trim(tnom(iq))//"_surf"
     718      DO iq=1,nqtot
     719          txt=trim(tname(iq))//"_surf"
    718720          if (txt.eq."h2o_vap") then
    719721            ! There is no surface tracer for h2o_vap;
    … …  
    748750        ENDIF
    749751
    750       ENDDO ! of DO iq=1,nqmx
     752      ENDDO ! of DO iq=1,nqtot
    751753
    752754
    … …  
    894896
    895897! Tracers:     
    896       do iq=1,nqmx
     898      do iq=1,nqtot
    897899         call initial0((jmold+1)*(imold+1)*lmold,qold(1,1,1,iq) )
    898900      enddo
    899901
    900       DO iq=1,nqmx
    901         txt=tnom(iq)
     902      DO iq=1,nqtot
     903        txt=tname(iq)
    902904        write(*,*)"lect_start_archive: loading tracer ",trim(txt)
    903905        ierr = NF_INQ_VARID (nid,txt,nvarid)
    … …  
    925927        ENDIF
    926928
    927       ENDDO ! of DO iq=1,nqmx
     929      ENDDO ! of DO iq=1,nqtot
    928930
    929931
    … …  
    12541256c     write(49,*) 'ucov',vcov
    12551257
    1256       if (nqmx .gt. 0) then
     1258      if (nqtot .gt. 0) then
    12571259c traceurs surface
    1258       do iq = 1, nqmx
     1260      do iq = 1, nqtot
    12591261            call interp_horiz(qsurfold(1,1,iq) ,qsurfs(1,1,iq),
    12601262     &                  imold,jmold,iim,jjm,1,
    … …  
    12621264      enddo
    12631265
    1264       call gr_dyn_fi (nqmx,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
     1266      call gr_dyn_fi (nqtot,iim+1,jjm+1,ngridmx,qsurfs,qsurf)
    12651267
    12661268c traceurs 3D
    1267       do  iq = 1, nqmx
     1269      do  iq = 1, nqtot
    12681270            call interp_vert(qold(1,1,1,iq),var,lmold,llm,
    12691271     &        apsold,bpsold,aps,bps,psold,(imold+1)*(jmold+1))
    … …  
    13051307
    13061308c     Periodicite :
    1307       do  iq = 1, nqmx
     1309      do  iq = 1, nqtot
    13081310         do l=1, llm
    13091311            do j = 1, jjp1
    … …  
    13161318! no need to transfer "co2ice" any more; it is in qsurf(igcm_co2_ice)
    13171319
    1318       endif !! if nqmx .ne. 0
     1320      endif !! if nqtot .ne. 0
    13191321
    13201322c-----------------------------------------------------------------------

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

    r965 r1216  
    11subroutine mkstats(ierr)
     2
    23 
    34!
    … …  
    910!  Yann W. july 2003
    1011
    11 #ifdef CPP_PARA
    1212use mod_phys_lmdz_para, only : is_master
    13 #endif
    1413
    1514implicit none
    … …  
    3534integer :: meanid,sdid
    3635!integer, dimension(4) :: dimout
    37 #ifndef CPP_PARA
    38 logical,parameter :: is_master=.true.
    39 #endif
    4036
    4137! Incrementation of count for the last step, which is not done in wstats

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

    r1023 r1216  
    1515c=======================================================================
    1616
    17       USE tracer_h
    18       USE comsoil_h
    19       USE surfdat_h
    20       USE comgeomfi_h
     17      use infotrac, only: iniadvtrac, nqtot, tname
     18      USE tracer_h, ONLY: igcm_co2_ice, igcm_h2o_vap, igcm_h2o_ice
     19      USE comsoil_h, ONLY: nsoilmx, layer, mlayer, inertiedat
     20      USE surfdat_h, ONLY: phisfi, albedodat,
     21     &                     zmea, zstd, zsig, zgam, zthe
     22      USE comgeomfi_h, ONLY: lati, long, area
    2123      use datafile_mod, only: datadir
    2224! to use  'getin'
    2325      USE ioipsl_getincom, only: getin
     26      use control_mod, only: day_step, iphysiq, anneeref
     27      use phyredem, only: physdem0, physdem1
     28      use iostart, only: open_startphy
     29      use comgeomphy, only: initcomgeomphy
    2430      implicit none
    2531
    … …  
    3137#include "comvert.h"
    3238#include "comgeom2.h"
    33 #include "control.h"
     39!#include "control.h"
    3440#include "logic.h"
    3541#include "ener.h"
    … …  
    3844#include "serre.h"
    3945#include "netcdf.inc"
    40 #include "advtrac.h"
     46!#include "advtrac.h"
    4147c=======================================================================
    4248c   Declarations
    … …  
    6470      REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
    6571      REAL phis(iip1,jjp1)
    66       REAL q(iip1,jjp1,llm,nqmx)               ! champs advectes
     72      REAL,ALLOCATABLE :: q(:,:,:,:)               ! champs advectes
    6773
    6874c autre variables dynamique nouvelle grille
    … …  
    9096      REAL emis(ngridmx)        ! surface emissivity
    9197      real emisread             ! added by RW
    92       REAL qsurf(ngridmx,nqmx)
     98      REAL,ALLOCATABLE :: qsurf(:,:)
    9399      REAL q2(ngridmx,nlayermx+1)
    94100!      REAL rnaturfi(ngridmx)
    … …  
    176182      pa     = 0. ! to ensure disaster rather than minor error if we don`t
    177183                  ! make deliberate choice of these values elsewhere.
     184
     185! initialize "serial/parallel" related stuff
     186      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     187      call initcomgeomphy
     188
     189! Load tracer number and names:
     190      call iniadvtrac(nqtot,numvanle)
     191! allocate arrays
     192      allocate(q(iip1,jjp1,llm,nqtot))
     193      allocate(qsurf(ngridmx,nqtot))
    178194
    179195c=======================================================================
    … …  
    243259c  INITIALISATIONS DIVERSES
    244260c=======================================================================
    245 ! Load tracer names:
    246       call iniadvtrac(nq,numvanle)
    247       ! tnom(:) now contains tracer names
    248 ! JL12 we will need the tracer names to read start in dyneta0
    249261
    250262! Initialize global tracer indexes (stored in tracer.h)
    251263! ... this has to be done before phyetat0
    252       call initracer(ngridmx,nqmx,tnom)
     264      call initracer(ngridmx,nqtot,tname)
    253265
    254266! Take care of arrays in common modules
    … …  
    320332        write(*,*) 'Reading file START'
    321333        fichnom = 'start.nc'
    322         CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
     334        CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
    323335     .       ps,phis,time)
    324336
    325337        write(*,*) 'Reading file STARTFI'
    326338        fichnom = 'startfi.nc'
    327         CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqmx,
     339        CALL phyetat0 (ngridmx,fichnom,tab0,Lmodif,nsoilmx,nqtot,
    328340     .        day_ini,time,
    329341     .        tsurf,tsoil,emis,q2,qsurf,   !) ! temporary modif by RDW
    … …  
    387399c-----------------------------------------------------------------------
    388400      if (choix_1.eq.0) then
     401         ! tabfi requires that input file be first opened by open_startphy(fichnom)
     402         fichnom = 'start_archive.nc'
     403         call open_startphy(fichnom)
    389404         call tabfi (ngridmx,nid,Lmodif,tab0,day_ini,lllm,p_rad,
    390405     .            p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    391406      else if (choix_1.eq.1) then
     407         fichnom = 'startfi.nc'
     408         call open_startphy(fichnom)
    392409         Lmodif=1 ! Lmodif set to 1 to allow modifications in phyeta0                           
    393410         call tabfi (ngridmx,nid_fi,Lmodif,tab0,day_ini,lllm,p_rad,
    … …  
    434451      airefi(1)=airefi(1)*iim
    435452      airefi(ngridmx)=airefi(ngridmx)*iim
     453
     454! also initialize various physics flags/settings which might be needed
     455!    (for instance initracer needs to know about some flags, and/or
     456!      'datafile' path may be changed by user)
     457      call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
     458     &                latfi,lonfi,airefi,rad,g,r,cpp)
    436459
    437460c=======================================================================
    … …  
    574597
    575598        read(*,fmt='(a20)') modif
    576         if (modif(1:1) .eq. ' ') exit ! exit loop on changes
     599        if (modif(1:1) .eq. ' ')then
     600         exit ! exit loop on changes
     601        endif
    577602
    578603        write(*,*)
    … …  
    774799         if (yes.eq.'y') then
    775800           write(*,*) 'OK : conservation of tracer total mass'
    776            DO iq =1, nqmx
     801           DO iq =1, nqtot
    777802             DO l=1,llm
    778803               DO j=1,jjp1
    … …  
    827852         do while (yes.eq.'y')
    828853          write(*,*) 'Which tracer name do you want to change ?'
    829           do iq=1,nqmx
    830             write(*,'(i3,a3,a20)')iq,' : ',trim(tnom(iq))
     854          do iq=1,nqtot
     855            write(*,'(i3,a3,a20)')iq,' : ',trim(tname(iq))
    831856          enddo
    832857          write(*,'(a35,i3)')
    833      &            '(enter tracer number; between 1 and ',nqmx
     858     &            '(enter tracer number; between 1 and ',nqtot
    834859          write(*,*)' or any other value to quit this option)'
    835860          read(*,*) iq
    836           if ((iq.ge.1).and.(iq.le.nqmx)) then
    837             write(*,*)'Change tracer name ',trim(tnom(iq)),' to ?'
     861          if ((iq.ge.1).and.(iq.le.nqtot)) then
     862            write(*,*)'Change tracer name ',trim(tname(iq)),' to ?'
    838863            read(*,*) txt
    839             tnom(iq)=txt
     864            tname(iq)=txt
    840865            write(*,*)'Do you want to change another tracer name (y/n)?'
    841866            read(*,'(a)') yes
    … …  
    843868! inapropiate value of iq; quit this option
    844869            yes='n'
    845           endif ! of if ((iq.ge.1).and.(iq.le.nqmx))
     870          endif ! of if ((iq.ge.1).and.(iq.le.nqtot))
    846871         enddo ! of do while (yes.ne.'y')
    847872
    … …  
    851876c          mise a 0 des q (traceurs)
    852877          write(*,*) 'Tracers set to 0 (1.E-30 in fact)'
    853            DO iq =1, nqmx
     878           DO iq =1, nqtot
    854879             DO l=1,llm
    855880               DO j=1,jjp1
    … …  
    862887
    863888c          set surface tracers to zero
    864            DO iq =1, nqmx
     889           DO iq =1, nqtot
    865890             DO ig=1,ngridmx
    866891                 qsurf(ig,iq)=0.
    … …  
    872897        else if (trim(modif).eq.'q=x') then
    873898             write(*,*) 'Which tracer do you want to modify ?'
    874              do iq=1,nqmx
    875                write(*,*)iq,' : ',trim(tnom(iq))
     899             do iq=1,nqtot
     900               write(*,*)iq,' : ',trim(tname(iq))
    876901             enddo
    877              write(*,*) '(choose between 1 and ',nqmx,')'
     902             write(*,*) '(choose between 1 and ',nqtot,')'
    878903             read(*,*) iq
    879              write(*,*)'mixing ratio of tracer ',trim(tnom(iq)),
     904             write(*,*)'mixing ratio of tracer ',trim(tname(iq)),
    880905     &                 ' ? (kg/kg)'
    881906             read(*,*) val
    … …  
    887912               ENDDO
    888913             ENDDO
    889              write(*,*) 'SURFACE value of tracer ',trim(tnom(iq)),
     914             write(*,*) 'SURFACE value of tracer ',trim(tname(iq)),
    890915     &                   ' ? (kg/m2)'
    891916             read(*,*) val
    … …  
    942967             write(*,*) "followed by 2nd, etc. up to top of atmosphere)"
    943968             write(*,*) 'Which tracer do you want to set?'
    944              do iq=1,nqmx
    945                write(*,*)iq,' : ',trim(tnom(iq))
     969             do iq=1,nqtot
     970               write(*,*)iq,' : ',trim(tname(iq))
    946971             enddo
    947              write(*,*) '(choose between 1 and ',nqmx,')'
     972             write(*,*) '(choose between 1 and ',nqtot,')'
    948973             read(*,*) iq
    949              if ((iq.lt.1).or.(iq.gt.nqmx)) then
     974             if ((iq.lt.1).or.(iq.gt.nqtot)) then
    950975               ! wrong value for iq, go back to menu
    951976               write(*,*) "wrong input value:",iq
    … …  
    953978             endif
    954979             ! look for input file 'profile_tracer'
    955              txt="profile_"//trim(tnom(iq))
     980             txt="profile_"//trim(tname(iq))
    956981             open(41,file=trim(txt),status='old',form='formatted',
    957982     &            iostat=ierr)
    … …  
    970995                   q(:,:,l,iq)=profile(l+1)
    971996                 enddo
    972                  write(*,*)'OK, tracer ',trim(tnom(iq)),' initialized ',
    973      &                     'using values from file ',trim(txt)
     997                 write(*,*)'OK, tracer ',trim(tname(iq)),' initialized '
     998     &                    ,'using values from file ',trim(txt)
    974999               else
    9751000                 write(*,*)'problem reading file ',trim(txt),' !'
    976                  write(*,*)'No modifications to tracer ',trim(tnom(iq))
     1001                 write(*,*)'No modifications to tracer ',trim(tname(iq))
    9771002               endif
    9781003             else
    9791004               write(*,*)'Could not find file ',trim(txt),' !'
    980                write(*,*)'No modifications to tracer ',trim(tnom(iq))
     1005               write(*,*)'No modifications to tracer ',trim(tname(iq))
    9811006             endif
    9821007             
    … …  
    15041529
    15051530         
    1506       CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqmx)
    1507       CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqmx,masse,ps)
     1531      CALL dynredem0("restart.nc",day_ini,anneeref,phis,nqtot)
     1532      CALL dynredem1("restart.nc",0.0,vcov,ucov,teta,q,nqtot,masse,ps)
    15081533C
    15091534C Ecriture etat initial physique
    15101535C
    15111536
    1512 !      do ig=1,ngridmx
    1513 !         print*,'surface ice in newstart=',qsurf(ig,igcm_h2o_ice)
    1514 !         print*,'surface liq in newstart=',qsurf(ig,igcm_h2o_vap)
    1515 !      enddo
    1516 !      stop
    1517 
    1518 
    1519       call physdem1(ngridmx,"restartfi.nc",lonfi,latfi,nsoilmx,nqmx,
    1520      .              dtphys,float(day_ini),
    1521      .              time,tsurf,tsoil,emis,q2,qsurf,
    1522      .              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe,
    1523      .              cloudfrac,totalfrac,hice,tnom)
     1537      call physdem0("restartfi.nc",lonfi,latfi,nsoilmx,ngridmx,llm,
     1538     &              nqtot,dtphys,real(day_ini),0.0,
     1539     &              airefi,albfi,ithfi,zmea,zstd,zsig,zgam,zthe)
     1540      call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nq,
     1541     &                dtphys,real(day_ini),
     1542     &                tsurf,tsoil,emis,q2,qsurf,
     1543     &                cloudfrac,totalfrac,hice)
    15241544
    15251545c=======================================================================
    … …  
    15341554     *rage est differente de la valeur parametree llm =',i4//)
    15351555
     1556      write(*,*) "newstart: All is well that ends well."
     1557
    15361558      end
    15371559

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

    r1214 r1216  
    1       SUBROUTINE phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq,
    2      .           day_ini,time,
    3      .           tsurf,tsoil,emis,q2,qsurf,cloudfrac,totcloudfrac,hice)
    4 
    5       USE surfdat_h
    6       USE comgeomfi_h
    7       USE tracer_h
    8 
    9       implicit none
    10 
    11 c======================================================================
    12 c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
    13 c  Adaptation à Mars : Yann Wanherdrick
    14 c Objet: Lecture de l etat initial pour la physique
    15 c======================================================================
     1subroutine phyetat0 (ngrid,fichnom,tab0,Lmodif,nsoil,nq, &
     2                     day_ini,time,tsurf,tsoil, &
     3                     emis,q2,qsurf,cloudfrac,totcloudfrac,hice)
     4
     5  USE infotrac, ONLY: tname
     6  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
     7  use iostart, only: nid_start, open_startphy, close_startphy, &
     8                     get_field, get_var, inquire_field, &
     9                     inquire_dimension, inquire_dimension_length
     10
     11  implicit none
     12
     13!======================================================================
     14! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
     15!  Adaptation à Mars : Yann Wanherdrick
     16! Objet: Lecture de l etat initial pour la physique
     17!======================================================================
    1618#include "netcdf.inc"
    1719#include "dimensions.h"
    … …  
    2022#include "comcstfi.h"
    2123
    22       INTEGER ngrid
    23 c======================================================================
    24       INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
    25       PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
     24!======================================================================
     25!  INTEGER nbsrf !Mars nbsrf a 1 au lieu de 4
     26!  PARAMETER (nbsrf=1) ! nombre de sous-fractions pour une maille
    2627!======================================================================
    2728!  Arguments:
    2829!  ---------
    2930!  inputs:
    30       character*(*) fichnom ! "startfi.nc" file
    31       integer tab0
    32       integer Lmodif
    33       integer nsoil ! # of soil layers
    34       integer nq
    35       integer day_ini
    36       real time
     31  integer,intent(in) :: ngrid
     32  character*(*),intent(in) :: fichnom ! "startfi.nc" file
     33  integer,intent(in) :: tab0
     34  integer,intent(in) :: Lmodif
     35  integer,intent(in) :: nsoil ! # of soil layers
     36  integer,intent(in) :: nq
     37  integer,intent(in) :: day_ini
     38  real,intent(in) :: time
    3739
    3840!  outputs:
    39       real tsurf(ngrid,nbsrf) ! surface temperature
    40       real tsoil(ngrid,nsoil,nbsrf) ! soil temperature
    41       real emis(ngrid) ! surface emissivity
    42       real q2(ngrid, llm+1) !
    43       real qsurf(ngrid,nq) ! tracers on surface
    44 !      real co2ice(ngrid) ! co2 ice cover
    45       real cloudfrac(ngrid,nlayermx)
    46       real hice(ngrid), totcloudfrac(ngrid)
    47 
     41  real,intent(out) :: tsurf(ngrid) ! surface temperature
     42  real,intent(out) :: tsoil(ngrid,nsoil) ! soil temperature
     43  real,intent(out) :: emis(ngrid) ! surface emissivity
     44  real,intent(out) :: q2(ngrid, llm+1) !
     45  real,intent(out) :: qsurf(ngrid,nq) ! tracers on surface
     46! real co2ice(ngrid) ! co2 ice cover
     47  real,intent(out) :: cloudfrac(ngrid,nlayermx)
     48  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
    4849
    4950!======================================================================
    … …  
    5556
    5657      real xmin,xmax ! to display min and max of a field
    57 c
     58!
    5859      INTEGER ig,iq,lmax
    5960      INTEGER nid, nvarid
    … …  
    6566      CHARACTER*2 str2
    6667      CHARACTER*1 yes
    67 c
     68!
    6869      REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec
    6970      INTEGER nqold
    7071
    7172! flag which identifies if 'startfi.nc' file is using old names (qsurf01,...)
    72       logical :: oldtracernames=.false.
     73!      logical :: oldtracernames=.false.
    7374      integer :: count
    7475      character(len=30) :: txt ! to store some text
     76     
     77      INTEGER :: indextime=1 ! index of selected time, default value=1
     78      logical :: found
    7579
    7680!! added variables by AS to allow to read only slices of startfi
    77       real :: toto(ngrid)
    78       integer :: sta   !! subscript in starti where we start to read
    79       integer, dimension(2) :: sta2d
    80       integer, dimension(2) :: siz2d
    81 
    82 c AS: get the cursor that is stored in dimphys.h
    83 c ... this allows to read only a part of startfi horiz grid
    84       sta = cursor
    85 
    86 c
    87 c ALLOCATE ARRAYS IN surfdat_h
    88 c
    89       IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
    90       IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
    91       IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
    92       IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
    93       IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
    94       IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
    95       IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
    96 
    97 c
    98 c Ouvrir le fichier contenant l etat initial:
    99 c
    100 
    101       ierr = NF_OPEN (fichnom, NF_NOWRITE,nid)
    102       IF (ierr.NE.NF_NOERR) THEN
    103         write(6,*)' Pb d''ouverture du fichier '//fichnom
    104         CALL ABORT
    105       ENDIF
    106 
    107 ! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
    108 !                    qsurf02, ...)
    109       count=0
    110       do iq=1,nq
    111         txt= " "
    112         write(txt,'(a5,i2.2)')'qsurf',iq
    113         ierr=NF_INQ_VARID(nid,txt,nvarid)
    114         if (ierr.ne.NF_NOERR) then
    115           ! did not find old tracer name
    116           exit ! might as well stop here
    117         else
    118           ! found old tracer name
    119           count=count+1
    120         endif
    121       enddo
    122       if (count.eq.nq) then
    123         write(*,*) "phyetat0:tracers seem to follow old naming ",
    124      &             "convention (qsurf01,qsurf02,...)"
    125         write(*,*) "   => will work for now ... "
    126         write(*,*) "      but you should run newstart to rename them"
    127         oldtracernames=.true.
    128       endif
    129 
    130 c modifications possibles des variables de tab_cntrl
    131       write(*,*) 'TABFI de phyeta0',Lmodif,tab0
    132       call tabfi (ngrid,nid,Lmodif,tab0,day_ini,lmax,p_rad,
    133      .              p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
    134 
    135 c
    136 c Lecture des latitudes (coordonnees):
    137 c
    138       ierr = NF_INQ_VARID (nid, "latitude", nvarid)
    139       IF (ierr.NE.NF_NOERR) THEN
    140          PRINT*, 'phyetat0: Le champ <latitude> est absent'
    141          CALL abort
    142       ENDIF
    143 #ifdef NC_DOUBLE
    144       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
    145 #else
    146       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
    147 #endif
    148       IF (ierr.NE.NF_NOERR) THEN
    149          PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
    150          CALL abort
    151       ENDIF
    152 c
    153 c Lecture des longitudes (coordonnees):
    154 c
    155       ierr = NF_INQ_VARID (nid, "longitude", nvarid)
    156       IF (ierr.NE.NF_NOERR) THEN
    157          PRINT*, 'phyetat0: Le champ <longitude> est absent'
    158          CALL abort
    159       ENDIF
    160 #ifdef NC_DOUBLE
    161       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
    162 #else
    163       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
    164 #endif
    165       IF (ierr.NE.NF_NOERR) THEN
    166          PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
    167          CALL abort
    168       ENDIF
    169 c
    170 c Lecture des aires des mailles:
    171 c
    172       ierr = NF_INQ_VARID (nid, "area", nvarid)
    173       IF (ierr.NE.NF_NOERR) THEN
    174          PRINT*, 'phyetat0: Le champ <area> est absent'
    175          CALL abort
    176       ENDIF
    177 #ifdef NC_DOUBLE
    178       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
    179 #else
    180       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
    181 #endif
    182       IF (ierr.NE.NF_NOERR) THEN
    183          PRINT*, 'phyetat0: Lecture echouee pour <area>'
    184          CALL abort
    185       ENDIF
    186       xmin = 1.0E+20
    187       xmax = -1.0E+20
    188       xmin = MINVAL(area)
    189       xmax = MAXVAL(area)
    190       PRINT*,'Aires des mailles <area>:', xmin, xmax
    191 c
    192 c Lecture du geopotentiel au sol:
    193 c
    194       ierr = NF_INQ_VARID (nid, "phisfi", nvarid)
    195       IF (ierr.NE.NF_NOERR) THEN
    196          PRINT*, 'phyetat0: Le champ <phisfi> est absent'
    197          CALL abort
    198       ENDIF
    199 #ifdef NC_DOUBLE
    200       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,phisfi)
    201 #else
    202       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,phisfi)
    203 #endif
    204       IF (ierr.NE.NF_NOERR) THEN
    205          PRINT*, 'phyetat0: Lecture echouee pour <phisfi>'
    206          CALL abort
    207       ENDIF
    208       xmin = 1.0E+20
    209       xmax = -1.0E+20
    210       xmin = MINVAL(phisfi)
    211       xmax = MAXVAL(phisfi)
    212       PRINT*,'Geopotentiel au sol <phisfi>:', xmin, xmax
    213 c
    214 c Lecture de l''albedo du sol nu:
    215 c
    216       ierr = NF_INQ_VARID (nid, "albedodat", nvarid)
    217       IF (ierr.NE.NF_NOERR) THEN
    218          PRINT*, 'phyetat0: Le champ <albedodat> est absent'
    219          CALL abort
    220       ENDIF
    221 #ifdef NC_DOUBLE
    222       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,albedodat)
    223 #else
    224       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,albedodat)
    225 #endif
    226       IF (ierr.NE.NF_NOERR) THEN
    227          PRINT*, 'phyetat0: Lecture echouee pour <albedodat>'
    228          CALL abort
    229       ENDIF
    230       xmin = 1.0E+20
    231       xmax = -1.0E+20
    232       xmin = MINVAL(albedodat)
    233       xmax = MAXVAL(albedodat)
    234       PRINT*,'Albedo du sol nu <albedodat>:', xmin, xmax
    235 c
    236 c ZMEA
    237 c
    238       ierr = NF_INQ_VARID (nid, "ZMEA", nvarid)
    239       IF (ierr.NE.NF_NOERR) THEN
    240          PRINT*, 'phyetat0: Le champ <ZMEA> est absent'
    241          CALL abort
    242       ENDIF
    243 #ifdef NC_DOUBLE
    244       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zmea)
    245 #else
    246       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zmea)
    247 #endif
    248       IF (ierr.NE.NF_NOERR) THEN
    249          PRINT*, 'phyetat0: Lecture echouee pour <ZMEA>'
    250          CALL abort
    251       ENDIF
    252       xmin = 1.0E+20
    253       xmax = -1.0E+20
    254       DO i = 1, ngrid
    255          xmin = MIN(zmea(i),xmin)
    256          xmax = MAX(zmea(i),xmax)
    257       ENDDO
    258       PRINT*,'<zmea>:', xmin, xmax
    259 c
    260 c ZSTD
    261 c
    262       ierr = NF_INQ_VARID (nid, "ZSTD", nvarid)
    263       IF (ierr.NE.NF_NOERR) THEN
    264          PRINT*, 'phyetat0: Le champ <ZSTD> est absent'
    265          CALL abort
    266       ENDIF
    267 #ifdef NC_DOUBLE
    268       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zstd)
    269 #else
    270       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zstd)
    271 #endif
    272       IF (ierr.NE.NF_NOERR) THEN
    273          PRINT*, 'phyetat0: Lecture echouee pour <ZSTD>'
    274          CALL abort
    275       ENDIF
    276       xmin = 1.0E+20
    277       xmax = -1.0E+20
    278       DO i = 1, ngrid
    279          xmin = MIN(zstd(i),xmin)
    280          xmax = MAX(zstd(i),xmax)
    281       ENDDO
    282       PRINT*,'<zstd>:', xmin, xmax
    283 c
    284 c ZSIG
    285 c
    286       ierr = NF_INQ_VARID (nid, "ZSIG", nvarid)
    287       IF (ierr.NE.NF_NOERR) THEN
    288          PRINT*, 'phyetat0: Le champ <ZSIG> est absent'
    289          CALL abort
    290       ENDIF
    291 #ifdef NC_DOUBLE
    292       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zsig)
    293 #else
    294       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zsig)
    295 #endif
    296       IF (ierr.NE.NF_NOERR) THEN
    297          PRINT*, 'phyetat0: Lecture echouee pour <ZSIG>'
    298          CALL abort
    299       ENDIF
    300       xmin = 1.0E+20
    301       xmax = -1.0E+20
    302       DO i = 1, ngrid
    303          xmin = MIN(zsig(i),xmin)
    304          xmax = MAX(zsig(i),xmax)
    305       ENDDO
    306       PRINT*,'<zsig>:', xmin, xmax
    307 c
    308 c ZGAM
    309 c
    310       ierr = NF_INQ_VARID (nid, "ZGAM", nvarid)
    311       IF (ierr.NE.NF_NOERR) THEN
    312          PRINT*, 'phyetat0: Le champ <ZGAM> est absent'
    313          CALL abort
    314       ENDIF
    315 #ifdef NC_DOUBLE
    316       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zgam)
    317 #else
    318       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zgam)
    319 #endif
    320       IF (ierr.NE.NF_NOERR) THEN
    321          PRINT*, 'phyetat0: Lecture echouee pour <ZGAM>'
    322          CALL abort
    323       ENDIF
    324       xmin = 1.0E+20
    325       xmax = -1.0E+20
    326       DO i = 1, ngrid
    327          xmin = MIN(zgam(i),xmin)
    328          xmax = MAX(zgam(i),xmax)
    329       ENDDO
    330       PRINT*,'<zgam>:', xmin, xmax
    331 c
    332 c ZTHE
    333 c
    334       ierr = NF_INQ_VARID (nid, "ZTHE", nvarid)
    335       IF (ierr.NE.NF_NOERR) THEN
    336          PRINT*, 'phyetat0: Le champ <ZTHE> est absent'
    337          CALL abort
    338       ENDIF
    339 #ifdef NC_DOUBLE
    340       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,zthe)
    341 #else
    342       ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,zthe)
    343 #endif
    344       IF (ierr.NE.NF_NOERR) THEN
    345          PRINT*, 'phyetat0: Lecture echouee pour <ZTHE>'
    346          CALL abort
    347       ENDIF
    348       xmin = 1.0E+20
    349       xmax = -1.0E+20
    350       DO i = 1, ngrid
    351          xmin = MIN(zthe(i),xmin)
    352          xmax = MAX(zthe(i),xmax)
    353       ENDDO
    354       PRINT*,'<zthe>:', xmin, xmax
    355 c
    356 c CO2 ice cover
    357 c
    358 ! Ehouarn: from now on, there is no "co2ice" standalone field; it is supposed
    359 ! to be with stored in qsurf(i_co2_ice)
    360 !      ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
    361 !      IF (ierr.NE.NF_NOERR) THEN
    362 !         PRINT*, 'phyetat0: Le champ <co2ice> est absent'
     81!      real :: toto(ngrid)
     82!      integer :: sta   !! subscript in starti where we start to read
     83!      integer, dimension(2) :: sta2d
     84!      integer, dimension(2) :: siz2d
     85
     86! AS: get the cursor that is stored in dimphys.h
     87! ... this allows to read only a part of startfi horiz grid
     88!sta = cursor
     89
     90!
     91! ALLOCATE ARRAYS IN surfdat_h
     92!
     93IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(ngrid))
     94IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(ngrid))
     95IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(ngrid))
     96IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(ngrid))
     97IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(ngrid))
     98IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(ngrid))
     99IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(ngrid))
     100
     101! open physics initial state file:
     102call open_startphy(fichnom)
     103
     104
     105! possibility to modify tab_cntrl in tabfi
     106write(*,*)
     107write(*,*) 'TABFI in phyeta0: Lmodif=',Lmodif," tab0=",tab0
     108call tabfi (ngrid,nid_start,Lmodif,tab0,day_ini,lmax,p_rad, &
     109                   p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time)
     110
     111!c
     112!c Lecture des latitudes (coordonnees):
     113!c
     114!      ierr = NF_INQ_VARID (nid, "latitude", nvarid)
     115!      IF (ierr.NE.NF_NOERR) THEN
     116!         PRINT*, 'phyetat0: Le champ <latitude> est absent'
    363117!         CALL abort
    364118!      ENDIF
    365119!#ifdef NC_DOUBLE
    366 !      ierr = NF_GET_VAR_DOUBLE(nid, nvarid, co2ice)
     120!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,lati)
    367121!#else
    368 !      ierr = NF_GET_VAR_REAL(nid, nvarid, co2ice)
     122!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,lati)
    369123!#endif
    370124!      IF (ierr.NE.NF_NOERR) THEN
    371 !         PRINT*, 'phyetat0: Lecture echouee pour <co2ice>'
     125!         PRINT*, 'phyetat0: Lecture echouee pour <latitude>'
     126!         CALL abort
     127!      ENDIF
     128!c
     129!c Lecture des longitudes (coordonnees):
     130!c
     131!      ierr = NF_INQ_VARID (nid, "longitude", nvarid)
     132!      IF (ierr.NE.NF_NOERR) THEN
     133!         PRINT*, 'phyetat0: Le champ <longitude> est absent'
     134!         CALL abort
     135!      ENDIF
     136!#ifdef NC_DOUBLE
     137!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,long)
     138!#else
     139!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,long)
     140!#endif
     141!      IF (ierr.NE.NF_NOERR) THEN
     142!         PRINT*, 'phyetat0: Lecture echouee pour <longitude>'
     143!         CALL abort
     144!      ENDIF
     145!c
     146!c Lecture des aires des mailles:
     147!c
     148!      ierr = NF_INQ_VARID (nid, "area", nvarid)
     149!      IF (ierr.NE.NF_NOERR) THEN
     150!         PRINT*, 'phyetat0: Le champ <area> est absent'
     151!         CALL abort
     152!      ENDIF
     153!#ifdef NC_DOUBLE
     154!      ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,area)
     155!#else
     156!      ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,area)
     157!#endif
     158!      IF (ierr.NE.NF_NOERR) THEN
     159!         PRINT*, 'phyetat0: Lecture echouee pour <area>'
    372160!         CALL abort
    373161!      ENDIF
    374162!      xmin = 1.0E+20
    375163!      xmax = -1.0E+20
    376 !      xmin = MINVAL(co2ice)
    377 !      xmax = MAXVAL(co2ice)
    378 !      PRINT*,'CO2 ice cover <co2ice>:', xmin, xmax
    379 c
    380 c Lecture des temperatures du sol:
    381 c
    382       ierr = NF_INQ_VARID (nid, "tsurf", nvarid)
    383       IF (ierr.NE.NF_NOERR) THEN
    384          PRINT*, 'phyetat0: Le champ <tsurf> est absent'
    385          PRINT*, '          Mais je vais essayer de lire TS**'
    386          IF (nbsrf.GT.99) THEN
    387             PRINT*, "Trop de sous-mailles"
    388             CALL abort
    389          ENDIF
    390          DO nsrf = 1, nbsrf
    391            WRITE(str2,'(i2.2)') nsrf
    392            ierr = NF_INQ_VARID (nid, "TS"//str2, nvarid)
    393            IF (ierr.NE.NF_NOERR) THEN
    394            PRINT*, "phyetat0: Le champ <TS"//str2//"> est absent"
    395               CALL abort
    396            ENDIF
    397 #ifdef NC_DOUBLE
    398            ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsurf(1,nsrf))
    399 #else
    400            ierr = NF_GET_VAR_REAL(nid, nvarid, tsurf(1,nsrf))
    401 #endif
    402            IF (ierr.NE.NF_NOERR) THEN
    403              PRINT*, "phyetat0: Lecture echouee pour <TS"//str2//">"
    404              CALL abort
    405            ENDIF
    406            xmin = 1.0E+20
    407            xmax = -1.0E+20
    408            xmin = MINVAL(tsurf)
    409            xmax = MAXVAL(tsurf)
    410            PRINT*,'Temperature du sol TS**:', nsrf, xmin, xmax
    411          ENDDO
    412       ELSE
    413          PRINT*, 'phyetat0: Le champ <tsurf> est present'
    414          PRINT*, '          J ignore donc les autres temperatures TS**'
    415 #ifdef NC_DOUBLE
    416          ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, tsurf)
    417 #else
    418          ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, tsurf)
    419 #endif
    420          IF (ierr.NE.NF_NOERR) THEN
    421             PRINT*, "phyetat0: Lecture echouee pour <TSURF>"
    422             CALL abort
    423          ENDIF
    424          xmin = 1.0E+20
    425          xmax = -1.0E+20
    426          xmin = MINVAL(tsurf)
    427          xmax = MAXVAL(tsurf)
    428          PRINT*,'Temperature du sol <tsurf>', xmin, xmax
    429          IF (nbsrf >= 2) THEN
    430             DO nsrf = 2, nbsrf
    431                DO i = 1, ngrid
    432                   tsurf(i,nsrf) = tsurf(i,1)
    433                ENDDO
    434             ENDDO
    435          ENDIF
    436       ENDIF
    437 c
    438 c Lecture des temperatures du sol profond:
    439 c
    440 !      IF (nsoil.GT.99 .OR. nbsrf.GT.99) THEN
    441 !         PRINT*, "Trop de couches ou sous-mailles"
    442 !         CALL abort
    443 !      ENDIF
    444 !      DO nsrf = 1, nbsrf
    445 !         DO isoil=1, nsoil
    446 !            WRITE(str7,'(i2.2,"srf",i2.2)') isoil, nsrf
    447 !            ierr = NF_INQ_VARID (nid, 'tsoil', nvarid)
    448 !            IF (ierr.NE.NF_NOERR) THEN
    449 !               PRINT*, "phyetat0: Le champ <tsoil> est absent"
    450 !               PRINT*, "          Il prend donc la valeur de surface"
    451 !               DO i=1, ngrid
    452 !                  tsoil(i,isoil,nsrf)=tsurf(i,nsrf)
    453 !               ENDDO
    454 !            ELSE
    455 !#ifdef NC_DOUBLE
    456 !              ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tsoil(1,1,nsrf))
    457 !#else
    458 !              ierr = NF_GET_VAR_REAL(nid, nvarid, tsoil(1,1,nsrf))
    459 !#endif
    460 !               IF (ierr.NE.NF_NOERR) THEN
    461 !                  PRINT*, "Lecture echouee pour <tsoil>"
    462 !                  CALL abort
    463 !               ENDIF
    464 !            ENDIF
    465 !         ENDDO
    466 !      ENDDO
    467 !      xmin = 1.0E+20
    468 !      xmax = -1.0E+20
    469 !      xmin = MINVAL(tsoil)
    470 !      xmax = MAXVAL(tsoil)
    471 !      PRINT*,'Temperatures du sol profond <tsoil>', xmin, xmax
    472 c
    473 c Surface emissivity
    474 c
    475       ierr = NF_INQ_VARID (nid, "emis", nvarid)
    476       IF (ierr.NE.NF_NOERR) THEN
    477          PRINT*, 'phyetat0: Le champ <emis> est absent'
    478          CALL abort
    479       ENDIF
    480 #ifdef NC_DOUBLE
    481       ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta,ngrid, emis)
    482 #else
    483       ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, emis)
    484 #endif
    485       IF (ierr.NE.NF_NOERR) THEN
    486          PRINT*, 'phyetat0: Lecture echouee pour <emis>'
    487          CALL abort
    488       ENDIF
    489       xmin = 1.0E+20
    490       xmax = -1.0E+20
    491       xmin = MINVAL(emis)
    492       xmax = MAXVAL(emis)
    493       PRINT*,'Surface emissivity <emis>:', xmin, xmax
    494 
    495 c
    496 c Cloud fraction (added by BC 2010)
    497 c
    498       ierr = NF_INQ_VARID (nid, "cloudfrac", nvarid)
    499       IF (ierr.NE.NF_NOERR) THEN
    500          PRINT*, 'phyetat0: Le champ <cloudfrac> est absent'
    501       cloudfrac(:,:)=0.5
    502 !         CALL abort
    503       ENDIF
    504       sta2d = (/sta,1/)
    505       siz2d = (/ngrid, nlayermx/)
    506 #ifdef NC_DOUBLE
    507       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,cloudfrac)
    508 #else
    509       ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,cloudfrac)
    510 #endif
    511       IF (ierr.NE.NF_NOERR) THEN
    512          PRINT*, 'phyetat0: Lecture echouee pour <cloudfrac>'
    513          CALL abort
    514       ENDIF
    515       xmin = 1.0E+20
    516       xmax = -1.0E+20
    517       xmin = MINVAL(cloudfrac)
    518       xmax = MAXVAL(cloudfrac)
    519       PRINT*,'Cloud fraction <cloudfrac>:', xmin, xmax
    520 
    521 
    522 c
    523 c Total cloud fraction (added by BC 2010)
    524 c
    525       ierr = NF_INQ_VARID (nid, "totcloudfrac", nvarid)
    526 !      ierr = NF_INQ_VARID (nid, "totalfrac", nvarid)
    527       IF (ierr.NE.NF_NOERR) THEN
    528          PRINT*, 'phyetat0: Le champ <totcloudfrac> est absent'
    529       totcloudfrac(:)=0.5
    530 !         CALL abort
    531       ENDIF
    532 #ifdef NC_DOUBLE
    533       ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, totcloudfrac)
    534 #else
    535       ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, totcloudfrac)
    536 #endif
    537       IF (ierr.NE.NF_NOERR) THEN
    538          PRINT*, 'phyetat0: Lecture echouee pour <totcloudfrac>'
    539          !CALL abort
    540       ENDIF
    541       xmin = 1.0E+20
    542       xmax = -1.0E+20
    543       xmin = MINVAL(totcloudfrac)
    544       xmax = MAXVAL(totcloudfrac)
    545       PRINT*,'Cloud fraction <totcloudfrac>:', xmin, xmax
    546 
    547 
    548 
    549 c
    550 c Height of oceanic ice (added by BC 2010)
    551 c
    552       ierr = NF_INQ_VARID (nid, "hice", nvarid)
    553       IF (ierr.NE.NF_NOERR) THEN
    554          PRINT*, 'phyetat0: Le champ <hice> est absent'
    555          CALL abort
    556       ENDIF
    557 #ifdef NC_DOUBLE
    558       ierr = NF_GET_VARA_DOUBLE(nid, nvarid,sta,ngrid, hice)
    559 #else
    560       ierr = NF_GET_VARA_REAL(nid, nvarid,sta,ngrid, hice)
    561 #endif
    562       IF (ierr.NE.NF_NOERR) THEN
    563          PRINT*, 'phyetat0: Lecture echouee pour <hice>'
    564          CALL abort
    565       ENDIF
    566       xmin = 1.0E+20
    567       xmax = -1.0E+20
    568       xmin = MINVAL(hice)
    569       xmax = MAXVAL(hice)
    570       PRINT*,'Height of oceanic ice <hice>:', xmin, xmax
    571 
    572 
    573 c
    574 c pbl wind variance
    575 c
    576       ierr = NF_INQ_VARID (nid, "q2", nvarid)
    577       IF (ierr.NE.NF_NOERR) THEN
    578          PRINT*, 'phyetat0: Le champ <q2> est absent'
    579          CALL abort
    580       ENDIF
    581       sta2d = (/sta,1/)
    582       siz2d = (/ngrid, llm+1/)
    583 #ifdef NC_DOUBLE
    584       ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,q2)
    585 #else
    586       ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,q2)
    587 #endif
    588       IF (ierr.NE.NF_NOERR) THEN
    589          PRINT*, 'phyetat0: Lecture echouee pour <q2>'
    590          CALL abort
    591       ENDIF
    592       xmin = 1.0E+20
    593       xmax = -1.0E+20
    594       xmin = MINVAL(q2)
    595       xmax = MAXVAL(q2)
    596       PRINT*,'pbl wind variance <q2>:', xmin, xmax
    597 c
    598 c tracer on surface
    599 c
    600 
    601       IF(nq.GE.1) THEN
    602          nqold=nq
    603          DO iq=1,nq
    604 !            str7(1:5)='qsurf'
    605 !            WRITE(str7(6:7),'(i2.2)') iq
    606 !            ierr = NF_INQ_VARID (nid,str7,nvarid)
    607            IF (oldtracernames) THEN
    608              txt=" "
    609              write(txt,'(a5,i2.2)')'qsurf',iq
    610            ELSE
    611              txt=noms(iq)
    612 !             if (txt.eq."h2o_vap") then
    613                ! There is no surface tracer for h2o_vap;
    614                ! "h2o_ice" should be loaded instead
    615 !               txt="h2o_ice"
    616 !               write(*,*) 'phyetat0: loading surface tracer',
    617 !     &                     ' h2o_ice instead of h2o_vap'
    618 !             endif
    619            ENDIF ! of IF (oldtracernames) THEN
    620            ierr=NF_INQ_VARID(nid,txt,nvarid)
    621            IF (ierr.NE.NF_NOERR) THEN
    622              write(*,*) 'PHYETAT0: WARNING : surface tracer',trim(txt),
    623      &                  ' not found in file'
    624              write(*,*) trim(txt), ' set to 0'
    625              do ig=1,ngrid
    626                qsurf(ig,iq)=0.
    627              end do
    628              nqold=min(iq-1,nqold)
    629            ELSE
    630            toto(:) = 0.
    631 #ifdef NC_DOUBLE
    632            ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta,ngrid,toto)
    633 #else
    634            ierr = NF_GET_VARA_REAL(nid,nvarid,sta,ngrid,toto)
    635 #endif
    636              IF (ierr.NE.NF_NOERR) THEN
    637                PRINT*, 'phyetat0: Lecture echouee pour <',trim(txt),'>'
    638                CALL abort
    639              ENDIF
    640              qsurf(1:ngrid,iq) = toto(1:ngrid)
    641            ENDIF
    642            xmin = 1.0E+20
    643            xmax = -1.0E+20
    644            xmin = MINVAL(qsurf(1:ngrid,iq))
    645            xmax = MAXVAL(qsurf(1:ngrid,iq))
    646            PRINT*,'tracer on surface <',trim(txt),'>:',xmin,xmax
    647          ENDDO
    648          if ((nqold.lt.nq).and.(nqold.ge.1)) then
    649 c        case when new tracer are added in addition to old ones
    650              write(*,*)'qsurf 1 to ', nqold,'were already present'
    651              write(*,*)'qsurf ', nqold+1,' to ', nq,'are new'
    652              write(*,*)' and initialized to zero'
    653              qsurf(:,nqold+1:nq)=0.0
    654 !            yes=' '
    655 !            do while ((yes.ne.'y').and.(yes.ne.'n'))
    656 !             write(*,*) 'Would you like to reindex qsurf # 1 ->',nqold
    657 !             write(*,*) 'to #',nq-nqold+1,'->', nq,'   (y or n) ?'
    658 !             read(*,fmt='(a)') yes
    659 !            end do
    660 !            if (yes.eq.'y') then
    661 !              write(*,*) 'OK, let s reindex qsurf'
    662 !                 do ig=1,ngrid
    663 !                    do iq=nq,nq-nqold+1,-1
    664 !                       qsurf(ig,iq)=qsurf(ig,iq-nq+nqold)
    665 !                    end do
    666 !                    do iq=nq-nqold,1,-1
    667 !                       qsurf(ig,iq)= 0.
    668 !                    end do
    669 !                 end do
    670 !            end if
    671          end if ! of if ((nqold.lt.nq).and.(nqold.ge.1))
    672       ENDIF ! of IF(nq.GE.1)
     164!      xmin = MINVAL(area)
     165!      xmax = MAXVAL(area)
     166!      PRINT*,'Aires des mailles <area>:', xmin, xmax
     167
     168! Load surface geopotential:
     169call get_field("phisfi",phisfi,found)
     170if (.not.found) then
     171  write(*,*) "phyetat0: Failed loading <phisfi>"
     172  call abort
     173else
     174  write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
     175             minval(phisfi), maxval(phisfi)
     176endif
     177
     178! Load bare ground albedo:
     179call get_field("albedodat",albedodat,found)
     180if (.not.found) then
     181  write(*,*) "phyetat0: Failed loading <albedodat>"
     182  call abort
     183else
     184  write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
     185             minval(albedodat), maxval(albedodat)
     186endif
     187
     188! ZMEA
     189call get_field("ZMEA",zmea,found)
     190if (.not.found) then
     191  write(*,*) "phyetat0: Failed loading <ZMEA>"
     192  call abort
     193else
     194  write(*,*) "phyetat0: <ZMEA> range:", &
     195             minval(zmea), maxval(zmea)
     196endif
     197
     198! ZSTD
     199call get_field("ZSTD",zstd,found)
     200if (.not.found) then
     201  write(*,*) "phyetat0: Failed loading <ZSTD>"
     202  call abort
     203else
     204  write(*,*) "phyetat0: <ZSTD> range:", &
     205             minval(zstd), maxval(zstd)
     206endif
     207
     208! ZSIG
     209call get_field("ZSIG",zsig,found)
     210if (.not.found) then
     211  write(*,*) "phyetat0: Failed loading <ZSIG>"
     212  call abort
     213else
     214  write(*,*) "phyetat0: <ZSIG> range:", &
     215             minval(zsig), maxval(zsig)
     216endif
     217
     218! ZGAM
     219call get_field("ZGAM",zgam,found)
     220if (.not.found) then
     221  write(*,*) "phyetat0: Failed loading <ZGAM>"
     222  call abort
     223else
     224  write(*,*) "phyetat0: <ZGAM> range:", &
     225             minval(zgam), maxval(zgam)
     226endif
     227
     228! ZTHE
     229call get_field("ZTHE",zthe,found)
     230if (.not.found) then
     231  write(*,*) "phyetat0: Failed loading <ZTHE>"
     232  call abort
     233else
     234  write(*,*) "phyetat0: <ZTHE> range:", &
     235             minval(zthe), maxval(zthe)
     236endif
     237
     238! Surface temperature :
     239call get_field("tsurf",tsurf,found,indextime)
     240if (.not.found) then
     241  write(*,*) "phyetat0: Failed loading <tsurf>"
     242  call abort
     243else
     244  write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
     245             minval(tsurf), maxval(tsurf)
     246endif
     247
     248! Surface emissivity
     249call get_field("emis",emis,found,indextime)
     250if (.not.found) then
     251  write(*,*) "phyetat0: Failed loading <emis>"
     252  call abort
     253else
     254  write(*,*) "phyetat0: Surface emissivity <emis> range:", &
     255             minval(emis), maxval(emis)
     256endif
     257
     258! Cloud fraction (added by BC 2010)
     259call get_field("cloudfrac",cloudfrac,found,indextime)
     260if (.not.found) then
     261  write(*,*) "phyetat0: Failed loading <cloudfrac>"
     262  call abort
     263else
     264  write(*,*) "phyetat0: Cloud fraction <cloudfrac> range:", &
     265             minval(cloudfrac), maxval(cloudfrac)
     266endif
     267
     268! Total cloud fraction (added by BC 2010)
     269call get_field("totcloudfrac",totcloudfrac,found,indextime)
     270if (.not.found) then
     271  write(*,*) "phyetat0: Failed loading <totcloudfrac>"
     272  call abort
     273else
     274  write(*,*) "phyetat0: Total cloud fraction <totcloudfrac> range:", &
     275             minval(totcloudfrac), maxval(totcloudfrac)
     276endif
     277
     278! Height of oceanic ice (added by BC 2010)
     279call get_field("hice",hice,found,indextime)
     280if (.not.found) then
     281  write(*,*) "phyetat0: Failed loading <hice>"
     282  call abort
     283else
     284  write(*,*) "phyetat0: Height of oceanic ice <hice> range:", &
     285             minval(hice), maxval(hice)
     286endif
     287
     288! pbl wind variance
     289call get_field("q2",q2,found,indextime)
     290if (.not.found) then
     291  write(*,*) "phyetat0: Failed loading <q2>"
     292  call abort
     293else
     294  write(*,*) "phyetat0: PBL wind variance <q2> range:", &
     295             minval(q2), maxval(q2)
     296endif
     297
     298! tracer on surface
     299if (nq.ge.1) then
     300  do iq=1,nq
     301    txt=tname(iq)
     302    if (txt.eq."h2o_vap") then
     303      ! There is no surface tracer for h2o_vap;
     304      ! "h2o_ice" should be loaded instead
     305      txt="h2o_ice"
     306      write(*,*) 'phyetat0: loading surface tracer', &
     307                           ' h2o_ice instead of h2o_vap'
     308    endif
     309    call get_field(txt,qsurf(:,iq),found,indextime)
     310    if (.not.found) then
     311      write(*,*) "phyetat0: Failed loading <",trim(txt),">"
     312      write(*,*) "         ",trim(txt)," is set to zero"
     313    else
     314      write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
     315                 minval(qsurf(:,iq)), maxval(qsurf(:,iq))
     316    endif
     317  enddo
     318endif ! of if (nq.ge.1)
    673319
    674320! Call to soil_settings, in order to read soil temperatures,
    675321! as well as thermal inertia and volumetric heat capacity
    676322
    677       PRINT*,'SOIL INIT'
    678       call soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
    679 c
    680 c Fermer le fichier:
    681 c
    682       ierr = NF_CLOSE(nid)
    683 c
    684 
    685       RETURN
    686       END
     323call soil_settings(nid_start,ngrid,nsoil,tsurf,tsoil,indextime)
     324
     325!
     326! close file:
     327!
     328call close_startphy
     329
     330END SUBROUTINE phyetat0

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

    r1214 r1216  
    1       subroutine physdem1(ngrid,filename,lonfi,latfi,nsoil,nq,
    2      .                   phystep,day_ini,
    3      .                   time,tsurf,tsoil,emis,q2,qsurf,
    4      .                   airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe,
    5      .                   cloudfrac,totcloudfrac,hice,nametrac)
    6 
    7       use surfdat_h
    8       use comsoil_h
    9       USE comgeomfi_h
    10  
    11       implicit none
    12 
    13 !==================================================================
    14 !     
    15 !     Purpose
    16 !     -------
    17 !     Create physics (re-)start data file "restartfi.nc"
    18 !     
    19 !     Called by
    20 !     ---------
    21 !     rcm1d.F90
    22 !     physiq.F90
    23 !     newstart.F
    24 !     
    25 !     Calls
    26 !     -----
    27 !     none
    28 !
    29 !==================================================================
    30 
    31 #include "dimensions.h"
    32 #include "paramet.h"
    33 #include "comvert.h"
    34 #include "comgeom2.h"
    35 #include "control.h"
    36 #include "comdissnew.h"
    37 #include "logic.h"
    38 #include "ener.h"
    39 #include "netcdf.inc"
    40 #include "dimphys.h"
    41 #include"callkeys.h"
    42 
    43       INTEGER :: ngrid
    44 
    45       INTEGER nid,iq
    46       INTEGER, parameter :: ivap=1
    47       REAL, parameter :: qsolmax= 150.0
    48       character (len=*) :: filename
    49       character (len=7) :: str7
    50 
    51       REAL day_ini
    52       INTEGER nsoil,nq
    53       integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid
    54 
    55       REAL phystep,time
    56       REAL latfi(ngrid), lonfi(ngrid)
    57 !      REAL champhys(ngrid)
    58       REAL tsurf(ngrid)
    59       INTEGER length
    60       PARAMETER (length=100)
    61       REAL tab_cntrl(length)
    62 
    63 !     added by BC
    64       REAL hice(ngrid),cloudfrac(ngrid,nlayermx)
    65       REAL totcloudfrac(ngrid)
    66 
    67 !     AS: name of tracers added as an argument to avoid using nqmx in commons (i.e. adv trac)
    68 !       nota: physdem1 is used both in physiq and newstart.
    69       character*20 nametrac(nq)   ! name of the tracer
    70 
    71 
    72 !      EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
    73 !      EXTERNAL exner_hyb , SSUM
    74 
    75 #include "serre.h"
    76 #include "fxyprim.h"
     1module phyredem
     2
     3implicit none
     4
     5contains
     6
     7subroutine physdem0(filename,lonfi,latfi,nsoil,ngrid,nlay,nq, &
     8                         phystep,day_ini,time,airefi, &
     9                         alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe)
     10! create physics restart file and write time-independent variables
     11  use comsoil_h, only: inertiedat, volcapa, mlayer
     12  use comgeomfi_h, only: area
     13  use surfdat_h, only: albedodat, zmea, zstd, zsig, zgam, zthe, &
     14                       albedice, emisice, emissiv, &
     15                       iceradius, dtemisice, phisfi
     16  use iostart, only : open_restartphy, close_restartphy, &
     17                      put_var, put_field, length
     18  use mod_grid_phy_lmdz, only : klon_glo
     19
     20  implicit none
    7721#include "planete.h"
    7822#include "comcstfi.h"
    79 
    80       real tsoil(ngrid,nsoil),emis(ngrid)
    81       real q2(ngrid, llm+1),qsurf(ngrid,nq)
    82       real airefi(ngrid)
    83       real alb(ngrid),ith(ngrid,nsoil)
    84       real pzmea(ngrid),pzstd(ngrid)
    85       real pzsig(ngrid),pzgam(ngrid),pzthe(ngrid)
    86       integer ig
    87 
    88 ! flag which identifies if we are using old tracer names (qsurf01,...)
    89       logical :: oldtracernames=.false.
    90       integer :: count
    91       character(len=30) :: txt ! to store some text
    92 ! indexes of water vapour & water ice tracers (if any):
    93       integer :: i_h2o_vap=0
    94       integer :: i_h2o_ice=0
    95 c-----------------------------------------------------------------------
    96 
    97       ! copy airefi(:) to area(:)
    98       CALL SCOPY(ngrid,airefi,1,area,1)
    99       ! note: area() is defined in comgeomfi.h
    100 
    101       DO ig=1,ngrid
    102          albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat.h
    103          zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat.h
    104          zstd(ig)=pzstd(ig) ! note: zstd() is defined in surfdat.h
    105          zsig(ig)=pzsig(ig) ! note: zsig() is defined in surfdat.h
    106          zgam(ig)=pzgam(ig) ! note: zgam() is defined in surfdat.h
    107          zthe(ig)=pzthe(ig) ! note: zthe() is defined in surfdat.h
    108       ENDDO
    109 
    110       inertiedat(:,:)=ith(:,:) ! note inertiedat() is defined in comsoil.h
    111 c
    112 c  things to store in the physics start file:
    113 c
    114       ierr = NF_CREATE(adjustl(filename),NF_CLOBBER, nid)
    115       IF (ierr.NE.NF_NOERR) THEN
    116         WRITE(6,*)'physdem1: Problem creating file ',adjustl(filename)
    117         write(6,*) NF_STRERROR(ierr)
    118         CALL ABORT
    119       ENDIF
    120 c
    121       print*,'we got this far'
    122 
    123       ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 18,
    124      .                       "Physics start file")
    125 c
    126       ierr = NF_DEF_DIM (nid,"index",length,idim1)
    127       if (ierr.ne.NF_NOERR) then
    128         WRITE(6,*)'physdem1: Problem defining index dimension'
    129         write(6,*) NF_STRERROR(ierr)
    130         call abort
    131       endif
    132 c
    133       ierr = NF_DEF_DIM (nid,"physical_points",ngrid,idim2)
    134       if (ierr.ne.NF_NOERR) then
    135         WRITE(6,*)'physdem1: Problem defining physical_points dimension'
    136         write(6,*) NF_STRERROR(ierr)
    137         call abort
    138       endif
    139 c
    140       ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoil,idim3)
    141       if (ierr.ne.NF_NOERR) then
    142       WRITE(6,*)'physdem1: Problem defining subsurface_layers dimension'
    143         write(6,*) NF_STRERROR(ierr)
    144         call abort
    145       endif
    146 c
    147       ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
    148       if (ierr.ne.NF_NOERR) then
    149         WRITE(6,*)'physdem1: Problem defining nlayer+1 dimension'
    150         write(6,*) NF_STRERROR(ierr)
    151         call abort
    152       endif
    153 c
    154       ierr = NF_DEF_DIM (nid,"number_of_advected_fields",nq,idim5)
    155       if (ierr.ne.NF_NOERR) then
    156         WRITE(6,*)'physdem1: Problem defining advected fields dimension'
    157         WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
    158         write(6,*) NF_STRERROR(ierr)
    159       endif
    160 
    161       ierr = NF_DEF_DIM (nid,"nlayer",llm,idim6)
    162       if (ierr.ne.NF_NOERR) then
    163         WRITE(6,*)'physdem1: Problem defining nlayer dimension'
    164         write(6,*) NF_STRERROR(ierr)
    165         call abort
    166       endif
    167 
    168       ierr = NF_ENDDEF(nid) ! exit NetCDF define mode
    169 
    170 c clear tab_cntrl(:) array
    171       DO ierr = 1, length
    172          tab_cntrl(ierr) = 0.0
    173       ENDDO
    174 
    175       write(*,*) "physdem1: ngrid: ",ngrid
    176 
    177 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    178 c Fill control array tab_cntrl(:) with paramleters for this run
    179 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    180 c Informations on the physics grid
    181       tab_cntrl(1) = float(ngrid)  ! number of nodes on physics grid
    182       tab_cntrl(2) = float(nlayermx) ! number of atmospheric layers
    183       tab_cntrl(3) = day_ini + int(time)         ! initial day
    184       tab_cntrl(4) = time -int(time)            ! initiale time of day
    185 
    186 c Informations about Mars, used by dynamics and physics
    187       tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
    188       tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
    189       tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
    190       tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
    191       tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
    192       tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
    193 
    194       tab_cntrl(11) = phystep  ! time step in the physics
    195       tab_cntrl(12) = 0.
    196       tab_cntrl(13) = 0.
    197 
    198 c Informations about Mars, only for physics
    199       tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
    200       tab_cntrl(15) = periastr  ! min. star-planet distance (AU)
    201       tab_cntrl(16) = apoastr   ! max. star-planet distance (AU)
    202       tab_cntrl(17) = peri_day  ! date of periastron (sols since N. spring)
    203       tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98
    204 
    205 c Boundary layer and turbulence
    206       tab_cntrl(19) = z0        ! surface roughness (m) ~0.01
    207       tab_cntrl(20) = lmixmin   ! mixing length ~100
    208       tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
    209 
    210 c Optical properties of polar caps and ground emissivity
    211       tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
    212       tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
    213       tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
    214       tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
    215       tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
    216       tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
    217       tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
    218       tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
    219       tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
    220 
    221       tab_cntrl(28) = 0.
    222       tab_cntrl(29) = 0.
    223       tab_cntrl(30) = 0.
     23  character(len=*), intent(in) :: filename
     24  real,intent(in) :: lonfi(ngrid)
     25  real,intent(in) :: latfi(ngrid)
     26  integer,intent(in) :: nsoil
     27  integer,intent(in) :: ngrid
     28  integer,intent(in) :: nlay
     29  integer,intent(in) :: nq
     30  real,intent(in) :: phystep
     31  real,intent(in) :: day_ini
     32  real,intent(in) :: time
     33  real,intent(in) :: airefi(ngrid)
     34  real,intent(in) :: alb(ngrid)
     35  real,intent(in) :: ith(ngrid,nsoil)
     36  real,intent(in) :: pzmea(ngrid)
     37  real,intent(in) :: pzstd(ngrid)
     38  real,intent(in) :: pzsig(ngrid)
     39  real,intent(in) :: pzgam(ngrid)
     40  real,intent(in) :: pzthe(ngrid)
     41 
     42  real :: tab_cntrl(length) ! nb "length=100" defined in iostart module
     43 
     44  ! Create physics start file
     45  call open_restartphy(filename)
     46
     47! tab_cntrl() contains run parameters
     48  tab_cntrl(:)=0 ! initialization
     49!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     50! Fill control array tab_cntrl(:) with paramleters for this run
     51!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     52! Informations on the physics grid
     53  tab_cntrl(1) = float(klon_glo)  ! number of nodes on physics grid
     54  tab_cntrl(2) = float(nlay) ! number of atmospheric layers
     55  tab_cntrl(3) = day_ini + int(time)         ! final day
     56  tab_cntrl(4) = time -int(time)            ! final time of day
     57
     58! Informations about Mars, used by dynamics and physics
     59  tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
     60  tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
     61  tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
     62  tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
     63  tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
     64  tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775
     65
     66  tab_cntrl(11) = phystep  ! time step in the physics
     67  tab_cntrl(12) = 0.
     68  tab_cntrl(13) = 0.
     69
     70! Informations about Mars, only for physics
     71  tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
     72  tab_cntrl(15) = periastr  ! min. star-planet distance (AU)
     73  tab_cntrl(16) = apoastr   ! max. star-planet distance (AU)
     74  tab_cntrl(17) = peri_day  ! date of periastron (sols since N. spring)
     75  tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98
     76
     77! Boundary layer and turbulence
     78  tab_cntrl(19) = z0        ! surface roughness (m) ~0.01
     79  tab_cntrl(20) = lmixmin   ! mixing length ~100
     80  tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8
     81
     82! Optical properties of polar caps and ground emissivity
     83  tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
     84  tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
     85  tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
     86  tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
     87  tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
     88  tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
     89  tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
     90  tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
     91  tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
     92
     93  tab_cntrl(28) = 0.
     94  tab_cntrl(29) = 0.
     95  tab_cntrl(30) = 0.
    22496
    22597! Soil properties:
    226       tab_cntrl(35) = volcapa ! soil volumetric heat capacity
    227      
    228 
    229 !      write(*,*) "physdem1: tab_cntrl():",tab_cntrl
    230      
    231       ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
    232       IF (ierr.NE.NF_NOERR) THEN
    233          PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
    234          CALL abort
    235       ENDIF
    236       ! define variable
    237 #ifdef NC_DOUBLE
    238       ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
    239 #else
    240       ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
    241 #endif
    242       IF (ierr.NE.NF_NOERR) THEN
    243          PRINT*, 'physdem1: Failed to define controle'
    244          CALL abort
    245       ENDIF
    246       ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
    247      .                        "Control parameters")
    248       IF (ierr.NE.NF_NOERR) THEN
    249          PRINT*, 'physdem1: Failed to define controle title attribute'
    250          CALL abort
    251       ENDIF
    252       ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
    253       IF (ierr.NE.NF_NOERR) THEN
    254          PRINT*, 'physdem1: Failed to swich out of NetCDF define mode'
    255          CALL abort
    256       ENDIF
    257       ! write variable
    258 #ifdef NC_DOUBLE
    259       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
    260 #else
    261       ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
    262 #endif
    263       IF (ierr.NE.NF_NOERR) THEN
    264          PRINT*, 'physdem1: Failed to write controle data'
    265          CALL abort
    266       ENDIF
    267 
    268 ! write mid-layer depths mlayer() !known from comsoil.h
    269 
    270       ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
    271       ! define variable
    272 #ifdef NC_DOUBLE
    273       ierr = NF_DEF_VAR (nid,"soildepth",NF_DOUBLE,1,idim3,nvarid)
    274 #else
    275       ierr = NF_DEF_VAR (nid,"soildepth",NF_FLOAT,1,idim3,nvarid)
    276 #endif
    277       ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
    278      .                        "Soil mid-layer depth")
    279       ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
    280       ! write variable
    281 #ifdef NC_DOUBLE
    282       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
    283 #else
    284       ierr = NF_PUT_VAR_REAL (nid,nvarid,mlayer)
    285 #endif
    286 
    287 c
    288 
    289       ierr = NF_REDEF (nid)
    290 #ifdef NC_DOUBLE
    291       ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
    292 #else
    293       ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
    294 #endif
    295       ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
    296      .                        "Longitudes of physics grid")
    297       ierr = NF_ENDDEF(nid)
    298 
    299 #ifdef NC_DOUBLE
    300       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
    301 #else
    302       ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
    303 #endif
    304 
    305 c
    306 
    307       ierr = NF_REDEF (nid)
    308 #ifdef NC_DOUBLE
    309       ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
    310 #else
    311       ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
    312 #endif
    313       ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
    314      .                        "Latitudes of physics grid")
    315       ierr = NF_ENDDEF(nid)
    316 #ifdef NC_DOUBLE
    317       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
    318 #else
    319       ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
    320 #endif
    321 
    322 c
    323 
    324       ierr = NF_REDEF (nid)
    325 #ifdef NC_DOUBLE
    326       ierr = NF_DEF_VAR (nid, "area", NF_DOUBLE, 1, idim2,nvarid)
    327 #else
    328       ierr = NF_DEF_VAR (nid, "area", NF_FLOAT, 1, idim2,nvarid)
    329 #endif
    330       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
    331      .                        "Mesh area")
    332       ierr = NF_ENDDEF(nid)
    333 #ifdef NC_DOUBLE
    334       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
    335 #else
    336       ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
    337 #endif
    338 
    339 c
    340 
    341       ierr = NF_REDEF (nid)
    342 #ifdef NC_DOUBLE
    343       ierr = NF_DEF_VAR (nid, "phisfi", NF_DOUBLE, 1, idim2,nvarid)
    344 #else
    345       ierr = NF_DEF_VAR (nid, "phisfi", NF_FLOAT, 1, idim2,nvarid)
    346 #endif
    347       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
    348      .                        "Geopotential at the surface")
    349       ierr = NF_ENDDEF(nid)
    350 #ifdef NC_DOUBLE
    351       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phisfi)
    352 #else
    353       ierr = NF_PUT_VAR_REAL (nid,nvarid,phisfi)
    354 #endif
    355 
    356 c
    357 
    358       ierr = NF_REDEF (nid)
    359 #ifdef NC_DOUBLE
    360       ierr = NF_DEF_VAR (nid, "albedodat", NF_DOUBLE, 1, idim2,nvarid)
    361 #else
    362       ierr = NF_DEF_VAR (nid, "albedodat", NF_FLOAT, 1, idim2,nvarid)
    363 #endif
    364       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
    365      .                        "Albedo of bare ground")
    366       ierr = NF_ENDDEF(nid)
    367 #ifdef NC_DOUBLE
    368       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedodat)
    369 #else
    370       ierr = NF_PUT_VAR_REAL (nid,nvarid,albedodat)
    371 #endif
    372 
    373 c
    374 c   some data for Francois Lott's programs
    375 c
    376 
    377       ierr = NF_REDEF (nid)
    378 #ifdef NC_DOUBLE
    379       ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
    380 #else
    381       ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
    382 #endif
    383       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
    384      .                        "Relief: mean relief")
    385       ierr = NF_ENDDEF(nid)
    386 #ifdef NC_DOUBLE
    387       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
    388 #else
    389       ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
    390 #endif
    391 c
    392       ierr = NF_REDEF (nid)
    393 #ifdef NC_DOUBLE
    394       ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
    395 #else
    396       ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
    397 #endif
    398       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
    399      .                        "Relief: standard deviation")
    400       ierr = NF_ENDDEF(nid)
    401 #ifdef NC_DOUBLE
    402       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
    403 #else
    404       ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
    405 #endif
    406 c
    407       ierr = NF_REDEF (nid)
    408 #ifdef NC_DOUBLE
    409       ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
    410 #else
    411       ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
    412 #endif
    413       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
    414      .                        "Relief: sigma parameter")
    415       ierr = NF_ENDDEF(nid)
    416 #ifdef NC_DOUBLE
    417       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
    418 #else
    419       ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
    420 #endif
    421 c
    422       ierr = NF_REDEF (nid)
    423 #ifdef NC_DOUBLE
    424       ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
    425 #else
    426       ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
    427 #endif
    428       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
    429      .                        "Relief: gamma parameter")
    430       ierr = NF_ENDDEF(nid)
    431 #ifdef NC_DOUBLE
    432       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
    433 #else
    434       ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
    435 #endif
    436 c
    437       ierr = NF_REDEF (nid)
    438 #ifdef NC_DOUBLE
    439       ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
    440 #else
    441       ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
    442 #endif
    443       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
    444      .                        "Relief: theta parameter")
    445       ierr = NF_ENDDEF(nid)
    446 #ifdef NC_DOUBLE
    447       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
    448 #else
    449       ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
    450 #endif
    451 
    452 c Write the physical fields
    453 
    454 
    455 ! Soil Thermal inertia
    456 
    457       ierr = NF_REDEF (nid)
    458 #ifdef NC_DOUBLE
    459       ierr = NF_DEF_VAR (nid,"inertiedat",NF_DOUBLE,
    460      &                   2,(/idim2,idim3/),nvarid)
    461 #else
    462       ierr = NF_DEF_VAR (nid,"inertiedat",NF_FLOAT,
    463      &                   2,(/idim2,idim3/),nvarid)
    464 #endif
    465       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 20,
    466      .                        "Soil thermal inertia")
    467       ierr = NF_ENDDEF(nid)
    468 #ifdef NC_DOUBLE
    469       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,inertiedat)
    470 #else
    471       ierr = NF_PUT_VAR_REAL (nid,nvarid,inertiedat)
    472 #endif
    473 
    474 
    475 
    476 
    477 !  Surface temperature
    478 
    479       ierr = NF_REDEF (nid)
    480 #ifdef NC_DOUBLE
    481       ierr = NF_DEF_VAR (nid, "tsurf", NF_DOUBLE, 1, idim2,nvarid)
    482 #else
    483       ierr = NF_DEF_VAR (nid, "tsurf", NF_FLOAT, 1, idim2,nvarid)
    484 #endif
    485       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
    486      .                        "Surface temperature")
    487       ierr = NF_ENDDEF(nid)
    488 #ifdef NC_DOUBLE
    489       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsurf)
    490 #else
    491       ierr = NF_PUT_VAR_REAL (nid,nvarid,tsurf)
    492 #endif
    493 
    494 ! Soil temperature
    495 
    496       ierr = NF_REDEF (nid)
    497 #ifdef NC_DOUBLE
    498       ierr = NF_DEF_VAR (nid,"tsoil",NF_DOUBLE,2,(/idim2,idim3/),nvarid)
    499 #else
    500 !      ierr = NF_DEF_VAR (nid, "tsoil", NF_FLOAT, 2, idim2,nvarid)
    501       ierr = NF_DEF_VAR (nid,"tsoil",NF_FLOAT,2,(/idim2,idim3/),nvarid)
    502 #endif
    503       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 16,
    504      .                        "Soil temperature")
    505       ierr = NF_ENDDEF(nid)
    506 #ifdef NC_DOUBLE
    507       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil)
    508 #else
    509       ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil)
    510 #endif
    511 
    512 c emissivity
    513 
    514       ierr = NF_REDEF (nid)
    515 #ifdef NC_DOUBLE
    516       ierr = NF_DEF_VAR (nid, "emis", NF_DOUBLE, 1, idim2,nvarid)
    517 #else
    518       ierr = NF_DEF_VAR (nid, "emis", NF_FLOAT, 1, idim2,nvarid)
    519 #endif
    520       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 18,
    521      .                        "Surface emissivity")
    522       ierr = NF_ENDDEF(nid)
    523 #ifdef NC_DOUBLE
    524       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,emis)
    525 #else
    526       ierr = NF_PUT_VAR_REAL (nid,nvarid,emis)
    527 #endif
    528 
    529 c planetary boundary layer
    530 
    531       ierr = NF_REDEF (nid)
    532 #ifdef NC_DOUBLE
    533       ierr = NF_DEF_VAR (nid,"q2", NF_DOUBLE, 2, (/idim2,idim4/),nvarid)
    534 #else
    535       ierr = NF_DEF_VAR (nid, "q2", NF_FLOAT, 2,(/idim2,idim4/),nvarid)
    536 #endif
    537       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
    538      .                        "pbl wind variance")
    539       ierr = NF_ENDDEF(nid)
    540 #ifdef NC_DOUBLE
    541       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q2)
    542 #else
    543       ierr = NF_PUT_VAR_REAL (nid,nvarid,q2)
    544 #endif
    545       IF (ierr.NE.NF_NOERR) THEN
    546         PRINT*, 'physdem1: Failed to write q2'
    547         CALL abort
    548       ENDIF
    549 
    550 
    551 
    552 
    553 
    554 
    555 
    556 c cloud fraction (added by BC 2010)
    557 
    558       ierr = NF_REDEF (nid)
    559 #ifdef NC_DOUBLE
    560       ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 2,
    561      & (/idim2,idim6/),nvarid)
    562 #else
    563       ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 2,
    564      & (/idim2,idim6/),nvarid)
    565 #endif
    566       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
    567      .                        "Cloud fraction")
    568       ierr = NF_ENDDEF(nid)
    569 #ifdef NC_DOUBLE
    570       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cloudfrac)
    571 #else
    572       ierr = NF_PUT_VAR_REAL (nid,nvarid,cloudfrac)
    573 #endif
    574 
    575 c total cloud fraction (added by BC 2010)
    576 
    577       ierr = NF_REDEF (nid)
    578 #ifdef NC_DOUBLE
    579       ierr = NF_DEF_VAR (nid,"totcloudfrac", NF_DOUBLE, 1, idim2,nvarid)
    580 #else
    581       ierr = NF_DEF_VAR (nid, "totcloudfrac", NF_FLOAT, 1, idim2,nvarid)
    582 #endif
    583       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
    584      .                        "Total fraction")
    585       ierr = NF_ENDDEF(nid)
    586 #ifdef NC_DOUBLE
    587       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,totcloudfrac)
    588 #else
    589       ierr = NF_PUT_VAR_REAL (nid,nvarid,totcloudfrac)
    590 #endif
    591 
    592 c oceanic ice (added by BC 2010)
    593 
    594       ierr = NF_REDEF (nid)
    595 #ifdef NC_DOUBLE
    596       ierr = NF_DEF_VAR (nid, "hice", NF_DOUBLE, 1, idim2,nvarid)
    597 #else
    598       ierr = NF_DEF_VAR (nid, "hice", NF_FLOAT, 1, idim2,nvarid)
    599 #endif
    600       ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
    601      .                        "Height of oceanic ice")
    602       ierr = NF_ENDDEF(nid)
    603 #ifdef NC_DOUBLE
    604       ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,hice)
    605 #else
    606       ierr = NF_PUT_VAR_REAL (nid,nvarid,hice)
    607 #endif
    608 
    609 
    610 
    611 
    612 
    613 
    614 
    615 
    616 
    617 
    618 
    619 c tracers
    620 
    621       if (nq > 0) then
    622 
    623 ! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
    624 !                    qsurf02, ...)
    625       count=0
    626       do iq=1,nq
    627         txt= " "
    628         write(txt,'(a1,i2.2)')'q',iq
    629         if (txt.ne.nametrac(iq)) then ! use tracer names stored in dynamics
    630           ! did not find old tracer name
    631           exit ! might as well stop here
    632         else
    633           ! found old tracer name
    634           count=count+1
    635         endif
    636       enddo
    637       if (count.eq.nq) then
    638         write(*,*) "physdem1:tracers seem to follow old naming ",
    639      &             "convention (qsurf01,qsurf02,...)"
    640 
    641         call abort
    642         !write(*,*) "   => will work for now ... "
    643         !write(*,*) "      but you should run newstart to rename them"
    644         !oldtracernames=.true.
    645         ! Moreover, if computing water cycle with ice, move surface ice
    646         ! back to qsurf(nq)
    647         !IF (iceparty) THEN
    648         !  write(*,*)'physdem1: moving surface water ice to index ',nq
    649         !  qsurf(1:ngrid,nq)=qsurf(1:ngrid,nq-1)
    650         !  qsurf(1:ngrid,nq-1)=0
    651         !ENDIF
    652       endif ! of if (count.eq.nq)
    653 
    654       IF(nq.GE.1) THEN
    655 ! preliminary stuff: look for water vapour & water ice tracers (if any)
    656         if (.not.oldtracernames) then
    657          do iq=1,nq
    658            if (nametrac(iq).eq."h2o_vap") then
    659              i_h2o_vap=iq
    660            endif
    661            if (nametrac(iq).eq."h2o_ice") then
    662              i_h2o_ice=iq
    663            endif
    664          enddo ! of iq=1,nq
    665          ! handle special case of only water vapour tracer (no ice)
    666          if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then
    667           ! then the index of (surface) ice is i_h2o_vap
    668             print*,'water vapour but no water ice, WTF?'
    669             call abort
    670             i_h2o_ice=i_h2o_vap
    671          endif
    672         endif ! of if (.not.oldtracernames)
    673 
    674          DO iq=1,nq
    675            IF (oldtracernames) THEN
    676              txt=" "
    677              write(txt,'(a5,i2.2)')'qsurf',iq
    678            ELSE
    679              txt=nametrac(iq)
    680 
    681 
    682 !     in new version, h2o_vap acts as liquid surface tracer
    683 !     so the section below is now redundant
    684 !     ------------------------------------------------------------------
    685 !             ! Exception: there is no water vapour surface tracer
    686 !             if (txt.eq."h2o_vap") then
    687 !               write(*,*)"physdem1: skipping water vapour tracer"
    688 !               if (i_h2o_ice.eq.i_h2o_vap) then
    689 !               ! then there is no "water ice" tracer; but still
    690 !               ! there is some water ice on the surface
    691 !                 write(*,*)"          writting water ice instead"
    692 !                 txt="h2o_ice"
    693 !               else
    694 !               ! there is a "water ice" tracer which has been / will be
    695 !               ! delt with in due time
    696 !                 cycle
    697 !               endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
    698 !             endif ! of if (txt.eq."h2o_vap")
    699 !     ------------------------------------------------------------------
    700 
    701 
    702 
    703            ENDIF ! of IF (oldtracernames)
    704 
    705            ierr=NF_REDEF(nid)
    706            IF (ierr.NE.NF_NOERR) THEN
    707              PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
    708              CALL abort
    709            ENDIF
    710 #ifdef NC_DOUBLE
    711            ierr=NF_DEF_VAR(nid,txt,NF_DOUBLE,1,idim2,nvarid)
    712 #else
    713            ierr=NF_DEF_VAR(nid,txt,NF_FLOAT,1,idim2,nvarid)
    714 #endif
    715            IF (ierr.NE.NF_NOERR) THEN
    716              PRINT*, 'physdem1: Failed to define ',trim(txt)
    717              CALL abort
    718            ENDIF
    719            ierr=NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
    720      &                        "tracer on surface")
    721            IF (ierr.NE.NF_NOERR) THEN
    722              PRINT*, 'physdem1: Failed to define ',trim(txt),
    723      &               ' title attribute'
    724              CALL abort
    725            ENDIF
    726            ierr=NF_ENDDEF(nid)
    727            IF (ierr.NE.NF_NOERR) THEN
    728              PRINT*, 'physdem1: Failed to swich out of define mode'
    729              CALL abort
    730            ENDIF
    731            
    732 #ifdef NC_DOUBLE
    733             ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,qsurf(1,iq))
    734 #else
    735             ierr=NF_PUT_VAR_REAL (nid,nvarid,qsurf(1,iq))
    736 #endif
    737            IF (ierr.NE.NF_NOERR) THEN
    738              PRINT*, 'physdem1: Failed to write ',trim(txt)
    739              CALL abort
    740            ENDIF
    741          ENDDO ! of DO iq=1,nq
    742       ENDIF ! of IF(nq.GE.1)
    743 
    744       endif ! of if tracer
    745 
    746 c close file
    747       ierr = NF_CLOSE(nid)
    748 
    749       RETURN
    750 
    751       END
     98  tab_cntrl(35) = volcapa ! soil volumetric heat capacity
     99
     100  call put_var("controle","Control parameters",tab_cntrl)
     101 
     102  ! Write the mid-layer depths
     103  call put_var("soildepth","Soil mid-layer depth",mlayer)
     104 
     105  ! Write longitudes
     106  call put_field("longitude","Longitudes of physics grid",lonfi)
     107 
     108  ! Write latitudes
     109  call put_field("latitude","Latitudes of physics grid",latfi)
     110 
     111  ! Write mesh areas
     112  call put_field("area","Mesh area",area)
     113 
     114  ! Write surface geopotential
     115  call put_field("phisfi","Geopotential at the surface",phisfi)
     116 
     117  ! Write surface albedo
     118  call put_field("albedodat","Albedo of bare ground",albedodat)
     119 
     120  ! Subgrid topogaphy variables
     121  call put_field("ZMEA","Relief: mean relief",zmea)
     122  call put_field("ZSTD","Relief: standard deviation",zstd)
     123  call put_field("ZSIG","Relief: sigma parameter",zsig)
     124  call put_field("ZGAM","Relief: gamma parameter",zgam)
     125  call put_field("ZTHE","Relief: theta parameter",zthe)
     126 
     127  ! Soil thermal inertia
     128  call put_field("inertiedat","Soil thermal inertia",inertiedat)
     129 
     130  ! Close file
     131  call close_restartphy
     132 
     133end subroutine physdem0
     134
     135subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
     136                    phystep,time,tsurf,tsoil,emis,q2,qsurf, &
     137                    cloudfrac,totcloudfrac,hice)
     138  ! write time-dependent variable to restart file
     139  use iostart, only : open_restartphy, close_restartphy, &
     140                      put_var, put_field
     141  use infotrac, only: tname
     142  implicit none
     143!======================================================================
     144!#include "temps.h"
     145#include "comcstfi.h"
     146#include "planete.h"
     147!======================================================================
     148  character(len=*),intent(in) :: filename
     149  integer,intent(in) :: nsoil
     150  integer,intent(in) :: ngrid
     151  integer,intent(in) :: nlay
     152  integer,intent(in) :: nq
     153  real,intent(in) :: phystep
     154  real,intent(in) :: time
     155  real,intent(in) :: tsurf(ngrid)
     156  real,intent(in) :: tsoil(ngrid,nsoil)
     157  real,intent(in) :: emis(ngrid)
     158  real,intent(in) :: q2(ngrid,nlay+1)
     159  real,intent(in) :: qsurf(ngrid,nq)
     160  real,intent(in) :: cloudfrac(ngrid,nlay)
     161  real,intent(in) :: totcloudfrac(ngrid)
     162  real,intent(in) :: hice(ngrid)
     163
     164  integer :: iq
     165 
     166  ! Open file
     167  call open_restartphy(filename)
     168
     169  ! First variable to write must be "Time", in order to correctly
     170  ! set time counter in file
     171  !call put_var("Time","Temps de simulation",time)
     172 
     173  ! Surface temperature
     174  call put_field("tsurf","Surface temperature",tsurf)
     175 
     176  ! Soil temperature
     177  call put_field("tsoil","Soil temperature",tsoil)
     178 
     179  ! Emissivity of the surface
     180  call put_field("emis","Surface emissivity",emis)
     181 
     182  ! Planetary Boundary Layer
     183  call put_field("q2","pbl wind variance",q2)
     184 
     185! cloud fraction and sea ice (NB: these should be optional... to be improved)
     186  call put_field("cloudfrac","Cloud fraction",cloudfrac)
     187  call put_field("totcloudfrac","Total fraction",totcloudfrac)
     188  call put_field("hice","Height of oceanic ice",hice)
     189
     190! tracers
     191  if (nq>0) then
     192    do iq=1,nq
     193      call put_field(tname(iq),"tracer on surface",qsurf(:,iq))
     194    enddo
     195  endif ! of if (nq>0)
     196
     197! close file
     198      CALL close_restartphy
     199!$OMP BARRIER
     200
     201end subroutine physdem1
     202
     203end module phyredem

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

    r1194 r1216  
    1010      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
    1111      use watercommon_h, only : RLVTT, Psat_water,epsi
    12       use gases_h
     12      use gases_h, only: gnom, gfrac
    1313      use radcommon_h, only: sigma, eclipse, glat, grav
    1414      use radii_mod, only: h2o_reffrad, co2_reffrad, h2o_cloudrad
    15       use aerosol_mod
    16       use surfdat_h
    17       use comdiurn_h
    18       use comsaison_h
    19       use comsoil_h
    20       USE comgeomfi_h
    21       USE tracer_h
     15      use aerosol_mod, only: iaero_co2, iaero_h2o
     16      use surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe, &
     17                           dryness, watercaptag
     18      use comdiurn_h, only: coslat, sinlat, coslon, sinlon
     19      use comsaison_h, only: mu0, fract, dist_star, declin
     20      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
     21      USE comgeomfi_h, only: long, lati, area, totarea
     22      USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
     23                          alpha_lift, alpha_devil, qextrhor, &
     24                          igcm_h2o_ice, igcm_h2o_vap, igcm_dustbin, &
     25                          igcm_co2_ice
     26      use control_mod, only: ecritphy, iphysiq, nday
     27      use phyredem, only: physdem0, physdem1
    2228
    2329      implicit none
    … …  
    8086!    pt(ngrid,nlayer)      Temperature (K)
    8187!    pq(ngrid,nlayer,nq)   Advected fields
    82 !    pudyn(ngrid,nlayer)    \ 
     88!    pudyn(ngrid,nlayer)    \
    8389!    pvdyn(ngrid,nlayer)     \ Dynamical temporal derivative for the
    8490!    ptdyn(ngrid,nlayer)     / corresponding variables
    … …  
    123129#include "comcstfi.h"
    124130#include "planete.h"
    125 #include "control.h"
     131!#include "control.h"
    126132#include "netcdf.inc"
    127133
    … …  
    621627                            ! need epsi for the wvp definitions in callcorrk.F
    622628
     629         if (ngrid.ne.1) then ! no need to create a restart file in 1d
     630           call physdem0("restartfi.nc",long,lati,nsoilmx,ngrid,nlayer,nq, &
     631                         ptimestep,pday+nday,time_phys,area, &
     632                         albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
     633         endif
     634         
    623635      endif        !  (end of "if firstcall")
    624636
    … …  
    626638! 1.2   Initializations done at every physical timestep:
    627639! ---------------------------------------------------
    628 
    629640!     Initialize various variables
    630641!     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    655666!    Compute variations of g with latitude (oblate case)
    656667!
    657         if (oblate .eq. .false.) then
     668        if (oblate .eqv. .false.) then
    658669           glat(:) = g
    659670        else if (flatten .eq. 0.0 .or. J2 .eq. 0.0 .or. Rmean .eq. 0.0 .or. MassPlanet .eq. 0.0) then
    … …  
    734745               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
    735746               ztim1,ztim2,ztim3,mu0,fract, flatten)
    736            else if(diurnal .eq. .false.) then
     747           else if(diurnal .eqv. .false.) then
    737748
    738749               call mucorr(ngrid,declin,lati,mu0,fract,10000.,rad,flatten)
    … …  
    745756           if(rings_shadow) then
    746757                write(*,*) 'Rings shadow activated'
    747                 if(diurnal .eq. .false.) then ! we need to compute the daily average insolation
     758                if(diurnal .eqv. .false.) then ! we need to compute the daily average insolation
    748759                    pas = 1./nb_hours
    749760                    ptime_day = 0.
    … …  
    784795            endif
    785796
    786 
    787797           if (corrk) then
    788798
    … …  
    15311541
    15321542         if(meanOLR)then
    1533             if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
     1543            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
    15341544               ! to record global radiative balance
    15351545               open(92,file="rad_bal.out",form='formatted',position='append')
    … …  
    16071617
    16081618         if(meanOLR)then
    1609             if((ngrid.gt.1) .or. (mod(icount-1,nint(ecritphy)).eq.0))then
     1619            if((ngrid.gt.1) .or. (mod(icount-1,ecritphy).eq.0))then
    16101620               ! to record global water balance
    16111621               open(98,file="h2o_bal.out",form='formatted',position='append')
    … …  
    16811691            endif
    16821692
    1683             write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
    1684 #ifdef CPP_PARA
    1685 ! for now in parallel we use a different routine to write restartfi.nc
    1686             call phyredem(ngrid,"restartfi.nc",           &
    1687                     ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
    1688                     cloudfrac,totcloudfrac,hice)
    1689 #else
    1690             call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
    1691                     ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
    1692                     area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
    1693                     cloudfrac,totcloudfrac,hice,noms)
    1694 #endif
    1695          endif
     1693            if (ngrid.ne.1) then
     1694              write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin
     1695!#ifdef CPP_PARA
     1696!! for now in parallel we use a different routine to write restartfi.nc
     1697!            call phyredem(ngrid,"restartfi.nc",           &
     1698!                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
     1699!                    cloudfrac,totcloudfrac,hice)
     1700!#else
     1701!            call physdem1(ngrid,"restartfi.nc",long,lati,nsoilmx,nq,            &
     1702!                    ptimestep,pday,ztime_fin,tsurf,tsoil,emis,q2,qsurf_hist, &
     1703!                    area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,      &
     1704!                    cloudfrac,totcloudfrac,hice,noms)
     1705!#endif
     1706              call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq, &
     1707                      ptimestep,ztime_fin, &
     1708                      tsurf,tsoil,emis,q2,qsurf_hist, &
     1709                      cloudfrac,totcloudfrac,hice)
     1710            endif
     1711         
     1712         endif ! of if (lastcall)
    16961713
    16971714

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

    r1150 r1216  
    22
    33! to use  'getin'
    4       use ioipsl_getincom
    5       use surfdat_h
    6       use comdiurn_h
    7       use comsaison_h
    8       use comsoil_h
    9       USE comgeomfi_h
    10 
     4      use ioipsl_getincom, only: getin
     5      use infotrac, only: nqtot, tname
     6      use surfdat_h, only: albedodat, phisfi, dryness, watercaptag,
     7     &                     zmea, zstd, zsig, zgam, zthe,
     8     &                     emissiv, emisice, albedice, iceradius,
     9     &                     dtemisice
     10      use comdiurn_h, only: sinlat, coslat, sinlon, coslon
     11!      use comsaison_h
     12      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
     13      USE comgeomfi_h, only: lati, long, area
     14      use control_mod, only: day_step, ecritphy
     15      use phyredem, only: physdem0,physdem1
     16      use comgeomphy, only: initcomgeomphy
    1117      implicit none
    1218
    … …  
    4147#include "comcstfi.h"
    4248#include "planete.h"
    43 #include "control.h"
     49!#include "control.h"
    4450#include "comvert.h"
    4551#include "netcdf.inc"
    … …  
    6167      REAL play(nlayermx)   ! Pressure at the middle of the layers (Pa)
    6268      REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
    63       REAL psurf,tsurf     
     69      REAL psurf,tsurf(1)     
    6470      REAL u(nlayermx),v(nlayermx)  ! zonal, meridional wind
    6571      REAL gru,grv   ! prescribed "geostrophic" background wind
    6672      REAL temp(nlayermx)   ! temperature at the middle of the layers
    67       REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
    68       REAL qsurf(nqmx)      ! tracer surface budget (e.g. kg.m-2)
     73      REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
     74      REAL,ALLOCATABLE :: qsurf(:)      ! tracer surface budget (e.g. kg.m-2)
    6975      REAL tsoil(nsoilmx)   ! subsurface soik temperature (K)
    7076!      REAL co2ice           ! co2ice layer (kg.m-2) !not used anymore
    … …  
    7278      integer :: i_h2o_ice=0  ! tracer index of h2o ice
    7379      integer :: i_h2o_vap=0  ! tracer index of h2o vapor
    74       REAL emis             ! surface layer
     80      REAL emis(1)             ! surface layer
    7581      REAL q2(nlayermx+1)   ! Turbulent Kinetic Energy
    7682      REAL zlay(nlayermx)   ! altitude estimee dans les couches (km)
    … …  
    8086      REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
    8187      REAL dpsurf   
    82       REAL dq(nlayermx,nqmx)
    83       REAL dqdyn(nlayermx,nqmx)
     88      REAL,ALLOCATABLE :: dq(:,:)
     89      REAL,ALLOCATABLE :: dqdyn(:,:)
    8490
    8591c   Various intermediate variables
    … …  
    8894      REAL phi(nlayermx),h(nlayermx),s(nlayermx)
    8995      REAL pks, ptif, w(nlayermx)
    90       REAL qtotinit, mqtot(nqmx),qtot
    9196      INTEGER ierr, aslun
    9297      REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
    … …  
    112117
    113118!     added by AS to avoid the use of adv trac common
    114       character*20 nametrac(nqmx)   ! name of the tracer (no need for adv trac common)
    115 
     119      character*20,allocatable :: nametrac(:)   ! name of the tracer (no need for adv trac common)
     120
     121      real :: latitude, longitude
    116122
    117123c=======================================================================
    118124c INITIALISATION
    119125c=======================================================================
     126! initialize "serial/parallel" related stuff
     127      CALL init_phys_lmdz(iim,jjm+1,llm,1,(/(jjm-1)*iim+2/))
     128      call initcomgeomphy
    120129
    121130      !! those are defined in surfdat_h.F90
    … …  
    227236          ! read number of tracers:
    228237          read(90,*,iostat=ierr) nq
     238          nqtot=nq ! set value of nqtot (in infotrac module) as nq
    229239          if (ierr.ne.0) then
    230240            write(*,*) "rcm1d: error reading number of tracers"
    231241            write(*,*) "   (first line of traceur.def) "
    232242            stop
     243          endif
     244          if (nq>0) then
     245            allocate(tname(nq))
     246            allocate(q(nlayermx,nq))
     247            allocate(qsurf(nq))
     248            allocate(dq(nlayermx,nq))
     249            allocate(dqdyn(nlayermx,nq))
    233250          else
    234             ! check that the number of tracers is indeed nqmx
    235             if (nq.ne.nqmx) then
    236               write(*,*) "rcm1d: error, wrong number of tracers:"
    237               write(*,*) "nq=",nq," whereas nqmx=",nqmx
    238               stop
    239             endif
     251            allocate(tname(1)) ! tname(1) is used below, even if nq=0
    240252          endif
    241253       
    242254          do iq=1,nq
    243255          ! minimal version, just read in the tracer names, 1 per line
    244             read(90,*,iostat=ierr) nametrac(iq)
     256            read(90,*,iostat=ierr) tname(iq)
    245257            if (ierr.ne.0) then
    246258              write(*,*) 'rcm1d: error reading tracer names...'
    … …  
    253265         i_h2o_vap=0
    254266         do iq=1,nq
    255            if (nametrac(iq)=="co2_ice") then
     267           if (tname(iq)=="co2_ice") then
    256268             i_co2_ice=iq
    257            elseif (nametrac(iq)=="h2o_ice") then
     269           elseif (tname(iq)=="h2o_ice") then
    258270             i_h2o_ice=iq
    259            elseif (nametrac(iq)=="h2o_vap") then
     271           elseif (tname(iq)=="h2o_vap") then
    260272             i_h2o_vap=iq
    261273           endif
    … …  
    275287
    276288      else
    277        nq=nqmx
     289       nq=nqtot
     290       if (nq>0) then
     291            allocate(tname(nq))
     292            allocate(q(nlayermx,nq))
     293            allocate(qsurf(nq))
     294            allocate(dq(nlayermx,nq))
     295            allocate(dqdyn(nlayermx,nq))
     296        else
     297            allocate(tname(1)) ! tname(1) is used below, even if nq=0
     298        endif
    278299       ! Check that tracer boolean is compliant with number of tracers
    279300       ! -- otherwise there is an error (and more generally we have to be consistent)
    … …  
    290311        do iq=1,nq
    291312          write(str7,'(a1,i2.2)')'q',iq
    292           nametrac(iq)=str7
     313          tname(iq)=str7
    293314        enddo
    294315        ! actually, we'll need at least one "co2_ice" tracer
    295316        ! (for surface CO2 ice)
    296         nametrac(1)="co2_ice"
    297317        i_co2_ice=1
     318        tname(i_co2_ice)="co2_ice"
    298319      endif ! of if (tracer)
    299320
    … …  
    315336      phisfi(1)=0.E+0
    316337     !!! LATITUDE. can be set.
    317       lati(1)=0 ! default value for lati(1)
     338      latitude=0 ! default value for latitude
    318339      PRINT *,'latitude (in degrees) ?'
    319       call getin("latitude",lati(1))
    320       write(*,*) " latitude = ",lati(1)
    321       lati(1)=lati(1)*pi/180.E+0
     340      call getin("latitude",latitude)
     341      write(*,*) " latitude = ",latitude
     342      latitude=latitude*pi/180.E+0
    322343     !!! LONGITUDE. useless in 1D.
    323       long(1)=0.E+0
    324       long(1)=long(1)*pi/180.E+0
     344      longitude=0.E+0
     345      longitude=longitude*pi/180.E+0
    325346
    326347!!! CALL INIFIS
    … …  
    332353!!! - physical frequency: nevermind, in inifis this is a simple print
    333354      cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
    334       CALL inifis(1,llm,day0,daysec,dtphys,
    335      .            lati,long,area,rad,g,r,cpp)
     355      CALL inifis(1,llm,nq,day0,daysec,dtphys,
     356     .            latitude,longitude,area,rad,g,r,cpp)
    336357!!! We check everything is OK.
    337358      PRINT *,"CHECK"
    … …  
    558579c  emissivity / surface co2 ice ( + h2o ice??)
    559580c  -------------------------------------------
    560       emis=emissiv ! default value for emissivity
     581      emis(1)=emissiv ! default value for emissivity
    561582      PRINT *,'Emissivity of bare ground ?'
    562       call getin("emis",emis)
    563       write(*,*) " emis = ",emis
    564       emissiv=emis ! we do this so that condense_co2 sets things to the right
     583      call getin("emis",emis(1))
     584      write(*,*) " emis = ",emis(1)
     585      emissiv=emis(1) ! we do this so that condense_co2 sets things to the right
    565586                   ! value if there is no snow
    566587
    … …  
    573594            ! if we have some CO2 ice on the surface, change emissivity
    574595            if (lati(1).ge.0) then ! northern hemisphere
    575               emis=emisice(1)
     596              emis(1)=emisice(1)
    576597            else ! southern hemisphere
    577               emis=emisice(2)
     598              emis(1)=emisice(2)
    578599            endif
    579600         ENDIF
    … …  
    673694      call profile(nlayer+1,tmp1,tmp2)
    674695
    675       tsurf=tmp2(0)
     696      tsurf(1)=tmp2(0)
    676697      DO ilayer=1,nlayer
    677698        temp(ilayer)=tmp2(ilayer)
    678699      ENDDO
    679700      print*,"check"
    680       PRINT*,"INPUT SURFACE TEMPERATURE",tsurf
     701      PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
    681702      PRINT*,"INPUT TEMPERATURE PROFILE",temp
    682703
    … …  
    698719      DO isoil=1,nsoil
    699720         inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
    700          tsoil(isoil)=tsurf  ! soil temperature
     721         tsoil(isoil)=tsurf(1)  ! soil temperature
    701722      ENDDO
    702723
    … …  
    711732
    712733
    713 c  Ecriture de "startfi"
     734c  Write a "startfi" file
    714735c  --------------------
    715 c  (Ce fichier sera aussitot relu au premier
    716 c   appel de "physiq", mais il est necessaire pour passer
    717 c   les variables purement physiques a "physiq"...
    718 
    719       call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
    720      &              dtphys,float(day0),
    721      &              time,tsurf,tsoil,emis,q2,qsurf,
    722      &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
    723      &              cloudfrac,totcloudfrac,hice,nametrac)
     736c  This file will be read during the first call to "physiq".
     737c  It is needed to transfert physics variables to "physiq"...
     738
     739      call physdem0("startfi.nc",long,lati,nsoilmx,1,llm,nq,
     740     &              dtphys,real(day0),time,area,
     741     &              albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe)
     742      call physdem1("startfi.nc",nsoilmx,1,llm,nq,
     743     &                dtphys,time,
     744     &                tsurf,tsoil,emis,q2,qsurf,
     745     &                cloudfrac,totcloudfrac,hice)
     746
     747!      call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
     748!     &              dtphys,float(day0),
     749!     &              time,tsurf,tsoil,emis,q2,qsurf,
     750!     &              area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
     751!     &              cloudfrac,totcloudfrac,hice,nametrac)
    724752
    725753c=======================================================================
    … …  
    775803
    776804      CALL physiq (1,llm,nq,
    777      .     nametrac,
     805     .     tname,
    778806     ,     firstcall,lastcall,
    779807     ,     day,time,dtphys,
    … …  
    858886      ENDDO   ! fin de la boucle temporelle
    859887
     888      write(*,*) "rcm1d: Everything is cool."
     889
    860890c    ========================================================
    861891      end                       !rcm1d
    … …  
    874904
    875905#include "../dyn3d/disvert.F"
     906#include "../dyn3d/abort_gcm.F"

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

    r787 r1216  
    44     &          capcal,fluxgrd)
    55
    6       use comsoil_h
     6      use comsoil_h, only: layer, mlayer, volcapa
    77
    88      implicit none
    … …  
    1717!-----------------------------------------------------------------------
    1818
    19 #include "dimensions.h"
    20 #include "dimphys.h"
     19      include "dimensions.h"
     20      include "dimphys.h"
    2121
    2222
    … …  
    2525!  ---------
    2626!  inputs:
    27       integer ngrid     ! number of (horizontal) grid-points
    28       integer nsoil     ! number of soil layers
    29       logical firstcall ! identifier for initialization call
    30       logical lastcall
    31       real therm_i(ngrid,nsoil) ! thermal inertia
    32       real timestep         ! time step
    33       real tsurf(ngrid)   ! surface temperature
     27      integer,intent(in) :: ngrid       ! number of (horizontal) grid-points
     28      integer,intent(in) :: nsoil       ! number of soil layers
     29      logical,intent(in) :: firstcall ! identifier for initialization call
     30      logical,intent(in) :: lastcall
     31      real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia
     32      real,intent(in) :: timestep           ! time step
     33      real,intent(in) :: tsurf(ngrid)   ! surface temperature
    3434! outputs:
    35       real tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
    36       real capcal(ngrid) ! surface specific heat
    37       real fluxgrd(ngrid) ! surface diffusive heat flux
     35      real,intent(out) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature
     36      real,intent(out) :: capcal(ngrid) ! surface specific heat
     37      real,intent(out) :: fluxgrd(ngrid) ! surface diffusive heat flux
    3838
    3939! local saved variables:
    … …  
    5050
    5151
    52 !      print*,'tsoil=',tsoil
    53 !      print*,'tsurf=',tsurf
    54 
    5552! 0. Initialisations and preprocessing step
    5653      if (firstcall) then
    … …  
    5855      !       and not changed by soil.F
    5956
    60       ALLOCATE(mthermdiff(ngrid,0:nsoilmx-1)) ! mid-layer thermal diffusivity
    61       ALLOCATE(thermdiff(ngrid,nsoilmx-1))    ! inter-layer thermal diffusivity
    62       ALLOCATE(coefq(0:nsoilmx-1))              ! q_{k+1/2} coefficients
    63       ALLOCATE(coefd(ngrid,nsoilmx-1))        ! d_k coefficients
    64       ALLOCATE(alph(ngrid,nsoilmx-1))         ! alpha_k coefficients
    65       ALLOCATE(beta(ngrid,nsoilmx-1))         ! beta_k coefficients
     57      ALLOCATE(mthermdiff(ngrid,0:nsoil-1)) ! mid-layer thermal diffusivity
     58      ALLOCATE(thermdiff(ngrid,nsoil-1))    ! inter-layer thermal diffusivity
     59      ALLOCATE(coefq(0:nsoil-1))              ! q_{k+1/2} coefficients
     60      ALLOCATE(coefd(ngrid,nsoil-1))        ! d_k coefficients
     61      ALLOCATE(alph(ngrid,nsoil-1))         ! alpha_k coefficients
     62      ALLOCATE(beta(ngrid,nsoil-1))         ! beta_k coefficients
    6663
    6764! 0.1 Build mthermdiff(:), the mid-layer thermal diffusivities
    … …  
    187184      enddo
    188185
    189       if (lastcall) then
    190         IF ( ALLOCATED(mthermdiff)) DEALLOCATE(mthermdiff)
    191         IF ( ALLOCATED(thermdiff)) DEALLOCATE(thermdiff)
    192         IF ( ALLOCATED(coefq)) DEALLOCATE(coefq)
    193         IF ( ALLOCATED(coefd)) DEALLOCATE(coefd)
    194         IF ( ALLOCATED(alph)) DEALLOCATE(alph)
    195         IF ( ALLOCATED(beta)) DEALLOCATE(beta)
    196       endif
    197 
    198186      end
    199187

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

    r787 r1216  
    1       subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil)
    2 
    3       use comsoil_h
    4 
     1      subroutine soil_settings(nid,ngrid,nsoil,tsurf,tsoil,indextime)
     2
     3!      use netcdf
     4      use comsoil_h, only: layer, mlayer, inertiedat, volcapa
     5      use iostart, only: inquire_field_ndims, get_var, get_field,
     6     &                   inquire_field, inquire_dimension_length
    57      implicit none
    68
    … …  
    911!
    1012!  Purpose: Read and/or initialise soil depths and properties
     13!
     14! Modifications: Aug.2010 EM : use NetCDF90 to load variables (enables using
     15!                      r4 or r8 restarts independently of having compiled
     16!                      the GCM in r4 or r8)
     17!                June 2013 TN : Possibility to read files with a time axis
     18!
    1119!
    1220!  This subroutine reads from a NetCDF file (opened by the caller)
    … …  
    2634!======================================================================
    2735
    28 #include "dimensions.h"
    29 #include "dimphys.h"
    30 
    31 #include "netcdf.inc"
    3236!======================================================================
    3337!  arguments
    3438!  ---------
    3539!  inputs:
    36       integer nid       ! Input Netcdf file ID
    37       integer ngrid     ! # of horizontal grid points
    38       integer nsoil     ! # of soil layers
    39       integer sta       ! # at which reading starts
    40       real tsurf(ngrid) ! surface temperature
     40      integer,intent(in) :: nid ! Input Netcdf file ID
     41      integer,intent(in) :: ngrid       ! # of horizontal grid points
     42      integer,intent(in) :: nsoil       ! # of soil layers
     43      real,intent(in) :: tsurf(ngrid)   ! surface temperature
     44      integer,intent(in) :: indextime   ! position on time axis
    4145!  output:
    42       real tsoil(ngrid,nsoilmx) ! soil temperature
     46      real,intent(out) :: tsoil(ngrid,nsoil)    ! soil temperature
    4347
    4448!======================================================================
    … …  
    5054      integer ndims     ! # of dimensions of read <inertiedat> data
    5155      integer ig,iloop  ! loop counters
     56     
     57      integer edges(3),corner(3) ! to read a specific time
    5258
    5359      logical :: olddepthdef=.false. ! flag
    … …  
    7076      real,parameter :: default_volcapa=1.e6
    7177
    72       integer, dimension(2) :: sta2d
    73       integer, dimension(2) :: siz2d
    74 
    75 !======================================================================
    76 
    77       !! this is defined in comsoil_h
    78       IF (.not.ALLOCATED(layer))
    79      . ALLOCATE(layer(nsoilmx))
    80       IF (.not.ALLOCATED(mlayer))
    81      . ALLOCATE(mlayer(0:nsoilmx-1))
    82       IF (.not.ALLOCATED(inertiedat))
    83      . ALLOCATE(inertiedat(ngrid,nsoilmx))
    84 
    85       !! this is defined in dimphys.h
    86       sta = cursor
     78      logical :: found,ok
     79     
     80!======================================================================
    8781
    8882! 1. Depth coordinate
    … …  
    9185! 1.1 Start by reading how many layers of soil there are
    9286
    93         ierr=NF_INQ_DIMID(nid,"subsurface_layers",dimid)
    94         if (ierr.ne.NF_NOERR) then
    95          write(*,*)'soil_settings: Problem reading <subsurface_layers>'
    96          call abort
    97         endif
    98 
    99         ierr=NF_INQ_DIMLEN(nid,dimid,dimlen)
    100         if (ierr.ne.NF_NOERR) then
    101          write(*,*)'soil_settings: Problem getting <subsurface_layers>',
    102      &             'length'
    103          call abort
    104         endif
     87        dimlen=inquire_dimension_length("subsurface_layers")
    10588
    10689        if (dimlen.ne.nsoil) then
    … …  
    11295        endif
    11396
    114       sta2d = (/sta,1/)
    115       siz2d = (/ngrid,dimlen/)
    116 
    117 
    11897! 1.2 Find out the # of dimensions <inertiedat> was defined as using
    119 !     (in ye old days, thermal inertia was only given at the "surface")
    120       ! Look for thermal inertia data
    121       ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)
    122       if (ierr.NE.NF_NOERR) then
    123          write(*,*)'soil_settings: Field <inertiedat> not found!'
    124          call abort
    125       endif
    126 
    127       ! Read the # of dimensions <inertidat> was defined as using
    128       ierr=NF_INQ_VARNDIMS(nid,nvarid,ndims)
    129       ! if (ndims.eq.1) then we have the "old 2D-surface" format
     98
     99      ndims=inquire_field_ndims("inertiedat")
    130100
    131101! 1.3 Read depths values or set olddepthdef flag and values
    … …  
    147117        enddo
    148118      else ! Look for depth
    149         ierr=NF_INQ_VARID(nid,"soildepth",nvarid)
    150         if (ierr.ne.NF_NOERR) then
    151           write(*,*)'soil_settings: Field <soildepth> not found!'
    152           call abort
    153         endif
    154119        ! read <depth> coordinate
    155120        if (interpol) then !put values in oldmlayer
    156           if (.not.allocated(oldmlayer)) then
    157              allocate(oldmlayer(dimlen),stat=ierr)
    158           endif
    159 #ifdef NC_DOUBLE
    160            ierr = NF_GET_VAR_DOUBLE(nid,nvarid,oldmlayer)
    161 #else
    162            ierr = NF_GET_VAR_REAL(nid,nvarid,oldmlayer)
    163 #endif
    164           if (ierr.ne.NF_NOERR) then
    165            write(*,*)'soil_settings: Problem while reading <soildepth>'
    166            call abort
     121          call get_var("soildepth",oldmlayer,found)
     122          if (.not.found) then
     123            write(*,*)'soil_settings: Problem while reading <soildepth>'
    167124          endif
    168125        else ! put values in mlayer
    169 #ifdef NC_DOUBLE
    170            ierr = NF_GET_VAR_DOUBLE(nid,nvarid,mlayer)
    171 #else
    172            ierr = NF_GET_VAR_REAL(nid,nvarid,mlayer)
    173 #endif
    174           if (ierr.ne.NF_NOERR) then
    175            write(*,*)'soil_settings: Problem while reading <soildepth>'
    176            call abort
     126          call get_var("soildepth",mlayer,found)
     127          if (.not.found) then
     128            write(*,*)'soil_settings: Problem while reading <soildepth>'
    177129          endif
    178130        endif !of if (interpol)
    … …  
    183135      ! default mlayer distribution, following a power law:
    184136      !  mlayer(k)=lay1*alpha**(k-1/2)
    185         lay1=2.e-4  !mars case (nsoilmx=18)
    186 !        lay1=3.e-2  !earth case (nsoilmx=13)
     137        lay1=2.e-4
    187138        alpha=2
    188139        do iloop=0,nsoil-1
    … …  
    200151      enddo
    201152
    202 ! 2. Volumetric heat capacity (note: it is declared in comsoil.h)
     153! 2. Volumetric heat capacity (note: it is declared in comsoil_h)
    203154! ---------------------------
    204155! "volcapa" is (so far) 0D and written in "controle" table of startfi file
    … …  
    214165        volcapa=default_volcapa
    215166      endif
    216 ! Look for it
    217 !      ierr=NF_INQ_VARID(nid,"volcapa",nvarid)
    218 !      if (ierr.NE.NF_NOERR) then
    219 !         write(*,*)'soil_settings: Field <volcapa> not found!'
    220 !         write(*,*)'Initializing Volumetric heat capacity to ',
    221 !     &             default_volcapa
    222 !         volcapa=default_volcapa
    223 !      else
    224 !#ifdef NC_DOUBLE
    225 !       ierr = NF_GET_VAR_DOUBLE(nid,nvarid,volcapa)
    226 !#else
    227 !       ierr = NF_GET_VAR_REAL(nid,nvarid,volcapa)
    228 !#endif
    229 !        if (ierr.ne.NF_NOERR) then
    230 !         write(*,*)'soil_settings: Problem while reading <volcapa>'
    231 !         call abort
    232 !       endif
    233 !      endif
    234 
    235 ! 3. Thermal inertia (note: it is declared in comsoil.h)
     167
     168! 3. Thermal inertia (note: it is declared in comsoil_h)
    236169! ------------------
    237170
    238171! 3.1 Look (again) for thermal inertia data (to reset nvarid)
    239       ierr=NF_INQ_VARID(nid, "inertiedat", nvarid)
    240       if (ierr.NE.NF_NOERR) then
    241          write(*,*)'soil_settings: Field <inertiedat> not found!'
    242          call abort
    243       endif
    244172
    245173! 3.2 Knowing the # of dimensions <inertidat> was defined as using,
    … …  
    251179       ! Read Surface thermal inertia
    252180       allocate(surfinertia(ngrid))
    253 #ifdef NC_DOUBLE
    254        ierr = NF_GET_VARA_DOUBLE(nid, nvarid, sta, ngrid, surfinertia)
    255 #else
    256        ierr = NF_GET_VARA_REAL(nid, nvarid, sta, ngrid, surfinertia)
    257 #endif
    258         if (ierr.NE.NF_NOERR) then
    259          write(*,*)'soil_settings: Problem while reading <inertiedat>'
     181       call get_field("inertiedat",surfinertia,found)
     182       if (.not.found) then
     183         write(*,*) "soil_settings: Failed loading <inertiedat>"
    260184         call abort
    261         endif
     185       endif
    262186       
    263187       write(*,*)' => Building soil thermal inertia (using reference sur
    … …  
    277201            stop
    278202           endif
    279          endif
    280 #ifdef NC_DOUBLE
    281         ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldinertiedat)
    282 #else
    283         ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldinertiedat)
    284 #endif
    285         if (ierr.NE.NF_NOERR) then
    286          write(*,*)'soil_settings: Problem while reading <inertiedat>'
    287          call abort
     203         endif ! of if (.not.allocated(oldinertiedat))
     204        call get_field("inertiedat",oldinertiedat,found)
     205        if (.not.found) then
     206          write(*,*) "soil_settings: Failed loading <inertiedat>"
     207          call abort
    288208        endif
    289209       else ! put values in therm_i
    290 #ifdef NC_DOUBLE
    291         ierr = NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,inertiedat)
    292 #else
    293         ierr = NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,inertiedat)
    294 #endif
    295         if (ierr.NE.NF_NOERR) then
    296          write(*,*)'soil_settings: Problem while reading <inertiedat>'
    297          call abort
    298         endif
     210         call get_field("inertiedat",inertiedat,found)
     211         if (.not.found) then
     212           write(*,*) "soil_settings: Failed loading <inertiedat>"
     213           call abort
     214         endif
     215!        endif
    299216       endif ! of if (interpol)
    300217      endif ! of if (ndims.eq.1)
    … …  
    303220! -------------------------
    304221
    305       ierr=NF_INQ_VARID(nid,"tsoil",nvarid)
    306       if (ierr.ne.NF_NOERR) then
     222!      ierr=nf90_inq_varid(nid,"tsoil",nvarid)
     223      ok=inquire_field("tsoil")
     224!      if (ierr.ne.nf90_noerr) then
     225      if (.not.ok) then
    307226        write(*,*)'soil_settings: Field <tsoil> not found!'
    308227        write(*,*)' => Building <tsoil> from surface values <tsurf>'
    … …  
    320239           endif
    321240         endif
    322 #ifdef NC_DOUBLE
    323         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,oldtsoil)
    324 #else
    325         ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,oldtsoil)
    326 #endif
    327         if (ierr.ne.NF_NOERR) then
    328          write(*,*)'soil_settings: Problem while reading <tsoil>'
    329          call abort
    330         endif
     241         call get_field("tsoil",oldtsoil,found)
     242         if (.not.found) then
     243           write(*,*) "soil_settings: Failed loading <tsoil>"
     244           call abort
     245         endif
    331246       else ! put values in tsoil
    332 #ifdef NC_DOUBLE
    333         ierr=NF_GET_VARA_DOUBLE(nid,nvarid,sta2d,siz2d,tsoil)
    334 #else
    335         ierr=NF_GET_VARA_REAL(nid,nvarid,sta2d,siz2d,tsoil)
    336 #endif
    337         if (ierr.ne.NF_NOERR) then
    338          write(*,*)'soil_settings: Problem while reading <tsoil>'
    339          call abort
    340         endif
     247         call get_field("tsoil",tsoil,found,timeindex=indextime)
     248         if (.not.found) then
     249           write(*,*) "soil_settings: Failed loading <tsoil>"
     250           call abort
     251         endif
    341252       endif ! of if (interpol)
    342       endif
     253      endif! of if (.not.ok)
    343254
    344255! 5. If necessary, interpolate soil temperatures and thermal inertias

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

    r993 r1216  
    1919c=======================================================================
    2020
     21      use infotrac, only: iniadvtrac, nqtot, tname
    2122      USE comsoil_h
    22       USE comgeomfi_h
    23 
     23      USE comgeomfi_h, ONLY: lati, long, area
     24!      use control_mod
     25      use comgeomphy, only: initcomgeomphy
    2426      implicit none
    2527
    … …  
    3234#include "logic.h"
    3335#include "temps.h"
    34 #include "control.h"
     36!#include "control.h"
    3537#include "ener.h"
    3638
    3739#include "dimphys.h"
    3840#include "planete.h"
    39 #include"advtrac.h"
     41!#include"advtrac.h"
    4042#include "netcdf.inc"
    4143
    … …  
    4850      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
    4951      REAL teta(ip1jmp1,llm)                    ! temperature potentielle
    50       REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
     52      REAL,ALLOCATABLE :: q(:,:,:)   ! champs advectes
    5153      REAL pks(ip1jmp1)                      ! exner (f pour filtre)
    5254      REAL pk(ip1jmp1,llm)
    … …  
    6365      REAL tsoil(ngridmx,nsoilmx) ! Soil temperature
    6466      REAL co2ice(ngridmx)      ! CO2 ice layer
    65       REAL q2(ngridmx,nlayermx+1),qsurf(ngridmx,nqmx)
     67      REAL q2(ngridmx,nlayermx+1)
     68      REAL,ALLOCATABLE :: qsurf(:,:)
    6669      REAL emis(ngridmx)
    6770      INTEGER start,length
    … …  
    8386      REAL ithS(ip1jmp1,nsoilmx) ! Soil Thermal Inertia
    8487      REAL co2iceS(ip1jmp1)
    85       REAL q2S(ip1jmp1,llm+1),qsurfS(ip1jmp1,nqmx)
     88      REAL q2S(ip1jmp1,llm+1)
     89      REAL,ALLOCATABLE :: qsurfS(:,:)
    8690      REAL emisS(ip1jmp1)
    8791
    … …  
    122126c-----------------------------------------------------------------------
    123127
     128      CALL defrun_new(99, .TRUE. )
    124129      grireg   = .TRUE.
    125130      cursor = 1 ! added by AS in dimphys.
     131
     132! initialize "serial/parallel" related stuff
     133      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
     134      call initcomgeomphy
    126135
    127136      ! ALLOCATE ARRAYS IN comgeomfi_h (usually done in inifis)
    … …  
    134143c Lecture des donnees
    135144c=======================================================================
    136 ! Load tracer names:
    137       call iniadvtrac(nq,numvanle)
    138 ! Initialize global tracer indexes (stored in tracer.h)
    139 ! ... this has to be done before phyetat0
    140       call initracer(ngridmx,nqmx,tnom)
     145! Load tracer number and names:
     146      call iniadvtrac(nqtot,numvanle)
     147
     148! allocate arrays:
     149      allocate(q(ip1jmp1,llm,nqtot))
     150      allocate(qsurf(ngridmx,nqtot))
     151      allocate(qsurfS(ip1jmp1,nqtot))
    141152
    142153      fichnom = 'start.nc'
    143       CALL dynetat0(fichnom,nqmx,vcov,ucov,teta,q,masse,
     154      CALL dynetat0(fichnom,nqtot,vcov,ucov,teta,q,masse,
    144155     .       ps,phis,timedyn)
    145156
    … …  
    173184      Lmodif=0
    174185
    175       CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqmx,day_ini_fi,
     186      CALL phyetat0 (ngridmx,fichnom,0,Lmodif,nsoilmx,nqtot,day_ini_fi,
    176187     .      timefi,
    177188     .      tsurf,tsoil,emis,q2,qsurf,
    … …  
    286297      call gr_fi_dyn(1,ngridmx,iip1,jjp1,emis,emisS)
    287298      call gr_fi_dyn(llm+1,ngridmx,iip1,jjp1,q2,q2S)
    288       call gr_fi_dyn(nqmx,ngridmx,iip1,jjp1,qsurf,qsurfS)
     299      call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
    289300      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,cloudfrac,cloudfracS)
    290301      call gr_fi_dyn(1,ngridmx,iip1,jjp1,hice,hiceS)
    … …  
    390401
    391402c-----------------------------------------------------------------------
    392 c Ecriture du champs  q  ( q[1,nqmx] )
    393 c-----------------------------------------------------------------------
    394       do iq=1,nqmx
    395         call write_archive(nid,ntime,tnom(iq),'tracer','kg/kg',
     403c Ecriture du champs  q  ( q[1,nqtot] )
     404c-----------------------------------------------------------------------
     405      do iq=1,nqtot
     406        call write_archive(nid,ntime,tname(iq),'tracer','kg/kg',
    396407     &         3,q(1,1,iq))
    397408      end do
    398409c-----------------------------------------------------------------------
    399 c Ecriture du champs  qsurf  ( qsurf[1,nqmx] )
    400 c-----------------------------------------------------------------------
    401       do iq=1,nqmx
    402         txt=trim(tnom(iq))//"_surf"
     410c Ecriture du champs  qsurf  ( qsurf[1,nqtot] )
     411c-----------------------------------------------------------------------
     412      do iq=1,nqtot
     413        txt=trim(tname(iq))//"_surf"
    403414        call write_archive(nid,ntime,txt,'Tracer on surface',
    404415     &  'kg.m-2',2,qsurfS(1,iq))
    … …  
    448459      ierr=NF_CLOSE(nid)
    449460
     461      write(*,*) "start2archive: All is well that ends well."
     462
    450463      end

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

    r1174 r1216  
    6565C             ---------
    6666C
    67       INTEGER kgrid
    68       REAL ptim1,ptim2,ptim3, pflat
    69       REAL psilon(kgrid),pcolon(kgrid),pmu0(kgrid),pfract(kgrid)
    70       REAL psilat(kgrid), pcolat(kgrid)
     67      INTEGER,INTENT(IN) :: kgrid
     68      REAL,INTENT(IN) :: ptim1,ptim2,ptim3, pflat
     69      REAL,INTENT(IN) :: psilon(kgrid),pcolon(kgrid)
     70      REAL,INTENT(IN) :: psilat(kgrid), pcolat(kgrid)
     71      REAL,INTENT(OUT) :: pmu0(kgrid),pfract(kgrid)
    7172C
    7273      INTEGER jl
    … …  
    100101        ztim3=pcolat(jl)*ptim3
    101102        pmu0(jl)=ztim1+ztim2*pcolon(jl)+ztim3*psilon(jl)
    102         pmu0(jl)=pmu0(jl)/SQRT(pcolat(jl)**2 + (rap**2) * (psilat(jl)**2))
     103        pmu0(jl)=pmu0(jl)/SQRT(pcolat(jl)**2+(rap**2)*(psilat(jl)**2))
    103104
    104105      ENDDO

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

    r837 r1216  
    11      SUBROUTINE surfini(ngrid,nq,qsurf,psolaralb)
    22
    3       USE surfdat_h
    4       USE tracer_h
     3      USE surfdat_h, only: albedodat, albedice
     4      USE tracer_h, only: igcm_co2_ice
     5      use comgeomfi_h, only: lati
     6      use planetwide_mod, only: planetwide_maxval, planetwide_minval
    57
    68      IMPLICIT NONE
    … …  
    1719#include "callkeys.h"
    1820c
    19       INTEGER ngrid,ig,icap
    20       INTEGER nq
    21       REAL  piceco2(ngrid),psolaralb(ngrid)
    22       REAL qsurf(ngrid,nq) !tracer on surface (kg/m2)
     21      INTEGER,INTENT(IN) :: ngrid
     22      INTEGER,INTENT(IN) :: nq
     23      REAL,INTENT(OUT) :: psolaralb(ngrid)
     24      REAL,INTENT(IN) :: qsurf(ngrid,nq) !tracer on surface (kg/m2)
    2325
    24       EXTERNAL ISMIN,ISMAX
    25       INTEGER ISMIN,ISMAX
     26      INTEGER :: ig,icap
     27      REAL :: min_albedo,max_albedo
    2628c
    2729c=======================================================================
    2830
    29 c
    30 c     calcul de piceco2 (kg/m2) a l'etat initial
    31 c     ------------------------------------------
    3231
    3332      DO ig=1,ngrid
    3433         psolaralb(ig)=albedodat(ig)
    35 !         psolaralb(ig,2)=albedodat(ig)
    3634      ENDDO
    3735
    38       PRINT*,'surfini: minimum des donnees albedo',
    39      s     albedodat(ISMIN(ngrid,albedodat,1))
    40       PRINT*,'surfini: maximum des donnees albedo',
    41      s     albedodat(ISMAX(ngrid,albedodat,1))
    42 
    43 c      calcul de psolaralb
    44 c      -------------------
    45 !      DO ig=1,ngrid
    46 c        IF (water) THEN
    47 c          if (qsurf(ig,nq).gt.0.005) then
    48 c             psolaralb(ig,1) = 0.4
    49 c             psolaralb(ig,2) = 0.4
    50 c           endif
    51 c         ENDIF
    52 !      ENDDO ! of DO ig=1,ngrid     
    53 c IF there is more than 5 pr. um of h2o ice but no C02 ice, surface albedo is set to 0.4.
     36      call planetwide_minval(albedodat,min_albedo)
     37      call planetwide_maxval(albedodat,max_albedo)
     38      write(*,*) 'surfini: minimum corrected albedo',min_albedo
     39      write(*,*) 'surfini: maximum corrected albedo',max_albedo
    5440
    5541      if (igcm_co2_ice.ne.0) then
    … …  
    5743        DO ig=1,ngrid
    5844          IF (qsurf(ig,igcm_co2_ice) .GT. 0.) THEN
    59              !!!! no good in parallel
    60              IF(ig.GT.ngrid/2+1) THEN
    61                 icap=2
     45             IF(lati(ig).LT.0.) THEN
     46                icap=2 ! Southern hemisphere
    6247             ELSE
    63                 icap=1
     48                icap=1 ! Northern hemisphere
    6449             ENDIF
    6550             psolaralb(ig) = albedice(icap)
    66 !             psolaralb(ig,2) =  albedice(icap)
    6751          END IF
    6852        ENDDO ! of DO ig=1,ngrid     
    … …  
    7256      endif
    7357
    74       PRINT*,'surfini: minimum des donnees albedo',
    75      s     psolaralb(ISMIN(ngrid,psolaralb,1))
    76       PRINT*,'surfini: maximum des donnees albedo',
    77      s     psolaralb(ISMAX(ngrid,psolaralb,1))
     58      call planetwide_minval(psolaralb,min_albedo)
     59      call planetwide_maxval(psolaralb,max_albedo)
     60      write(*,*) 'surfini: minimum corrected albedo',min_albedo
     61      write(*,*) 'surfini: maximum corrected albedo',max_albedo
    7862
    79       RETURN
    8063      END

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

    r787 r1216  
    4444c=======================================================================
    4545! to use  'getin'
    46       use ioipsl_getincom
    47 
    48       USE surfdat_h
    49       use comsoil_h
    50       USE comgeomfi_h
     46      use ioipsl_getincom , only: getin
     47
     48      use surfdat_h, only: albedice, emisice, iceradius, dtemisice,
     49     &                     emissiv
     50      use comsoil_h, only: volcapa
     51      use iostart, only: get_var
     52      use mod_phys_lmdz_para, only: is_parallel
    5153
    5254      implicit none
    5355 
    54 #include "dimensions.h"
    55 #include "dimphys.h"
    5656#include "comcstfi.h"
    5757#include "planete.h"
    … …  
    6363c-----------------------------------------------------------------------
    6464
    65       INTEGER ngrid
    66 
    6765c Arguments
    6866c ---------
    69       INTEGER nid,nvarid,tab0
    70       INTEGER*4 day_ini
    71       INTEGER Lmodif
    72       INTEGER lmax
    73       REAL p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
     67      INTEGER,INTENT(IN) :: ngrid,nid,tab0
     68      INTEGER*4,INTENT(OUT) :: day_ini
     69      INTEGER,INTENT(IN) :: Lmodif
     70      INTEGER,INTENT(OUT) :: lmax
     71      REAL,INTENT(OUT) :: p_rad,p_omeg,p_g,p_cpp,p_mugaz,p_daysec,time
    7472
    7573c Variables
    7674c ---------
    77       INTEGER length
    78       parameter (length = 100)
     75      INTEGER,PARAMETER :: length=100
    7976      REAL tab_cntrl(length) ! array in which are stored the run's parameters
    80       INTEGER  ierr
     77      INTEGER  ierr,nvarid
    8178      INTEGER size
    8279      CHARACTER modif*20
    83 
     80      LOGICAL :: found
     81     
     82      write(*,*)"tabfi: nid=",nid," tab0=",tab0," Lmodif=",Lmodif
     83
     84      IF (nid.eq.0) then
     85c-----------------------------------------------------------------------
     86c  Initialization of various physical constants to defaut values (nid = 0 case)
     87c-----------------------------------------------------------------------
     88      ELSE
    8489c-----------------------------------------------------------------------
    8590c  Initialization of physical constants by reading array tab_cntrl(:)
    … …  
    8893c Read 'controle' array
    8994c
    90       ierr = NF_INQ_VARID (nid, "controle", nvarid)
    91       IF (ierr .NE. NF_NOERR) THEN
    92          PRINT*, "Tabfi: Could not find <controle> data"
    93          CALL abort
    94       ENDIF
    95 #ifdef NC_DOUBLE
    96       ierr = NF_GET_VAR_DOUBLE(nid, nvarid, tab_cntrl)
    97 #else
    98       ierr = NF_GET_VAR_REAL(nid, nvarid, tab_cntrl)
    99 #endif
    100       IF (ierr .NE. NF_NOERR) THEN
    101          PRINT*, "Tabfi: Failed reading <controle> array"
    102          CALL abort
    103       ENDIF
    104 
    105        print*,'tabfi: tab_cntrl',tab_cntrl
     95
     96       call get_var("controle",tab_cntrl,found)
     97       if (.not.found) then
     98         write(*,*)"tabfi: Failed reading <controle> array"
     99         call abort
     100       else
     101         write(*,*)'tabfi: tab_cntrl',tab_cntrl
     102       endif
    106103c
    107104c  Initialization of some physical constants
    108105c informations on physics grid
    109       if(ngrid.ne.tab_cntrl(tab0+1)) then
    110          print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
    111          print*,tab_cntrl(tab0+1),ngrid
    112       endif
     106!      if(ngrid.ne.tab_cntrl(tab0+1)) then
     107!         print*,'tabfi: WARNING !!! tab_cntrl(tab0+1).ne.ngrid'
     108!         print*,tab_cntrl(tab0+1),ngrid
     109!      endif
    113110      lmax = nint(tab_cntrl(tab0+2))
    114111      day_ini = tab_cntrl(tab0+3)
    … …  
    146143c soil properties
    147144      volcapa = tab_cntrl(tab0+35) ! volumetric heat capacity
    148 
    149 
    150 
    151 
    152145c-----------------------------------------------------------------------
    153146c       Save some constants for later use (as routine arguments)
    … …  
    160153      p_rad=rad
    161154
     155      ENDIF    ! end of (nid = 0)
    162156
    163157c-----------------------------------------------------------------------
    … …  
    209203c-----------------------------------------------------------------------
    210204c        Modifications...
     205! NB: Modifying controls should only be done by newstart, and in seq mode
     206      if ((Lmodif.eq.1).and.is_parallel) then
     207        write(*,*) "tabfi: Error modifying tab_control should",
     208     &             " only happen in serial mode (eg: by newstart)"
     209        stop
     210      endif
    211211c-----------------------------------------------------------------------
    212212

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

    r993 r1216  
    3232c=======================================================================
    3333
     34      use comsoil_h, only: nsoilmx
    3435      implicit none
    3536
    3637#include "dimensions.h"
    37 #include "dimphys.h"
    3838#include "paramet.h"
    39 #include "control.h"
    4039#include "comvert.h"
    4140#include "comgeom.h"

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

    r993 r1216  
    3939!
    4040!=================================================================
    41 
    42       USE surfdat_h
    43 #ifdef CPP_PARA
     41      use surfdat_h, only: phisfi
     42      use control_mod, only: ecritphy, day_step, iphysiq
    4443      USE mod_phys_lmdz_para, only : is_parallel, is_mpi_root,
    4544     &                               is_master, gather
    4645      USE mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    47 #endif
    4846      implicit none
    4947
    5048! Commons
    51 #include "dimensions.h"
    52 #include "dimphys.h"
    53 #include "paramet.h"
    54 #include "control.h"
    55 #include "comvert.h"
    56 #include "comgeom.h"
    57 #include "netcdf.inc"
    58 #include "temps.h"
     49      include "dimensions.h"
     50      include "paramet.h"
     51      include "comvert.h"
     52      include "comgeom.h"
     53      include "netcdf.inc"
     54      include "temps.h"
    5955
    6056! Arguments on input:
    … …  
    7369      real*4,save :: date
    7470
    75       REAL phis(ip1jmp1) ! surface geopotential on extended lonxlat grid
     71      REAL phis(ip1jmp1)
    7672
    7773      integer irythme
    … …  
    9692      integer,save :: n_nom_def
    9793      integer :: n
    98       integer,parameter :: n_nom_def_max=99
    99       character(len=20),save :: nom_def(n_nom_def_max)
     94      integer,parameter :: n_nom_def_max=199
     95      character(len=120),save :: nom_def(n_nom_def_max)
    10096      logical,save :: firstcall=.true.
    101 
    102       integer dimvert
     97     
     98#ifndef MESOSCALE
    10399
    104100#ifdef CPP_PARA
    … …  
    113109      real phisfi_glo(klon_glo) ! surface geopotential on global physics grid
    114110#else
    115       logical,parameter :: is_parallel=.false.
    116       logical,parameter :: is_master=.true.
    117       logical,parameter :: is_mpi_root=.true.
    118111      real phisfi_glo(ngrid) ! surface geopotential on global physics grid
    119112#endif
     113
    120114!***************************************************************
    121115!Sortie des variables au rythme voulu
    122116
    123       irythme = int(ecritphy) ! sortie au rythme de ecritphy
     117      irythme = ecritphy ! sortie au rythme de ecritphy
    124118!     irythme = iconser  ! sortie au rythme des variables de controle
    125119!     irythme = iphysiq  ! sortie a tous les pas physique
    … …  
    129123!***************************************************************
    130124
    131 
    132125! At very first call, check if there is a "diagfi.def" to use and read it
    133126! -----------------------------------------------------------------------
    … …  
    135128         firstcall=.false.
    136129
    137 !      Open diagfi.def definition file if there is one:
     130  !      Open diagfi.def definition file if there is one:
    138131         open(99,file="diagfi.def",status='old',form='formatted',
    139132     s   iostat=ierr2)
    … …  
    170163      end if
    171164
    172 
    173165! Initialisation of 'firstnom' and create/open the "diagfi.nc" NetCDF file
    174166! ------------------------------------------------------------------------
    … …  
    185177           stop
    186178         endif
    187 
     179         
    188180         zitau = -1 ! initialize zitau
    189181
    … …  
    240232        ! are undefined; so set them to 1
    241233        iphysiq=1
    242         !JL11 irythme=1
     234        irythme=1
    243235        ! NB:
    244236      endif
    … …  
    546538      endif
    547539
     540#endif
    548541      end
    549 

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

    r965 r1216  
    1010! (yielding the sought time series of the variable)
    1111
    12 #ifdef CPP_PARA
     12! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
     13
     14use comsoil_h, only: nsoilmx
     15use control_mod, only: ecritphy, day_step, iphysiq
    1316use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather
    1417use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    15 #endif
    1618
    1719implicit none
    1820
    1921#include"dimensions.h"
    20 #include"dimphys.h"
    2122#include"paramet.h"
    22 #include"control.h"
    2323#include"netcdf.inc"
    2424
    … …  
    3131integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
    3232real,dimension(ngrid,nsoilmx),intent(in) :: px ! variable
    33 ! Note: nsoilmx is a common parameter set in 'dimphys.h'
     33! Note: nsoilmx is a parameter set in 'comsoil_h'
    3434
    3535! Local variables:
    36 real,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data on extended lonxlat grid
    37 ! Note iip1,jjp1 known from paramet.h; nsoilmx known from dimphys.h
    38 real,dimension(iip1,jjp1) :: data2 ! to store 2D data on extended lonxlat grid
    39 real :: data0 ! to store 0D data
     36real*4,dimension(iip1,jjp1,nsoilmx) :: data3 ! to store 3D data
     37! Note iip1,jjp1 known from paramet.h; nsoilmx known from comsoil_h
     38real*4,dimension(iip1,jjp1) :: data2 ! to store 2D data
     39real*4 :: data0 ! to store 0D data
    4040integer :: i,j,l ! for loops
    4141integer :: ig0
    4242
    43 real,save :: date ! time counter (in elapsed days)
     43real*4,save :: date ! time counter (in elapsed days)
    4444integer,save :: isample ! sample rate at which data is to be written to output
    4545integer,save :: ntime=0 ! counter to internally store time steps
    … …  
    6363real dx2_glo(iim,jjp1)     ! to store a global 2D (surface) data set
    6464real px2(ngrid)
    65 #else
    66 logical,parameter :: is_master=.true.
    67 logical,parameter :: is_mpi_root=.true.
    6865#endif
    6966
    … …  
    8380 
    8481  ! Set output sample rate
    85   isample=int(ecritphy) ! same as for diagfi outputs
     82  isample=ecritphy ! same as for diagfi outputs
    8683  ! Note ecritphy is known from control.h
    8784 
    … …  
    128125     ierr=NF_INQ_VARID(nid,"time",varid)
    129126     ! Add the current value of date to the "time" array
    130 #ifdef NC_DOUBLE
    131      ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
    132 #else
     127!#ifdef NC_DOUBLE
     128!     ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
     129!#else
    133130     ierr=NF_PUT_VARA_REAL(nid,varid,ntime,1,date)
    134 #endif
     131!#endif
    135132     if (ierr.ne.NF_NOERR) then
    136133      write(*,*)"writediagsoil: Failed writing date to time variable"
    … …  
    203200 
    204201  ! B.3. Write the slab of data
    205 #ifdef NC_DOUBLE
    206   ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
    207 #else
     202!#ifdef NC_DOUBLE
     203!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
     204!#else
    208205  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3)
    209 #endif
     206!#endif
    210207  if (ierr.ne.NF_NOERR) then
    211208    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
    … …  
    273270 
    274271  ! B.3. Write the slab of data
    275 #ifdef NC_DOUBLE
    276   ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
    277 #else
     272!#ifdef NC_DOUBLE
     273!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
     274!#else
    278275  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2)
    279 #endif
     276!#endif
    280277  if (ierr.ne.NF_NOERR) then
    281278    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
    … …  
    312309
    313310  ! B.3. Write the data
    314 #ifdef NC_DOUBLE
    315   ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0)
    316 #else
     311!#ifdef NC_DOUBLE
     312!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0)
     313!#else
    317314  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data0)
    318 #endif
     315!#endif
    319316  if (ierr.ne.NF_NOERR) then
    320317    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&

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

    r993 r1216  
    4848      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    4949#endif
     50      use control_mod, only: ecritphy, iphysiq, day_step
    5051
    5152      implicit none
    … …  
    5556#include "dimphys.h"
    5657#include "paramet.h"
    57 #include "control.h"
     58!#include "control.h"
    5859#include "comvert.h"
    5960#include "comgeom.h"
    … …  
    114115!Sortie des variables au rythme voulu
    115116
    116       irythme = int(ecritphy)*iradia ! sortie au rythme de ecritphy*iradia
     117      irythme = ecritphy*iradia ! sortie au rythme de ecritphy*iradia
    117118!EM+JL if the spetra need to be output more frequently, need to define a ecritSpec...
    118119!     irythme = iphysiq  ! sortie a tous les pas physique

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

    r993 r1216  
    4848      use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    4949#endif
     50      use control_mod, only: ecritphy, iphysiq, day_step
    5051
    5152      implicit none
    … …  
    5556#include "dimphys.h"
    5657#include "paramet.h"
    57 #include "control.h"
     58!#include "control.h"
    5859#include "comvert.h"
    5960#include "comgeom.h"
    … …  
    114115!Sortie des variables au rythme voulu
    115116
    116       irythme = int(ecritphy)*iradia ! sortie au rythme de ecritphy
     117      irythme = ecritphy*iradia ! sortie au rythme de ecritphy
    117118!EM+JL if the spetra need to be output more frequently, need to define a ecritSpec...
    118119!     irythme = iphysiq  ! sortie a tous les pas physique

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

    r965 r1216  
    11subroutine wstats(ngrid,nom,titre,unite,dim,px)
    22
    3 #ifdef CPP_PARA
    43use mod_phys_lmdz_para, only : is_mpi_root, is_master, gather, klon_mpi_begin
    54use mod_grid_phy_lmdz, only : klon_glo, Grid1Dto2D_glo
    6 #endif
    75
    86implicit none
    … …  
    4745real px2_glop(ngrid) ! to store a 2D data set on global physics grid
    4846real px2_glo(iim,jjp1) ! to store a 2D data set on global lonxlat grid
    49 logical,parameter :: is_master=.true.
    50 logical,parameter :: is_mpi_root=.true.
    51 integer,parameter :: klon_mpi_begin=1
    5247#endif
    5348
    … …  
    290285endif ! of if (is_master)
    291286
    292 end
     287end subroutine wstats
    293288
    294289!======================================================
    … …  
    367362endif
    368363
    369 end
     364end subroutine inivar
     365
     366!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     367
     368subroutine def_var_stats(nid,name,title,units,nbdim,dimids,nvarid,ierr)
     369
     370! This subroutine defines variable 'name' in a (pre-existing and opened)
     371! NetCDF file (known from its NetCDF ID 'nid').
     372! The number of dimensions 'nbdim' of the variable, as well as the IDs of
     373! corresponding dimensions must be set (in array 'dimids').
     374! Upon successfull definition of the variable, 'nvarid' contains the
     375! NetCDF ID of the variable.
     376! The variables' attributes 'title' (Note that 'long_name' would be more
     377! appropriate) and 'units' are also set.
     378
     379implicit none
     380
     381#include "netcdf.inc"
     382
     383integer,intent(in) :: nid ! NetCDF file ID
     384character(len=*),intent(in) :: name ! the variable's name
     385character(len=*),intent(in) :: title ! 'title' attribute of variable
     386character(len=*),intent(in) :: units ! 'units' attribute of variable
     387integer,intent(in) :: nbdim ! number of dimensions of the variable
     388integer,dimension(nbdim),intent(in) :: dimids ! NetCDF IDs of the dimensions
     389                                              ! the variable is defined along
     390integer,intent(out) :: nvarid ! NetCDF ID of the variable
     391integer,intent(out) :: ierr ! returned NetCDF staus code
     392
     393! 1. Switch to NetCDF define mode
     394ierr=NF_REDEF(nid)
     395
     396! 2. Define the variable
     397#ifdef NC_DOUBLE
     398ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dimids,nvarid)
     399#else
     400ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dimids,nvarid)
     401#endif
     402if(ierr/=NF_NOERR) then
     403   write(*,*) "def_var_stats: Failed defining variable "//trim(name)
     404   write(*,*) NF_STRERROR(ierr)
     405   stop ""
     406endif
     407
     408! 3. Write attributes
     409ierr=NF_PUT_ATT_TEXT(nid,nvarid,"title",&
     410                     len_trim(adjustl(title)),adjustl(title))
     411if(ierr/=NF_NOERR) then
     412   write(*,*) "def_var_stats: Failed writing title attribute for "//trim(name)
     413   write(*,*) NF_STRERROR(ierr)
     414   stop ""
     415endif
     416
     417ierr=NF_PUT_ATT_TEXT(nid,nvarid,"units",&
     418                     len_trim(adjustl(units)),adjustl(units))
     419if(ierr/=NF_NOERR) then
     420   write(*,*) "def_var_stats: Failed writing units attribute for "//trim(name)
     421   write(*,*) NF_STRERROR(ierr)
     422   stop ""
     423endif
     424
     425! 4. Switch out of NetCDF define mode
     426ierr = NF_ENDDEF(nid)
     427
     428end subroutine def_var_stats
     429

  • trunk/LMDZ.GENERIC/makegcm_g95

    r988 r1216  
    1515set bands="32x36"
    1616set scatterers="1"
     17set full=""
    1718########################################################################
    1819# path a changer contenant les sources et les objets du modele
    … …  
    273274             latitudes and vertical layers respectively.
    274275
    275 -t ntrac   Selects the number of tracers present in the model
    276 
    277              Options -d and -t overwrite file
    278              $LMDGCM/libf/grid/dimensions.h
    279              which contains the 3 dimensions of the
    280              horizontal grid
    281              im, jm, lm plus the number of tracers passively advected
    282              by the dynamics ntrac,
    283              in 4 PARAMETER FORTRAN format
    284              with a new file:
    285              $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
    286              If the file does not exist already
    287              it is created by the script
    288              $LMDGCM/libf/grid/dimension/makdim
    289 
    290276-s nscat   Number of radiatively active scatterers
    291277
    … …  
    331317           there is no need to specify  -Ldirn.
    332318
     319-full      Full (re)compilation (from scratch)
     320
    333321eod
    334322exit
    … …  
    367355     case -olddyn
    368356        set dyntype="olddyn" ; shift; goto top
     357
     358     case -full
     359        set full="full" ; shift ; goto top
    369360
    370361     case -filtre
    … …  
    527518# Build the appropriate 'dimensions.h' file
    528519cd dimension
    529 ./makdim $ntrac $dim
     520./makdim $dim
    530521# echo contents of dimensions.h to standard output
    531522cat $libf/grid/dimensions.h
    … …  
    559550echo dimension $dimension dim $dim
    560551if ( $dimension == 1 ) then
    561   echo pas de dynamique
    562   set dyn="L_DYN= DYN= L_FILTRE= "
     552  echo "No dynamics"
     553##  set dyn="L_DYN= DYN= L_FILTRE= "
     554## NB: we still need to have L_DYN=libdyn3d to reach routines and module
     555## objects which are located in dyn3d
     556  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
    563557endif
    564558endif
    … …  
    588582
    589583echo "dimc $dimc"
     584
     585#cleanup for a full recompilation, if requested
     586if ("$full" == "full") then
     587# remove makefile and $libo
     588  cd $model
     589  \rm -f makefile
     590  \rm -rf $libo/*
     591endif
    590592
    591593########################################################################

  • trunk/LMDZ.GENERIC/makegcm_gfortran

    r989 r1216  
    1515set bands="32x36"
    1616set scatterers="1"
     17set full=""
    1718########################################################################
    1819# path a changer contenant les sources et les objets du modele
    … …  
    184185   set optim90="-O3 -funroll-loops "
    185186   set optimtru90="-O3 -funroll-loops "
    186    set opt_link=" -L$NCDFLIB -lnetcdf -lnetcdff "
     187   set opt_link=" -L$NCDFLIB -lnetcdf "
    187188
    188189   #NB: on gnome -O3 ==> NaNs ...
    … …  
    274275             latitudes and vertical layers respectively.
    275276
    276 -t ntrac   Selects the number of tracers present in the model
    277 
    278              Options -d and -t overwrite file
    279              $LMDGCM/libf/grid/dimensions.h
    280              which contains the 3 dimensions of the
    281              horizontal grid
    282              im, jm, lm plus the number of tracers passively advected
    283              by the dynamics ntrac,
    284              in 4 PARAMETER FORTRAN format
    285              with a new file:
    286              $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
    287              If the file does not exist already
    288              it is created by the script
    289              $LMDGCM/libf/grid/dimension/makdim
    290 
    291277-s nscat   Number of radiatively active scatterers
    292278
    … …  
    332318           there is no need to specify  -Ldirn.
    333319
     320-full      Full (re)compilation (from scratch)
     321
    334322eod
    335323exit
    … …  
    368356     case -olddyn
    369357        set dyntype="olddyn" ; shift; goto top
     358
     359     case -full
     360        set full="full" ; shift ; goto top
    370361
    371362     case -filtre
    … …  
    532523# Build the appropriate 'dimensions.h' file
    533524cd dimension
    534 ./makdim $ntrac $dim
     525./makdim $dim
    535526# echo contents of dimensions.h to standard output
    536527cat $libf/grid/dimensions.h
    … …  
    564555echo dimension $dimension dim $dim
    565556if ( $dimension == 1 ) then
    566   echo pas de dynamique
    567   set dyn="L_DYN= DYN= L_FILTRE= "
     557  echo "No dynamics"
     558##  set dyn="L_DYN= DYN= L_FILTRE= "
     559## NB: we still need to have L_DYN=libdyn3d to reach routines and module
     560## objects which are located in dyn3d
     561  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
    568562endif
    569563endif
    … …  
    593587
    594588echo "dimc $dimc"
     589
     590#cleanup for a full recompilation, if requested
     591if ("$full" == "full") then
     592# remove makefile and $libo
     593  cd $model
     594  \rm -f makefile
     595  \rm -rf $libo/*
     596endif
    595597
    596598########################################################################
    … …  
    657659   set f77="gfortran -ffree-line-length-264"
    658660   set f90="gfortran -ffree-line-length-264"
    659    set opt_link=" -L$LIBOGCM -L$NCDFLIB -lnetcdff -lnetcdf "
     661#   set opt_link=" -L$LIBOGCM -L$NCDFLIB -lnetcdff -lnetcdf "
     662   set opt_link=" -L$LIBOGCM -L$NCDFLIB -lnetcdf "
    660663else if $SUN then
    661664   set f77=f90

  • trunk/LMDZ.GENERIC/makegcm_ifort

    r988 r1216  
    1515set bands="32x36"
    1616set scatterers="1"
     17set full=""
    1718########################################################################
    1819# path a changer contenant les sources et les objets du modele
    … …  
    279280             latitudes and vertical layers respectively.
    280281
    281 -t ntrac   Selects the number of tracers present in the model
    282 
    283              Options -d and -t overwrite file
    284              $LMDGCM/libf/grid/dimensions.h
    285              which contains the 3 dimensions of the
    286              horizontal grid
    287              im, jm, lm plus the number of tracers passively advected
    288              by the dynamics ntrac,
    289              in 4 PARAMETER FORTRAN format
    290              with a new file:
    291              $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
    292              If the file does not exist already
    293              it is created by the script
    294              $LMDGCM/libf/grid/dimension/makdim
    295 
    296282-s nscat   Number of radiatively active scatterers
    297283
    … …  
    337323           there is no need to specify  -Ldirn.
    338324
     325-full      Full (re)compilation (from scratch)
     326
    339327eod
    340328exit
    … …  
    373361     case -olddyn
    374362        set dyntype="olddyn" ; shift; goto top
     363
     364     case -full
     365        set full="full" ; shift ; goto top
    375366
    376367     case -filtre
    … …  
    534525# Build the appropriate 'dimensions.h' file
    535526cd dimension
    536 ./makdim $ntrac $dim
     527./makdim $dim
    537528# echo contents of dimensions.h to standard output
    538529cat $libf/grid/dimensions.h
    … …  
    566557echo dimension $dimension dim $dim
    567558if ( $dimension == 1 ) then
    568   echo pas de dynamique
    569   set dyn="L_DYN= DYN= L_FILTRE= "
     559  echo "No dynamics"
     560##  set dyn="L_DYN= DYN= L_FILTRE= "
     561## NB: we still need to have L_DYN=libdyn3d to reach routines and module
     562## objects which are located in dyn3d
     563  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
    570564endif
    571565endif
    … …  
    595589
    596590echo "dimc $dimc"
     591
     592#cleanup for a full recompilation, if requested
     593if ("$full" == "full") then
     594# remove makefile and $libo
     595  cd $model
     596  \rm -f makefile
     597  \rm -rf $libo/*
     598endif
    597599
    598600########################################################################

  • trunk/LMDZ.GENERIC/makegcm_pgf90

    r988 r1216  
    1515set bands="32x36"
    1616set scatterers="1"
     17set full=""
    1718########################################################################
    1819# path a changer contenant les sources et les objets du modele
    … …  
    271272             latitudes and vertical layers respectively.
    272273
    273 -t ntrac   Selects the number of tracers present in the model
    274 
    275              Options -d and -t overwrite file
    276              $LMDGCM/libf/grid/dimensions.h
    277              which contains the 3 dimensions of the
    278              horizontal grid
    279              im, jm, lm plus the number of tracers passively advected
    280              by the dynamics ntrac,
    281              in 4 PARAMETER FORTRAN format
    282              with a new file:
    283              $LMDGCM/libf/grid/dimension/dimensions.im.jm.lm.tntrac
    284              If the file does not exist already
    285              it is created by the script
    286              $LMDGCM/libf/grid/dimension/makdim
    287 
    288274-s nscat   Number of radiatively active scatterers
    289275
    … …  
    329315           there is no need to specify  -Ldirn.
    330316
     317-full      Full (re)compilation (from scratch)
     318
    331319eod
    332320exit
    … …  
    365353     case -olddyn
    366354        set dyntype="olddyn" ; shift; goto top
     355
     356     case -full
     357        set full="full" ; shift ; goto top
    367358
    368359     case -filtre
    … …  
    529520# Build the appropriate 'dimensions.h' file
    530521cd dimension
    531 ./makdim $ntrac $dim
     522./makdim $dim
    532523# echo contents of dimensions.h to standard output
    533524cat $libf/grid/dimensions.h
    … …  
    561552echo dimension $dimension dim $dim
    562553if ( $dimension == 1 ) then
    563   echo pas de dynamique
    564   set dyn="L_DYN= DYN= L_FILTRE= "
     554  echo "No dynamics"
     555##  set dyn="L_DYN= DYN= L_FILTRE= "
     556## NB: we still need to have L_DYN=libdyn3d to reach routines and module
     557## objects which are located in dyn3d
     558  set dyn="L_DYN=-ldyn3d DYN= L_FILTRE= DIRMAIN=phy$physique "
    565559endif
    566560endif
    … …  
    590584
    591585echo "dimc $dimc"
     586
     587#cleanup for a full recompilation, if requested
     588if ("$full" == "full") then
     589# remove makefile and $libo
     590  cd $model
     591  \rm -f makefile
     592  \rm -rf $libo/*
     593endif
    592594
    593595########################################################################
Note: See TracChangeset for help on using the changeset viewer.