source: trunk/LMDZ.GENERIC/libf/phystd/physdem1.F @ 787

      subroutine physdem1(ngrid,filename,lonfi,latfi,nsoil,nq,
     .                   phystep,day_ini,
     .                   time,tsurf,tsoil,emis,q2,qsurf,
     .                   airefi,alb,ith,pzmea,pzstd,pzsig,pzgam,pzthe,
     .                   cloudfrac,totcloudfrac,hice,nametrac)

      use surfdat_h
      use comsoil_h
      USE comgeomfi_h
 
      implicit none

!==================================================================
!     
!     Purpose
!     -------
!     Create physics (re-)start data file "restartfi.nc"
!     
!     Called by
!     ---------
!     rcm1d.F90
!     physiq.F90
!     newstart.F
!     
!     Calls
!     -----
!     none
!
!==================================================================

#include "dimensions.h"
#include "paramet.h"
#include "comvert.h"
#include "comgeom2.h"
#include "control.h"
#include "comdissnew.h"
#include "logic.h"
#include "ener.h"
#include "netcdf.inc"
#include "dimphys.h"
#include"callkeys.h"

      INTEGER :: ngrid

      INTEGER nid,iq
      INTEGER, parameter :: ivap=1
      REAL, parameter :: qsolmax= 150.0
      character (len=*) :: filename
      character (len=7) :: str7

      REAL day_ini
      INTEGER nsoil,nq
      integer ierr,idim1,idim2,idim3,idim4,idim5,idim6,nvarid

      REAL phystep,time
      REAL latfi(ngrid), lonfi(ngrid)
!      REAL champhys(ngrid)
      REAL tsurf(ngrid)
      INTEGER length
      PARAMETER (length=100)
      REAL tab_cntrl(length)

!     added by BC
      REAL hice(ngrid),cloudfrac(ngrid,nlayermx)
      REAL totcloudfrac(ngrid)

!     AS: name of tracers added as an argument to avoid using nqmx in commons (i.e. adv trac)
!       nota: physdem1 is used both in physiq and newstart.
      character*20 nametrac(nq)   ! name of the tracer


!      EXTERNAL defrun_new,iniconst,geopot,inigeom,massdair,pression
!      EXTERNAL exner_hyb , SSUM

#include "serre.h"
#include "clesph0.h"
#include "fxyprim.h"
#include "planete.h"
#include "comcstfi.h"

      real tsoil(ngrid,nsoil),emis(ngrid)
      real q2(ngrid, llm+1),qsurf(ngrid,nq)
      real airefi(ngrid)
      real alb(ngrid),ith(ngrid,nsoil)
      real pzmea(ngrid),pzstd(ngrid)
      real pzsig(ngrid),pzgam(ngrid),pzthe(ngrid)
      integer ig

! flag which identifies if we are using old tracer names (qsurf01,...)
      logical :: oldtracernames=.false.
      integer :: count
      character(len=30) :: txt ! to store some text
! indexes of water vapour & water ice tracers (if any):
      integer :: i_h2o_vap=0
      integer :: i_h2o_ice=0
c-----------------------------------------------------------------------

      ! copy airefi(:) to area(:)
      CALL SCOPY(ngrid,airefi,1,area,1)
      ! note: area() is defined in comgeomfi.h

      DO ig=1,ngrid
         albedodat(ig)=alb(ig) ! note: albedodat() is defined in surfdat.h
         zmea(ig)=pzmea(ig) ! note: zmea() is defined in surfdat.h
         zstd(ig)=pzstd(ig) ! note: zstd() is defined in surfdat.h
         zsig(ig)=pzsig(ig) ! note: zsig() is defined in surfdat.h
         zgam(ig)=pzgam(ig) ! note: zgam() is defined in surfdat.h
         zthe(ig)=pzthe(ig) ! note: zthe() is defined in surfdat.h
      ENDDO

      inertiedat(:,:)=ith(:,:) ! note inertiedat() is defined in comsoil.h
c
c  things to store in the physics start file:
c
      ierr = NF_CREATE(adjustl(filename),NF_CLOBBER, nid)
      IF (ierr.NE.NF_NOERR) THEN
        WRITE(6,*)'physdem1: Problem creating file ',adjustl(filename)
        write(6,*) NF_STRERROR(ierr)
        CALL ABORT
      ENDIF
c
      print*,'we got this far'

      ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 18,
     .                       "Physics start file")
c
      ierr = NF_DEF_DIM (nid,"index",length,idim1)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining index dimension'
        write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"physical_points",ngrid,idim2)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining physical_points dimension'
        write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoil,idim3)
      if (ierr.ne.NF_NOERR) then
      WRITE(6,*)'physdem1: Problem defining subsurface_layers dimension'
        write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"nlayer_plus_1",llm+1,idim4)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining nlayer+1 dimension'
        write(6,*) NF_STRERROR(ierr)
        call abort
      endif
c
      ierr = NF_DEF_DIM (nid,"number_of_advected_fields",nq,idim5)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining advected fields dimension'
        WRITE(6,*)' nq = ',nq,' and ierr = ', ierr
        write(6,*) NF_STRERROR(ierr)
      endif

      ierr = NF_DEF_DIM (nid,"nlayer",llm,idim6)
      if (ierr.ne.NF_NOERR) then
        WRITE(6,*)'physdem1: Problem defining nlayer dimension'
        write(6,*) NF_STRERROR(ierr)
        call abort
      endif

      ierr = NF_ENDDEF(nid) ! exit NetCDF define mode

c clear tab_cntrl(:) array
      DO ierr = 1, length
         tab_cntrl(ierr) = 0.0
      ENDDO

      write(*,*) "physdem1: ngrid: ",ngrid

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Fill control array tab_cntrl(:) with paramleters for this run
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c Informations on the physics grid
      tab_cntrl(1) = float(ngrid)  ! number of nodes on physics grid
      tab_cntrl(2) = float(nlayermx) ! number of atmospheric layers
      tab_cntrl(3) = day_ini + int(time)         ! initial day 
      tab_cntrl(4) = time -int(time)            ! initiale time of day

c Informations about Mars, used by dynamics and physics
      tab_cntrl(5) = rad      ! radius of Mars (m) ~3397200
      tab_cntrl(6) = omeg     ! rotation rate (rad.s-1)
      tab_cntrl(7) = g        ! gravity (m.s-2) ~3.72
      tab_cntrl(8) = mugaz    ! Molar mass of the atmosphere (g.mol-1) ~43.49
      tab_cntrl(9) = rcp      !  = r/cp  ~0.256793 (=kappa dans dynamique)
      tab_cntrl(10) = daysec  ! length of a sol (s)  ~88775

      tab_cntrl(11) = phystep  ! time step in the physics
      tab_cntrl(12) = 0.
      tab_cntrl(13) = 0.

c Informations about Mars, only for physics
      tab_cntrl(14) = year_day  ! length of year (sols) ~668.6
      tab_cntrl(15) = periastr  ! min. star-planet distance (AU)
      tab_cntrl(16) = apoastr   ! max. star-planet distance (AU) 
      tab_cntrl(17) = peri_day  ! date of periastron (sols since N. spring)
      tab_cntrl(18) = obliquit  ! Obliquity of the planet (deg) ~23.98

c Boundary layer and turbulence
      tab_cntrl(19) = z0        ! surface roughness (m) ~0.01
      tab_cntrl(20) = lmixmin   ! mixing length ~100
      tab_cntrl(21) = emin_turb ! minimal energy ~1.e-8

c Optical properties of polar caps and ground emissivity
      tab_cntrl(22) = albedice(1)  ! Albedo of northern cap ~0.5
      tab_cntrl(23) = albedice(2)  ! Albedo of southern cap ~0.5
      tab_cntrl(24) = emisice(1)   ! Emissivity of northern cap ~0.95
      tab_cntrl(25) = emisice(2)   ! Emissivity of southern cap ~0.95
      tab_cntrl(26) = emissiv      ! Emissivity of martian soil ~.95
      tab_cntrl(31) = iceradius(1) ! mean scat radius of CO2 snow (north)
      tab_cntrl(32) = iceradius(2) ! mean scat radius of CO2 snow (south)
      tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
      tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)

      tab_cntrl(28) = 0. 
      tab_cntrl(29) = 0.
      tab_cntrl(30) = 0.

! Soil properties:
      tab_cntrl(35) = volcapa ! soil volumetric heat capacity
      

!      write(*,*) "physdem1: tab_cntrl():",tab_cntrl
      
      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
         CALL abort
      ENDIF
      ! define variable
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "controle", NF_DOUBLE, 1, idim1,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "controle", NF_FLOAT, 1, idim1,nvarid)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to define controle'
         CALL abort
      ENDIF
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 18,
     .                        "Control parameters")
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to define controle title attribute'
         CALL abort
      ENDIF
      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to swich out of NetCDF define mode'
         CALL abort
      ENDIF
      ! write variable
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tab_cntrl)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,tab_cntrl)
#endif
      IF (ierr.NE.NF_NOERR) THEN
         PRINT*, 'physdem1: Failed to write controle data'
         CALL abort
      ENDIF

! write mid-layer depths mlayer() !known from comsoil.h

      ierr = NF_REDEF (nid) ! Enter NetCDF (re-)define mode
      ! define variable
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"soildepth",NF_DOUBLE,1,idim3,nvarid)
#else
      ierr = NF_DEF_VAR (nid,"soildepth",NF_FLOAT,1,idim3,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 20,
     .                        "Soil mid-layer depth")
      ierr = NF_ENDDEF(nid) ! Leave NetCDF define mode
      ! write variable
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,mlayer)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,mlayer)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "longitude", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "longitude", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 26,
     .                        "Longitudes of physics grid")
      ierr = NF_ENDDEF(nid)

#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,lonfi)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,lonfi)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "latitude", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "latitude", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid,nvarid,"title", 25,
     .                        "Latitudes of physics grid")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,latfi)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,latfi)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "area", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "area", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 9,
     .                        "Mesh area")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,area)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,area)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "phisfi", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "phisfi", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 27,
     .                        "Geopotential at the surface")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,phisfi)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,phisfi)
#endif

c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "albedodat", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "albedodat", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
     .                        "Albedo of bare ground")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,albedodat)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,albedodat)
#endif

c
c   some data for Francois Lott's programs
c

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZMEA", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZMEA", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
     .                        "Relief: mean relief")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zmea)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zmea)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZSTD", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZSTD", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 26,
     .                        "Relief: standard deviation")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zstd)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zstd)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZSIG", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZSIG", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
     .                        "Relief: sigma parameter")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zsig)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zsig)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZGAM", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZGAM", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
     .                        "Relief: gamma parameter")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zgam)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zgam)
#endif
c
      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "ZTHE", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "ZTHE", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 23,
     .                        "Relief: theta parameter")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,zthe)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,zthe)
#endif

c Write the physical fields


! Soil Thermal inertia

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"inertiedat",NF_DOUBLE,
     &                   2,(/idim2,idim3/),nvarid)
#else
      ierr = NF_DEF_VAR (nid,"inertiedat",NF_FLOAT,
     &                   2,(/idim2,idim3/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 20,
     .                        "Soil thermal inertia")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,inertiedat)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,inertiedat)
#endif




!  Surface temperature

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "tsurf", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "tsurf", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 19,
     .                        "Surface temperature")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsurf)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsurf)
#endif

! Soil temperature

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"tsoil",NF_DOUBLE,2,(/idim2,idim3/),nvarid)
#else
!      ierr = NF_DEF_VAR (nid, "tsoil", NF_FLOAT, 2, idim2,nvarid)
      ierr = NF_DEF_VAR (nid,"tsoil",NF_FLOAT,2,(/idim2,idim3/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 16,
     .                        "Soil temperature")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,tsoil)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,tsoil)
#endif

c emissivity

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "emis", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "emis", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 18,
     .                        "Surface emissivity")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,emis)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,emis)
#endif

c planetary boundary layer

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"q2", NF_DOUBLE, 2, (/idim2,idim4/),nvarid)
#else
      ierr = NF_DEF_VAR (nid, "q2", NF_FLOAT, 2,(/idim2,idim4/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
     .                        "pbl wind variance")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,q2)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,q2)
#endif
      IF (ierr.NE.NF_NOERR) THEN
        PRINT*, 'physdem1: Failed to write q2'
        CALL abort
      ENDIF







c cloud fraction (added by BC 2010)

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_DOUBLE, 2,
     & (/idim2,idim6/),nvarid)
#else
      ierr = NF_DEF_VAR (nid, "cloudfrac", NF_FLOAT, 2,
     & (/idim2,idim6/),nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
     .                        "Cloud fraction")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,cloudfrac)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,cloudfrac)
#endif

c total cloud fraction (added by BC 2010)

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid,"totcloudfrac", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "totcloudfrac", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 14,
     .                        "Total fraction")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,totcloudfrac)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,totcloudfrac)
#endif

c oceanic ice (added by BC 2010)

      ierr = NF_REDEF (nid)
#ifdef NC_DOUBLE
      ierr = NF_DEF_VAR (nid, "hice", NF_DOUBLE, 1, idim2,nvarid)
#else
      ierr = NF_DEF_VAR (nid, "hice", NF_FLOAT, 1, idim2,nvarid)
#endif
      ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", 21,
     .                        "Height of oceanic ice")
      ierr = NF_ENDDEF(nid)
#ifdef NC_DOUBLE
      ierr = NF_PUT_VAR_DOUBLE (nid,nvarid,hice)
#else
      ierr = NF_PUT_VAR_REAL (nid,nvarid,hice)
#endif











c tracers

! Preliminary stuff: check if tracers follow old naming convention (qsurf01,
!                    qsurf02, ...)
      count=0
      do iq=1,nq
        txt= " "
        write(txt,'(a1,i2.2)')'q',iq
        if (txt.ne.nametrac(iq)) then ! use tracer names stored in dynamics
          ! did not find old tracer name
          exit ! might as well stop here
        else
          ! found old tracer name
          count=count+1
        endif
      enddo
      if (count.eq.nq) then
        write(*,*) "physdem1:tracers seem to follow old naming ",
     &             "convention (qsurf01,qsurf02,...)"

        call abort
        !write(*,*) "   => will work for now ... "
        !write(*,*) "      but you should run newstart to rename them"
        !oldtracernames=.true.
        ! Moreover, if computing water cycle with ice, move surface ice
        ! back to qsurf(nq)
        !IF (iceparty) THEN
        !  write(*,*)'physdem1: moving surface water ice to index ',nq
        !  qsurf(1:ngrid,nq)=qsurf(1:ngrid,nq-1)
        !  qsurf(1:ngrid,nq-1)=0
        !ENDIF
      endif ! of if (count.eq.nq)

      IF(nq.GE.1) THEN
! preliminary stuff: look for water vapour & water ice tracers (if any)
        if (.not.oldtracernames) then
         do iq=1,nq
           if (nametrac(iq).eq."h2o_vap") then
             i_h2o_vap=iq
           endif
           if (nametrac(iq).eq."h2o_ice") then
             i_h2o_ice=iq
           endif
         enddo ! of iq=1,nq
         ! handle special case of only water vapour tracer (no ice)
         if ((i_h2o_vap.ne.0).and.(i_h2o_ice.eq.0)) then
          ! then the index of (surface) ice is i_h2o_vap
            print*,'water vapour but no water ice, WTF?'
            call abort
            i_h2o_ice=i_h2o_vap
         endif
        endif ! of if (.not.oldtracernames) 

         DO iq=1,nq
           IF (oldtracernames) THEN
             txt=" "
             write(txt,'(a5,i2.2)')'qsurf',iq
           ELSE
             txt=nametrac(iq)


!     in new version, h2o_vap acts as liquid surface tracer
!     so the section below is now redundant
!     ------------------------------------------------------------------
!             ! Exception: there is no water vapour surface tracer
!             if (txt.eq."h2o_vap") then
!               write(*,*)"physdem1: skipping water vapour tracer"
!               if (i_h2o_ice.eq.i_h2o_vap) then
!               ! then there is no "water ice" tracer; but still
!               ! there is some water ice on the surface
!                 write(*,*)"          writting water ice instead"
!                 txt="h2o_ice"
!               else
!               ! there is a "water ice" tracer which has been / will be
!               ! delt with in due time
!                 cycle
!               endif ! of if (igcm_h2o_ice.eq.igcm_h2o_vap)
!             endif ! of if (txt.eq."h2o_vap")
!     ------------------------------------------------------------------



           ENDIF ! of IF (oldtracernames)

           ierr=NF_REDEF(nid)
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to swich to NetCDF define mode'
             CALL abort
           ENDIF
#ifdef NC_DOUBLE
           ierr=NF_DEF_VAR(nid,txt,NF_DOUBLE,1,idim2,nvarid)
#else
           ierr=NF_DEF_VAR(nid,txt,NF_FLOAT,1,idim2,nvarid)
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to define ',trim(txt)
             CALL abort
           ENDIF
           ierr=NF_PUT_ATT_TEXT (nid, nvarid, "title", 17,
     &                        "tracer on surface")
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to define ',trim(txt),
     &               ' title attribute'
             CALL abort
           ENDIF
           ierr=NF_ENDDEF(nid)
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to swich out of define mode'
             CALL abort
           ENDIF
           
#ifdef NC_DOUBLE
            ierr=NF_PUT_VAR_DOUBLE (nid,nvarid,qsurf(1,iq))
#else
            ierr=NF_PUT_VAR_REAL (nid,nvarid,qsurf(1,iq))
#endif
           IF (ierr.NE.NF_NOERR) THEN
             PRINT*, 'physdem1: Failed to write ',trim(txt)
             CALL abort
           ENDIF
         ENDDO ! of DO iq=1,nq
      ENDIF ! of IF(nq.GE.1)

c close file
      ierr = NF_CLOSE(nid)

      RETURN

      END