source: LMDZ5/trunk/libf/dyn3d/dynredem.F90 @ 2293

SUBROUTINE dynredem0(fichnom,iday_end,phis)
!
!-------------------------------------------------------------------------------
! Write the NetCDF restart file (initialization).
!-------------------------------------------------------------------------------
#ifdef CPP_IOIPSL
  USE IOIPSL
#endif
  USE infotrac
  USE netcdf, ONLY:   NF90_CREATE, NF90_DEF_DIM, NF90_REDEF,  NF90_INQ_VARID, &
      NF90_CLOBBER,   NF90_CLOSE,  NF90_DEF_VAR, NF90_ENDDEF, NF90_PUT_ATT,   &
      NF90_UNLIMITED, NF90_GLOBAL, NF90_FLOAT,   NF90_DOUBLE
  USE netcdf95, ONLY: NF95_PUT_VAR
  IMPLICIT NONE
  include "dimensions.h"
  include "paramet.h"
  include "comconst.h"
  include "comvert.h"
  include "comgeom2.h"
  include "temps.h"
  include "ener.h"
  include "logic.h"
  include "netcdf.inc"
  include "description.h"
  include "serre.h"
  include "iniprint.h"
!===============================================================================
! Arguments:
  CHARACTER(LEN=*), INTENT(IN) :: fichnom          !--- FILE NAME
  INTEGER,          INTENT(IN) :: iday_end         !--- 
  REAL,             INTENT(IN) :: phis(iip1, jjp1) !--- GROUND GEOPOTENTIAL
!===============================================================================
! Local variables:
  INTEGER :: iq, l
  INTEGER, PARAMETER :: length=100
  REAL    :: tab_cntrl(length)                     !--- RUN PARAMETERS TABLE
  CHARACTER(LEN=20) :: modname
!   For NetCDF:
  CHARACTER(LEN=30) :: unites
  INTEGER :: indexID
  INTEGER :: rlonuID, rlonvID, rlatuID, rlatvID
  INTEGER :: sID, sigID, nID, vID, timID
  INTEGER :: yyears0, jjour0, mmois0
  REAL :: zan0, zjulian, hours
!===============================================================================
  modname='dynredem0'
#ifdef CPP_IOIPSL
  CALL ymds2ju(annee_ref, 1, iday_end, 0.0, zjulian)
  CALL ju2ymds(zjulian, yyears0, mmois0, jjour0, hours)
#else
! set yyears0, mmois0, jjour0 to 0,1,1 (hours is not used)
  yyears0=0
  mmois0=1
  jjour0=1
#endif        

  tab_cntrl(:)  = 0.
  tab_cntrl(1)  = REAL(iim)
  tab_cntrl(2)  = REAL(jjm)
  tab_cntrl(3)  = REAL(llm)
  tab_cntrl(4)  = REAL(day_ref)
  tab_cntrl(5)  = REAL(annee_ref)
  tab_cntrl(6)  = rad
  tab_cntrl(7)  = omeg
  tab_cntrl(8)  = g
  tab_cntrl(9)  = cpp
  tab_cntrl(10) = kappa
  tab_cntrl(11) = daysec
  tab_cntrl(12) = dtvr
  tab_cntrl(13) = etot0
  tab_cntrl(14) = ptot0
  tab_cntrl(15) = ztot0
  tab_cntrl(16) = stot0
  tab_cntrl(17) = ang0
  tab_cntrl(18) = pa
  tab_cntrl(19) = preff

!    .....    parameters for zoom    ......   
  tab_cntrl(20) = clon
  tab_cntrl(21) = clat
  tab_cntrl(22) = grossismx
  tab_cntrl(23) = grossismy
!
  IF ( fxyhypb )   THEN
    tab_cntrl(24) = 1.
    tab_cntrl(25) = dzoomx
    tab_cntrl(26) = dzoomy
    tab_cntrl(27) = 0.
    tab_cntrl(28) = taux
    tab_cntrl(29) = tauy
  ELSE
    tab_cntrl(24) = 0.
    tab_cntrl(25) = dzoomx
    tab_cntrl(26) = dzoomy
    tab_cntrl(27) = 0.
    tab_cntrl(28) = 0.
    tab_cntrl(29) = 0.
    IF( ysinus )  tab_cntrl(27) = 1.
  END IF
  tab_cntrl(30) = REAL(iday_end)
  tab_cntrl(31) = REAL(itau_dyn + itaufin)
! start_time: start_time of simulation (not necessarily 0.)
  tab_cntrl(32) = start_time
!.........................................................

!--- File creation
  CALL err(NF90_CREATE(fichnom,NF90_CLOBBER,nid))

!--- Some global attributes
  CALL err(NF90_PUT_ATT(nid,NF90_GLOBAL,"title","Fichier demarrage dynamique"))

!--- Dimensions
  CALL err(NF90_DEF_DIM(nid,"index", length, indexID))
  CALL err(NF90_DEF_DIM(nid,"rlonu", iip1,   rlonuID))
  CALL err(NF90_DEF_DIM(nid,"rlatu", jjp1,   rlatuID))
  CALL err(NF90_DEF_DIM(nid,"rlonv", iip1,   rlonvID))
  CALL err(NF90_DEF_DIM(nid,"rlatv", jjm,    rlatvID))
  CALL err(NF90_DEF_DIM(nid,"sigs",  llm,        sID))
  CALL err(NF90_DEF_DIM(nid,"sig",   llmp1,    sigID))
  CALL err(NF90_DEF_DIM(nid,"temps", NF90_UNLIMITED, timID))

!--- Define and save invariant fields
  CALL put_var1("controle","Parametres de controle" ,[indexID],tab_cntrl)
  CALL put_var1("rlonu"   ,"Longitudes des points U",[rlonuID],rlonu)
  CALL put_var1("rlatu"   ,"Latitudes des points U" ,[rlatuID],rlatu)
  CALL put_var1("rlonv"   ,"Longitudes des points V",[rlonvID],rlonv)
  CALL put_var1("rlatv"   ,"Latitudes des points V" ,[rlatvID],rlatv)
  CALL put_var1("nivsigs" ,"Numero naturel des couches s"    ,[sID]  ,nivsigs)
  CALL put_var1("nivsig"  ,"Numero naturel des couches sigma",[sigID],nivsig)
  CALL put_var1("ap"      ,"Coefficient A pour hybride"      ,[sigID],ap)
  CALL put_var1("bp"      ,"Coefficient B pour hybride"      ,[sigID],bp)
  CALL put_var1("presnivs",""                                ,[sID]  ,presnivs)
! covariant <-> contravariant <-> natural conversion coefficients
  CALL put_var2("cu","Coefficient de passage pour U",[rlonuID,rlatuID],cu)
  CALL put_var2("cv","Coefficient de passage pour V",[rlonvID,rlatvID],cv)
  CALL put_var2("aire","Aires de chaque maille"     ,[rlonvID,rlatuID],aire)
  CALL put_var2("phisinit","Geopotentiel au sol"    ,[rlonvID,rlatuID],phis)

!--- Define fields saved later
  WRITE(unites,"('days since ',i4,'-',i2.2,'-',i2.2,' 00:00:00')"),&
               yyears0,mmois0,jjour0
  CALL put_var0("temps","Temps de simulation",[timID],unites)
  CALL put_var0("ucov" ,"Vitesse U"  ,[rlonuID,rlatuID,sID,timID])
  CALL put_var0("vcov" ,"Vitesse V"  ,[rlonvID,rlatvID,sID,timID])
  CALL put_var0("teta" ,"Temperature",[rlonvID,rlatuID,sID,timID])
  DO iq=1,nqtot
    CALL put_var0(tname(iq),ttext(iq) ,[rlonvID,rlatuID,sID,timID])
  END DO
  CALL put_var0("masse","Masse d air"    ,[rlonvID,rlatuID,sID,timID])
  CALL put_var0("ps"   ,"Pression au sol",[rlonvID,rlatuID    ,timID])
  CALL err(NF90_ENDDEF(nid))
  CALL err(NF90_CLOSE (nid))

  WRITE(lunout,*)TRIM(modname)//': iim,jjm,llm,iday_end',iim,jjm,llm,iday_end
  WRITE(lunout,*)TRIM(modname)//': rad,omeg,g,cpp,kappa',rad,omeg,g,cpp,kappa


CONTAINS


SUBROUTINE put_var0(var,title,did,units)
  CHARACTER(LEN=*),           INTENT(IN) :: var, title
  INTEGER,                    INTENT(IN) :: did(:)
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
#ifdef NC_DOUBLE
  CALL err(NF90_DEF_VAR(nid,var,NF90_DOUBLE,did,vID),var)
#else
  CALL err(NF90_DEF_VAR(nid,var,NF90_FLOAT,did,vID),var)
#endif
  IF(title/="")      CALL err(NF90_PUT_ATT(nid,vID,"title",title),var)
  IF(PRESENT(units)) CALL err(NF90_PUT_ATT(nid,vID,"units",units),var)
END SUBROUTINE put_var0


SUBROUTINE put_var1(var,title,did,v,units)
  CHARACTER(LEN=*),           INTENT(IN) :: var, title
  INTEGER,                    INTENT(IN) :: did(1)
#ifdef NC_DOUBLE
  DOUBLE PRECISION,           INTENT(IN) :: v(:)
#else
  REAL,                       INTENT(IN) :: v(:)
#endif
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
#ifdef NC_DOUBLE
  CALL err(NF90_DEF_VAR(nid,var,NF90_DOUBLE,did,vID),var)
#else
  CALL err(NF90_DEF_VAR(nid,var,NF90_FLOAT,did,vID),var)
#endif
  IF(title/="")      CALL err(NF90_PUT_ATT(nid,vID,"title",title),var)
  IF(PRESENT(units)) CALL err(NF90_PUT_ATT(nid,vID,"units",units),var)
  CALL err(NF90_ENDDEF(nid))
  CALL NF95_PUT_VAR(nid,vID,v)
  CALL err(NF90_REDEF(nid))
END SUBROUTINE put_var1


SUBROUTINE put_var2(var,title,did,v,units)
  CHARACTER(LEN=*),           INTENT(IN) :: var, title
  INTEGER,                    INTENT(IN) :: did(2)
#ifdef NC_DOUBLE
  DOUBLE PRECISION,           INTENT(IN) :: v(:,:)
#else
  REAL,                       INTENT(IN) :: v(:,:)
#endif
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: units
#ifdef NC_DOUBLE
  CALL err(NF90_DEF_VAR(nid,var,NF90_DOUBLE,did,vID),var)
#else
  CALL err(NF90_DEF_VAR(nid,var,NF90_FLOAT,did,vID),var)
#endif
  IF(title/="")      CALL err(NF90_PUT_ATT(nid,vID,"title",title),var)
  IF(PRESENT(units)) CALL err(NF90_PUT_ATT(nid,vID,"units",units),var)
  CALL err(NF90_ENDDEF(nid))
  CALL NF95_PUT_VAR(nid,vID,v)
  CALL err(NF90_REDEF(nid))

END SUBROUTINE put_var2


SUBROUTINE err(ierr,var)
  USE netcdf, ONLY: NF90_STRERROR, NF90_NOERR
  IMPLICIT NONE
  INTEGER,                    INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: var    !--- VARIABLE NAME
  CHARACTER(LEN=256) :: file, msg
  IF(ierr==NF90_NoERR) RETURN
  msg='Error in "'//TRIM(modname)//'" for file "'//TRIM(fichnom)//'"'
  IF(PRESENT(var)) msg=TRIM(msg)//'" and variable "'//TRIM(var)//'"'
  WRITE(lunout,*)TRIM(msg)//': '//NF90_STRERROR(ierr)

END SUBROUTINE err

END SUBROUTINE dynredem0
!
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!
SUBROUTINE dynredem1(fichnom,time,vcov,ucov,teta,q,masse,ps)
!
!-------------------------------------------------------------------------------
! Purpose: Write the NetCDF restart file (append).
!-------------------------------------------------------------------------------
  USE infotrac
  USE control_mod
  USE netcdf,   ONLY: NF90_OPEN,  NF90_NOWRITE, NF90_INQ_VARID, NF90_NoErr,    &
                      NF90_CLOSE, NF90_WRITE,   NF90_GET_VAR
  USE netcdf95, ONLY: NF95_PUT_VAR
  USE assert_eq_m, ONLY: assert_eq
  IMPLICIT NONE
#include "dimensions.h"
#include "paramet.h"
#include "description.h"
#include "comvert.h"
#include "comgeom.h"
#include "temps.h"
#include "iniprint.h"
!===============================================================================
! Arguments:
  CHARACTER(LEN=*), INTENT(IN) :: fichnom           !-- FILE NAME
  REAL, INTENT(IN) ::  time                         !-- TIME
  REAL, INTENT(IN) ::  vcov(iip1,jjm, llm)          !-- V COVARIANT WIND
  REAL, INTENT(IN) ::  ucov(iip1,jjp1,llm)          !-- U COVARIANT WIND
  REAL, INTENT(IN) ::  teta(iip1,jjp1,llm)          !-- POTENTIAL TEMPERATURE
  REAL, INTENT(IN) ::     q(iip1,jjp1,llm,nqtot)    !-- TRACERS
  REAL, INTENT(IN) :: masse(iip1,jjp1,llm)          !-- MASS PER CELL
  REAL, INTENT(IN) ::    ps(iip1,jjp1)              !-- GROUND PRESSURE
!===============================================================================
! Local variables:
  INTEGER :: l, iq, nid, vID, nid_trac, vID_trac
  INTEGER, SAVE :: nb=0
  INTEGER, PARAMETER :: length=100
#ifdef NC_DOUBLE
  DOUBLE PRECISION   :: trac_tmp(ip1jmp1,llm)
#else
  REAL               :: trac_tmp(ip1jmp1,llm)
#endif
  REAL               :: tab_cntrl(length) ! tableau des parametres du run
  CHARACTER(LEN=256) :: modname, var, fil
  LOGICAL            :: exist_file
!===============================================================================
  modname='dynredem1'
  fil=fichnom
  CALL err(NF90_OPEN(fil,NF90_WRITE,nid),"open",fil)

!--- Write/extend time coordinate
  nb = nb + 1
  CALL sav_var1("temps",[time],nb)
  WRITE(lunout,*)TRIM(modname)//": Saving for ", nb, time

!--- Rewrite control table (itaufin undefined in dynredem0)
  var="controle"
  CALL get_var1(var,tab_cntrl); tab_cntrl(31)=DBLE(itau_dyn + itaufin)
  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
  CALL NF95_PUT_VAR(nid,vID,tab_cntrl)

!--- Save fields
  CALL sav_var3("ucov",ucov)
  CALL sav_var3("vcov",vcov)
  CALL sav_var3("teta",teta)
  CALL sav_var3("masse",masse)
  CALL sav_var2("ps"   ,ps)

!--- Tracers in file "start_trac.nc" (added by Anne)
  IF (type_trac == 'inca') THEN
    fil="start_trac.nc"; INQUIRE(FILE=fil,EXIST=exist_file)
    IF(.NOT.exist_file) CALL war(-1,"open",fil)
  END IF
  DO iq=1,nqtot; var=tname(iq)

  !--- Usual case
    IF(type_trac/='inca') THEN
      CALL sav_var3(var,q(:,:,:,iq)); CYCLE
    END IF

  !--- Special case for INCA tracer read from "start_trac.nc"
    IF(NF90_INQ_VARID(nid_trac,var,vID_trac)/=NF90_NoErr) THEN
      CALL war(-1,"inq",var,fil)
      CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var,fil)
      CALL NF95_PUT_VAR(nid,vID,q(:,:,:,iq))
    ELSE
      WRITE(lunout,*)TRIM(modname)//": <"//TRIM(var)//"> found in "//fil
      CALL err(NF90_GET_VAR(nid_trac,vID_trac,trac_tmp),"get",var,fil)
    END IF
    CALL sav_var3(var,RESHAPE(trac_tmp,SHAPE=[iip1,jjp1,llm]))
  END DO
  CALL err(NF90_CLOSE(nid),"close",fichnom)


CONTAINS


SUBROUTINE get_var1(var,v)
  CHARACTER(LEN=*), INTENT(IN) :: var
  REAL,             INTENT(OUT) :: v(:)
  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
  CALL err(NF90_GET_VAR(nid,vID,v),"get",var)
END SUBROUTINE get_var1


SUBROUTINE sav_var1(var,v,start)
  CHARACTER(LEN=*),  INTENT(IN) :: var
#ifdef NC_DOUBLE
  DOUBLE PRECISION,  INTENT(IN) :: v(:)
#else
  REAL,              INTENT(IN) :: v(:)
#endif
  INTEGER, OPTIONAL, INTENT(IN) :: start
  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
  IF(PRESENT(start)) THEN
    CALL NF95_PUT_VAR(nid,vID,v,start=[start])
  ELSE
    CALL NF95_PUT_VAR(nid,vID,v)
  END IF
END SUBROUTINE sav_var1


SUBROUTINE sav_var2(var,v)
  CHARACTER(LEN=*), INTENT(IN) :: var
#ifdef NC_DOUBLE
  DOUBLE PRECISION, INTENT(IN) :: v(:,:)
#else
  REAL,             INTENT(IN) :: v(:,:)
#endif
  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
  CALL NF95_PUT_VAR(nid,vID,v)
END SUBROUTINE sav_var2


SUBROUTINE sav_var3(var,v)
  CHARACTER(LEN=*), INTENT(IN) :: var
#ifdef NC_DOUBLE
  DOUBLE PRECISION, INTENT(IN) :: v(:,:,:)
#else
  REAL,             INTENT(IN) :: v(:,:,:)
#endif

print*,'var='//TRIM(var)
print*,SHAPE(v)
  CALL err(NF90_INQ_VARID(nid,var,vID),"inq",var)
  CALL NF95_PUT_VAR(nid,vID,v)
END SUBROUTINE sav_var3


FUNCTION msg(typ,nam,fil)
  IMPLICIT NONE
  CHARACTER(LEN=256)           :: msg    !--- STANDARDIZED MESSAGE
  CHARACTER(LEN=*), INTENT(IN) :: typ    !--- TYPE OF OPERATION
  CHARACTER(LEN=*), INTENT(IN) :: nam    !--- FIELD NAME
  CHARACTER(LEN=*), INTENT(IN) :: fil    !--- FILE  NAME
  SELECT CASE(typ)
    CASE('inq');   msg="Missing field <"//TRIM(nam)//">"
    CASE('get');   msg="Reading failed for <"//TRIM(nam)//">"
    CASE('open');  msg="Opening failed for <"//TRIM(nam)//">"
    CASE('close'); msg="Closing failed for <"//TRIM(nam)//">"
  END SELECT
  msg=TRIM(modname)//": "//TRIM(msg)
  IF(typ=="inq".AND.fil/="") msg=TRIM(msg)//" in file <"//TRIM(fil)//">"

END FUNCTION msg


SUBROUTINE err(ierr,typ,nam,fil)
  IMPLICIT NONE
  INTEGER,                    INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
  CHARACTER(LEN=*),           INTENT(IN) :: typ    !--- TYPE OF OPERATION
  CHARACTER(LEN=*),           INTENT(IN) :: nam    !--- FIELD NAME
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fil    !--- FILE  NAME
  CHARACTER(LEN=256) :: file
  IF(ierr==NF90_NoERR) RETURN
  file=""; IF(PRESENT(fil)) file=fil
  CALL ABORT_gcm(modname,msg(typ,nam,file),ierr)
END SUBROUTINE err


SUBROUTINE war(ierr,typ,nam,fil)
  IMPLICIT NONE
  INTEGER,                    INTENT(IN) :: ierr   !--- NetCDF ERROR CODE
  CHARACTER(LEN=*),           INTENT(IN) :: typ    !--- TYPE OF OPERATION
  CHARACTER(LEN=*),           INTENT(IN) :: nam    !--- FIELD NAME
  CHARACTER(LEN=*), OPTIONAL, INTENT(IN) :: fil    !--- FILE  NAME
  CHARACTER(LEN=256) :: file
  IF(ierr==NF90_NoERR) RETURN
  file=""; IF(PRESENT(fil)) file=fil
  WRITE(lunout,*)msg(typ,nam,file)
END SUBROUTINE war


END SUBROUTINE dynredem1