source: trunk/LMDZ.COMMON/libf/evolution/writediagpem.F90 @ 3316

MODULE writediagpem_mod

implicit none

!=======================================================================
contains
!=======================================================================

SUBROUTINE writediagpem(ngrid,nom,titre,unite,dim,px)

!  Ecriture de variables diagnostiques au choix dans la physique
!  dans un fichier NetCDF nomme  'diagpem'. Ces variables peuvent etre
!  3d (ex : temperature), 2d (ex : temperature de surface), ou
!  0d (pour un scalaire qui ne depend que du temps : ex : la longitude
!  solaire)
!  (ou encore 1d, dans le cas de testphys1d, pour sortir une colonne)
!  La periode d'ecriture est donnee par
!  "ecritpem " regle dans le fichier de controle de run :  run.def
!
!    writediagpem peut etre appele de n'importe quelle subroutine
!    de la physique, plusieurs fois. L'initialisation et la creation du
!    fichier se fait au tout premier appel.
!
! WARNING : les variables dynamique (u,v,t,q,ps)
!  sauvees par writediagpem avec une
! date donnee sont legerement differentes que dans le fichier histoire car
! on ne leur a pas encore ajoute de la dissipation et de la physique !!!
! IL est  RECOMMANDE d'ajouter les tendance physique a ces variables
! avant l'ecriture dans diagpem (cf. physiq.F)
!
! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4
!         Oct 2011 Francois: enable having a 'diagpem.def' file to select
!                            at runtime, which variables to put in file
!         Oct 2023 JBC: conversion into Fortran 90 with module for the PEM
!
!  parametres (input) :
!  ----------
!      ngrid : nombres de point ou est calcule la physique
!                (ngrid = 2+(jjm-1)*iim - 1/jjm)
!                (= nlon ou klon dans la physique terrestre)
!      unit : unite logique du fichier de sortie (toujours la meme)
!      nom  : nom de la variable a sortir (chaine de caracteres)
!      titre: titre de la variable (chaine de caracteres)
!      unite : unite de la variable (chaine de caracteres)
!      px : variable a sortir (real 0, 1, 2, ou 3d)
!      dim : dimension de px : 0, 1, 2, ou 3 dimensions
!
!=================================================================
use surfdat_h,          only: phisfi
use geometry_mod,       only: cell_area
use mod_phys_lmdz_para, only: is_parallel, is_mpi_root, is_master, gather
use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat, nbp_lev, grid_type, unstructured
use time_evol_mod,      only: ecritpem, dt_pem

implicit none

! Commons
include "netcdf.inc"

! Arguments on input:
integer,                        intent(in) :: ngrid
character(len=*),               intent(in) :: nom, titre, unite
integer,                        intent(in) :: dim
real, dimension(ngrid,nbp_lev), intent(in) :: px

! Local variables:
real*4, dimension(nbp_lon + 1,nbp_lat,nbp_lev) :: dx3    ! to store a 3D data set
real*4, dimension(nbp_lon + 1,nbp_lat)         :: dx2    ! to store a 2D (surface) data set
real*4, dimension(nbp_lev)                     :: dx1    ! to store a 1D (column) data set
real*4                                         :: dx0
real*4, dimension(1,nbp_lev)                   :: dx3_1d ! to store a profile with 1D model
real*4                                         :: dx2_1d ! to store a surface value with 1D model
real*4, save                                   :: date

!$OMP THREADPRIVATE(date)
real, dimension((nbp_lon + 1),nbp_lat) :: phis
real, dimension((nbp_lon + 1),nbp_lat) :: area
integer                                :: irythme, ierr, ierr2, i, j, l, ig0
integer,                          save :: zitau = 0
character(27),                    save :: firstnom='1234567890'
!$OMP THREADPRIVATE(zitau,firstnom)

! Ajouts
integer,           save :: ntime = 0
!$OMP THREADPRIVATE(ntime)
integer                 :: idim, varid
integer                 :: nid
character(*), parameter :: fichnom = "diagpem.nc"
integer, dimension(4)   :: id
integer, dimension(4)   :: edges, corner

! Added to use diagpem.def to select output variable
logical,                                  save :: diagpem_def
logical                                        :: getout
integer,                                  save :: n_nom_def
integer                                        :: n
integer, parameter                             :: n_nom_def_max = 199
character(120), dimension(n_nom_def_max), save :: nom_def
logical,                                  save :: firstcall = .true.
!$OMP THREADPRIVATE(firstcall) !diagpem_def,n_nom_def,nom_def read in diagpem.def

#ifdef CPP_PARA
! Added to work in parallel mode
    real, dimension(klon_glo,nbp_lev)        :: dx3_glop
    real, dimension(nbp_lon,nbp_lat,nbp_lev) :: dx3_glo  ! to store a global 3D data set
    real, dimension(klon_glo)                :: dx2_glop
    real, dimension(nbp_lon,nbp_lat)         :: dx2_glo  ! to store a global 2D (surface) data set
    real, dimension(ngrid)                   :: px2
!    real, dimension(nbp_lev)                :: dx1_glo   ! to store a 1D (column) data set
!    real :: dx0_glo
    real, dimension(klon_glo) :: phisfi_glo ! surface geopotential on global physics grid
    real, dimension(klon_glo) :: areafi_glo ! mesh area on global physics grid
#else
    real, dimension(ngrid) :: phisfi_glo ! surface geopotential on global physics grid
    real, dimension(ngrid) :: areafi_glo ! mesh area on global physics grid
#endif

if (grid_type == unstructured) return

!***************************************************************
!Sortie des variables au rythme voulu
irythme = int(ecritpem) ! output rate set by ecritpem

!***************************************************************
! At very first call, check if there is a "diagpem.def" to use and read it
! ------------------------------------------------------------------------
IF (firstcall) THEN
    firstcall=.false.

!$OMP MASTER
    ! Open diagpem.def definition file if there is one:
        open(99,file="diagpem.def",status='old',form='formatted',iostat=ierr2)

        if (ierr2 == 0) then
            diagpem_def=.true.
            write(*,*) "*******************"
            write(*,*) "Reading diagpem.def"
            write(*,*) "*******************"
            do n=1,n_nom_def_max
                read(99,fmt='(a)',end=88) nom_def(n)
                write(*,*) 'Output in diagpem: ', trim(nom_def(n))
            enddo
 88         continue
            if (n.ge.n_nom_def_max) then
                write(*,*)"n_nom_def_max too small in writediagpem.F:",n
                call abort_physic("writediagpem","n_nom_def_max too small",1)
            end if
            n_nom_def=n-1
            close(99)
        else
            diagpem_def=.false.
        endif
!$OMP END MASTER
!$OMP BARRIER
ENDIF ! of IF (firstcall)

! Get out of write_diagpem if there is diagpem.def AND variable not listed
!  -----------------------------------------------------------------------
if (diagpem_def) then
    getout=.true.
    do n=1,n_nom_def
        if (trim(nom_def(n)).eq.nom) getout=.false.
    enddo
    if (getout) return
endif

! Initialisation of 'firstnom' and create/open the "diagpem.nc" NetCDF file
! -------------------------------------------------------------------------
! (at very first call to the subroutine, in accordance with diagpem.def)
if (firstnom.eq.'1234567890') then ! .true. for the very first valid
    ! call to this subroutine; now set 'firstnom'
    firstnom = nom
    ! just to be sure, check that firstnom is large enough to hold nom
    if (len_trim(firstnom).lt.len_trim(nom)) then
        write(*,*) "writediagpem: Error !!!"
        write(*,*) "   firstnom string not long enough!!"
        write(*,*) "   increase its size to at least ",len_trim(nom)
        call abort_physic("writediagpem","firstnom too short",1)
    endif

    zitau = -1 ! initialize zitau

#ifdef CPP_PARA
        ! Gather phisfi() geopotential on physics grid
        call Gather(phisfi,phisfi_glo)
        ! Gather cell_area() mesh area on physics grid
        call Gather(cell_area,areafi_glo)
#else
        phisfi_glo(:)=phisfi(:)
        areafi_glo(:)=cell_area(:)
#endif

    !! parallel: we cannot use the usual writediagpem method
    !! call iophys_ini
    if (is_master) then
        ! only the master is required to do this
        ! Create the NetCDF file
        ierr = NF_CREATE(fichnom,NF_CLOBBER,nid)
        ! Define the 'Time' dimension
        ierr = nf_def_dim(nid,"Time",NF_UNLIMITED,idim)
        ! Define the 'Time' variable
    !#ifdef NC_DOUBLE
    !    ierr = NF_DEF_VAR (nid, "Time", NF_DOUBLE, 1, idim,varid)
    !#else
        ierr = NF_DEF_VAR (nid, "Time", NF_FLOAT, 1, idim,varid)
    !#endif
        ! Add a long_name attribute
        ierr = NF_PUT_ATT_TEXT (nid, varid, "long_name",4,"Time")
        ! Add a units attribute
        ierr = NF_PUT_ATT_TEXT(nid, varid,'units',29,"days since 0000-00-0 00:00:00")
        ! Switch out of NetCDF Define mode
        ierr = NF_ENDDEF(nid)

        ! Build phis() and area()
        IF (klon_glo>1) THEN
            do i=1,nbp_lon+1 ! poles
            phis(i,1)=phisfi_glo(1)
            phis(i,nbp_lat)=phisfi_glo(klon_glo)
            ! for area, divide at the poles by nbp_lon
            area(i,1)=areafi_glo(1)/nbp_lon
            area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
            enddo
            do j=2,nbp_lat-1
            ig0= 1+(j-2)*nbp_lon
            do i=1,nbp_lon
                phis(i,j)=phisfi_glo(ig0+i)
                area(i,j)=areafi_glo(ig0+i)
            enddo
            ! handle redundant point in longitude
            phis(nbp_lon+1,j)=phis(1,j)
            area(nbp_lon+1,j)=area(1,j)
            enddo
        ENDIF

        ! write "header" of file (longitudes, latitudes, geopotential, ...)
        IF (klon_glo>1) THEN ! general 3D case
            call iniwrite(nid,0,phis,area,nbp_lon+1,nbp_lat)
        ELSE ! 1D model
            call iniwrite(nid,0,phisfi_glo(1),areafi_glo(1),1,1)
        ENDIF

    endif ! of if (is_master)

    else

    if (is_master) then
        ! only the master is required to do this
        ! Open the NetCDF file
        ierr = NF_OPEN(fichnom,NF_WRITE,nid)
    endif ! of if (is_master)

endif ! if (firstnom.eq.'1234567890')

! Increment time index 'zitau' if it is the "fist call" (at given time level)
! to writediagpem
!------------------------------------------------------------------------
if (nom == firstnom) zitau = zitau + dt_pem

!--------------------------------------------------------
! Write the variables to output file if it's time to do so
!--------------------------------------------------------
if (MOD(zitau+1,irythme) == 0.) then
! Compute/write/extend 'Time' coordinate (date given in days)
! (done every "first call" (at given time level) to writediagpem)
! Note: date is incremented as 1 step ahead of physics time
!--------------------------------------------------------
    if (is_master) then
        ! only the master is required to do this
        if (nom.eq.firstnom) then
        ! We have identified a "first call" (at given date)
            ntime=ntime+1 ! increment # of stored time steps
            ! compute corresponding date (in days and fractions thereof)
            date = float(zitau + 1)
            ! Get NetCDF ID of 'Time' variable
            ierr= NF_INQ_VARID(nid,"Time",varid)
            ! Write (append) the new date to the 'Time' array
!#ifdef NC_DOUBLE
!    ierr= NF_PUT_VARA_DOUBLE(nid,varid,[ntime],[1],[date])
!#else
            ierr= NF_PUT_VARA_REAL(nid,varid,[ntime],[1],[date])
!#endif
            if (ierr.ne.NF_NOERR) then
                write(*,*) "***** PUT_VAR matter in writediagpem_nc"
                write(*,*) "***** with time"
                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
!                call abort
            endif

            write(6,*)'WRITEDIAGPEM: date= ', date
        endif ! of if (nom.eq.firstnom)
    endif ! of if (is_master)

!Case of a 3D variable
!---------------------
    if (dim == 3) then
#ifdef CPP_PARA
        ! Gather field on a "global" (without redundant longitude) array
        call Gather(px,dx3_glop)
!$OMP MASTER
        if (is_mpi_root) then
            call Grid1Dto2D_glo(dx3_glop,dx3_glo)
            ! copy dx3_glo() to dx3(:) and add redundant longitude
            dx3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
            dx3(nbp_lon+1,:,:)=dx3(1,:,:)
        endif
!$OMP END MASTER
!$OMP BARRIER
#else
!         Passage variable physique -->  variable dynamique
!         recast (copy) variable from physics grid to dynamics grid
        IF (klon_glo>1) THEN ! General case
            DO l=1,nbp_lev
                DO i=1,nbp_lon+1
                    dx3(i,1,l)=px(1,l)
                    dx3(i,nbp_lat,l)=px(ngrid,l)
                ENDDO
                DO j=2,nbp_lat-1
                    ig0= 1+(j-2)*nbp_lon
                    DO i=1,nbp_lon
                        dx3(i,j,l)=px(ig0+i,l)
                    ENDDO
                    dx3(nbp_lon+1,j,l)=dx3(1,j,l)
                ENDDO
            ENDDO
        ELSE ! 1D model case
            dx3_1d(1,1:nbp_lev)=px(1,1:nbp_lev)
        ENDIF
#endif
!       Ecriture du champs
        if (is_master) then
! only the master writes to output
! name of the variable
            ierr= NF_INQ_VARID(nid,nom,varid)
            if (ierr /= NF_NOERR) then
! corresponding dimensions
                ierr = NF_INQ_DIMID(nid,"longitude",id(1))
                ierr = NF_INQ_DIMID(nid,"latitude",id(2))
                ierr = NF_INQ_DIMID(nid,"altitude",id(3))
                ierr = NF_INQ_DIMID(nid,"Time",id(4))

! Create the variable if it doesn't exist yet
                write (*,*) "==========================="
                write (*,*) "DIAGPEM: creating variable ",trim(nom)
                call def_var(nid,nom,titre,unite,4,id,varid,ierr)

            else
                if (ntime==0) then
                    write(*,*) "DIAGPEM Error: failed creating variable ",trim(nom)
                    write(*,*) "it seems it already exists!"
                    call abort_physic("writediagpem",trim(nom)//" already exists",1)
                endif
            endif

            corner(1)=1
            corner(2)=1
            corner(3)=1
            corner(4)=ntime

            IF (klon_glo==1) THEN
                edges(1)=1
            ELSE
                edges(1)=nbp_lon+1
            ENDIF
            edges(2)=nbp_lat
            edges(3)=nbp_lev
            edges(4)=1
!#ifdef NC_DOUBLE
!            ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx3)
!#else
!            write(*,*)"test:  nid=",nid," varid=",varid
!            write(*,*)"       corner()=",corner
!            write(*,*)"       edges()=",edges
!            write(*,*)"       dx3()=",dx3
            IF (klon_glo>1) THEN ! General case
                ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3)
            ELSE
                ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx3_1d)
            ENDIF
!#endif

            if (ierr.ne.NF_NOERR) then
                write(*,*) "***** PUT_VAR problem in writediagpem"
                write(*,*) "***** with dx3: ",trim(nom)
                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
                call abort_physic("writediagpem","failed writing "//trim(nom),1)
            endif

        endif !of if (is_master)

!Case of a 2D variable
!---------------------
    else if (dim == 2) then

#ifdef CPP_PARA
        ! Gather field on a "global" (without redundant longitude) array
        px2(:)=px(:,1)
        call Gather(px2,dx2_glop)
!$OMP MASTER
        if (is_mpi_root) then
            call Grid1Dto2D_glo(dx2_glop,dx2_glo)
            ! copy dx2_glo() to dx2(:) and add redundant longitude
            dx2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
            dx2(nbp_lon+1,:)=dx2(1,:)
        endif
!$OMP END MASTER
!$OMP BARRIER
#else
!         Passage variable physique -->  physique dynamique
!         recast (copy) variable from physics grid to dynamics grid
        IF (klon_glo>1) THEN ! General case
            DO i=1,nbp_lon+1
                dx2(i,1)=px(1,1)
                dx2(i,nbp_lat)=px(ngrid,1)
            ENDDO
            DO j=2,nbp_lat-1
                ig0= 1+(j-2)*nbp_lon
                DO i=1,nbp_lon
                    dx2(i,j)=px(ig0+i,1)
                ENDDO
                dx2(nbp_lon+1,j)=dx2(1,j)
            ENDDO
        ELSE ! 1D model case
            dx2_1d=px(1,1)
        ENDIF
#endif

        if (is_master) then
        ! only the master writes to output
!            write (*,*) 'In  writediagpem, on sauve:  ' , nom
!            write (*,*) 'In  writediagpem. Estimated date = ' ,date
            ierr= NF_INQ_VARID(nid,nom,varid)
            if (ierr /= NF_NOERR) then
! corresponding dimensions
                ierr= NF_INQ_DIMID(nid,"longitude",id(1))
                ierr= NF_INQ_DIMID(nid,"latitude",id(2))
                ierr= NF_INQ_DIMID(nid,"Time",id(3))

! Create the variable if it doesn't exist yet
                write (*,*) "==========================="
                write (*,*) "DIAGPEM: creating variable ",trim(nom)

                call def_var(nid,nom,titre,unite,3,id,varid,ierr)

            else
                if (ntime==0) then
                    write(*,*) "DIAGPEM Error: failed creating variable ",trim(nom)
                    write(*,*) "it seems it already exists!"
                    call abort_physic("writediagpem",trim(nom)//" already exists",1)
                endif
            endif

            corner(1)=1
            corner(2)=1
            corner(3)=ntime
            IF (klon_glo==1) THEN
                edges(1)=1
            ELSE
                edges(1)=nbp_lon+1
            ENDIF
            edges(2)=nbp_lat
            edges(3)=1


!#ifdef NC_DOUBLE
!           ierr = NF_PUT_VARA_DOUBLE (nid,varid,corner,edges,dx2)
!#else
            IF (klon_glo>1) THEN ! General case
                ierr = NF_PUT_VARA_REAL(nid,varid,corner,edges,dx2)
            ELSE
                ierr = NF_PUT_VARA_REAL(nid,varid,corner,edges,[dx2_1d])
            ENDIF
!#endif

            if (ierr.ne.NF_NOERR) then
                write(*,*) "***** PUT_VAR matter in writediagpem"
                write(*,*) "***** with dx2: ",trim(nom)
                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
                call abort_physic("writediagpem","failed writing "//trim(nom),1)
            endif

        endif !of if (is_master)

!Case of a 1D variable (ie: a column)
!---------------------------------------------------
    else if (dim.eq.1) then
        if (is_parallel) then
            write(*,*) "writediagpem error: dim=1 not implemented ","in parallel mode. Problem for ",trim(nom)
                call abort_physic("writediagpem","failed writing "//trim(nom),1)
        endif
!       Passage variable physique -->  physique dynamique
!       recast (copy) variable from physics grid to dynamics grid
        do l=1,nbp_lev
            dx1(l)=px(1,l)
        enddo

        ierr= NF_INQ_VARID(nid,nom,varid)
        if (ierr /= NF_NOERR) then
! corresponding dimensions
            ierr= NF_INQ_DIMID(nid,"altitude",id(1))
            ierr= NF_INQ_DIMID(nid,"Time",id(2))

! Create the variable if it doesn't exist yet

            write (*,*) "==========================="
            write (*,*) "DIAGPEM: creating variable ",trim(nom)

            call def_var(nid,nom,titre,unite,2,id,varid,ierr)

        else
            if (ntime==0) then
                write(*,*) "DIAGPEM Error: failed creating variable ",trim(nom)
                write(*,*) "it seems it already exists!"
                call abort_physic("writediagpem",trim(nom)//" already exists",1)
            endif
        endif

        corner(1)=1
        corner(2)=ntime
        edges(1)=nbp_lev
        edges(2)=1
!#ifdef NC_DOUBLE
!           ierr= NF_PUT_VARA_DOUBLE(nid,varid,corner,edges,dx1)
!#else
        ierr= NF_PUT_VARA_REAL(nid,varid,corner,edges,dx1)
!#endif

        if (ierr /= NF_NOERR) then
            write(*,*) "***** PUT_VAR problem in writediagpem"
            write(*,*) "***** with dx1: ",trim(nom)
            write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
            call abort_physic("writediagpem","failed writing "//trim(nom),1)
        endif

!Case of a 0D variable (ie: a time-dependent scalar)
!---------------------------------------------------
    else if (dim == 0) then

        dx0 = px (1,1)

        if (is_master) then
            ! only the master writes to output
            ierr= NF_INQ_VARID(nid,nom,varid)
            if (ierr /= NF_NOERR) then
! corresponding dimensions
                ierr= NF_INQ_DIMID(nid,"Time",id(1))

! Create the variable if it doesn't exist yet
                write (*,*) "==========================="
                write (*,*) "DIAGPEM: creating variable ",trim(nom)
                call def_var(nid,nom,titre,unite,1,id,varid,ierr)

            else
                if (ntime==0) then
                    write(*,*) "DIAGPEM Error: failed creating variable ",trim(nom)
                    write(*,*) "it seems it already exists!"
                    call abort_physic("writediagpem",trim(nom)//" already exists",1)
                endif
            endif

            corner(1)=ntime
            edges(1)=1

!#ifdef NC_DOUBLE
!            ierr = NF_PUT_VARA_DOUBLE (nid,varid,[corner(1)],[edges(1)],[dx0])
!#else
            ierr= NF_PUT_VARA_REAL(nid,varid,[corner(1)],[edges(1)],[dx0])
!#endif
            if (ierr.ne.NF_NOERR) then
                write(*,*) "***** PUT_VAR matter in writediagpem"
                write(*,*) "***** with dx0: ",trim(nom)
                write(*,*) 'ierr=', ierr,": ",NF_STRERROR(ierr)
                call abort_physic("writediagpem","failed writing "//trim(nom),1)
            endif
        endif !of if (is_master)
    endif ! of if (dim.eq.3) elseif(dim.eq.2)...
endif ! of if ( MOD(zitau+1,irythme) .eq.0.)

if (is_master) ierr= NF_CLOSE(nid)

END SUBROUTINE writediagpem

!=================================================================

SUBROUTINE writediagsoilpem(ngrid,name,title,units,dimpx,px)

! Write variable 'name' to NetCDF file 'diagsoil.nc'.
! The variable may be 3D (lon,lat,depth) subterranean field,
! a 2D (lon,lat) surface field, or a simple scalar (0D variable).
!
! Calls to 'writediagsoilpem' can originate from anywhere in the program;
! An initialisation of variable 'name' is done if it is the first time
! that this routine is called with given 'name'; otherwise data is appended
! (yielding the sought time series of the variable)

! Modifs: Aug.2010 Ehouarn: enforce outputs to be real*4

use comsoil_h_PEM,      only: nsoilmx_PEM, inertiedat_PEM
use geometry_mod,       only: cell_area
use mod_phys_lmdz_para, only: is_mpi_root, is_master, gather
use mod_grid_phy_lmdz,  only: klon_glo, Grid1Dto2D_glo, nbp_lon, nbp_lat
use mod_grid_phy_lmdz,  only: grid_type, unstructured
use time_evol_mod,      only: ecritpem, dt_pem

implicit none

include"netcdf.inc"

! Arguments:
integer,intent(in) :: ngrid ! number of (horizontal) points of physics grid
! i.e. ngrid = 2+(jjm-1)*iim - 1/jjm
character(len=*),intent(in) :: name ! 'name' of the variable
character(len=*),intent(in) :: title ! 'long_name' attribute of the variable
character(len=*),intent(in) :: units ! 'units' attribute of the variable
integer,intent(in) :: dimpx ! dimension of the variable (3,2 or 0)
real,dimension(ngrid,nsoilmx_PEM),intent(in) :: px ! variable

! Local variables:
real*4,dimension(nbp_lon+1,nbp_lat,nsoilmx_PEM) :: data3 ! to store 3D data
real*4,dimension(nbp_lon+1,nbp_lat) :: data2 ! to store 2D data
real*4 :: data0 ! to store 0D data
real*4 :: data3_1d(1,nsoilmx_PEM) ! to store a profile in 1D model
real*4 :: data2_1d ! to store surface value with 1D model
integer :: i,j,l ! for loops
integer :: ig0

real*4,save :: date ! time counter (in elapsed days)

real :: inertia((nbp_lon+1),nbp_lat,nsoilmx_PEM)
real :: area((nbp_lon+1),nbp_lat)

real :: inertiafi_glo(klon_glo,nsoilmx_PEM)
real :: areafi_glo(klon_glo)

integer,save :: isample ! sample rate at which data is to be written to output
integer,save :: ntime=0 ! counter to internally store time steps
character(len=20),save :: firstname="1234567890"
integer,save :: zitau=0
!$OMP THREADPRIVATE(date,isample,ntime,firstname,zitau)

character(len=30) :: filename="diagsoilpem.nc"

! NetCDF stuff:
integer :: nid ! NetCDF output file ID
integer :: varid ! NetCDF ID of a variable
integer :: ierr ! NetCDF routines return code
integer,dimension(4) :: id ! NetCDF IDs of the dimensions of the variable
integer,dimension(4) :: edges,corners

#ifdef CPP_PARA
! Added to work in parallel mode
real dx3_glop(klon_glo,nsoilmx_PEM)
real dx3_glo(nbp_lon,nbp_lat,nsoilmx_PEM) ! to store a global 3D data set
real dx2_glop(klon_glo)
real dx2_glo(nbp_lon,nbp_lat)     ! to store a global 2D (surface) data set
real px2(ngrid)
#endif

! 0. Do we ouput a diagsoil.nc file? If not just bail out now.

! additional check: one can only output diagsoil.nc files
! in lon-lat case (or 1D)
if (grid_type==unstructured) then
  write(*,*) "writediagsoil: Error !!!"
  write(*,*) "diagsoil.nc outputs not possible on unstructured grids!!"
  call abort_physic("writediagsoil","impossible on unstructured grid",1)
endif

! 1. Initialization step
if (firstname.eq."1234567890") then
  ! Store 'name' as 'firstname'
  firstname=name
  ! From now on, if 'name'.eq.'firstname', then it is a new time cycle

  ! just to be sure, check that firstnom is large enough to hold nom
  if (len_trim(firstname).lt.len_trim(name)) then
    write(*,*) "writediagsoil: Error !!!"
    write(*,*) "   firstname string not long enough!!"
    write(*,*) "   increase its size to at least ",len_trim(name)
    call abort_physic("writediagsoil","firstname too short",1)
  endif
  
  ! Set output sample rate
  isample = int(ecritpem) ! same as for diagpem outputs
  
  ! Create output NetCDF file
  if (is_master) then
   ierr=NF_CREATE(filename,NF_CLOBBER,nid)
   if (ierr.ne.NF_NOERR) then
    write(*,*)'writediagsoil: Error, failed creating file '//trim(filename)
    call abort_physic("writediagsoil","failed creating"//trim(filename),1)
   endif
  endif

#ifdef CPP_PARA
   ! Gather inertiedat() soil thermal inertia on physics grid
   call Gather(inertiedat_PEM,inertiafi_glo)
   ! Gather cell_area() mesh area on physics grid
   call Gather(cell_area,areafi_glo)
#else
         inertiafi_glo(:,:)=inertiedat_PEM(:,:)
         areafi_glo(:)=cell_area(:)
#endif

  if (is_master) then
   ! build inertia() and area()
   if (klon_glo>1) then
    do i=1,nbp_lon+1 ! poles
     inertia(i,1,1:nsoilmx_PEM)=inertiafi_glo(1,1:nsoilmx_PEM)
     inertia(i,nbp_lat,1:nsoilmx_PEM)=inertiafi_glo(klon_glo,1:nsoilmx_PEM)
     ! for area, divide at the poles by nbp_lon
     area(i,1)=areafi_glo(1)/nbp_lon
     area(i,nbp_lat)=areafi_glo(klon_glo)/nbp_lon
    enddo
    do j=2,nbp_lat-1
     ig0= 1+(j-2)*nbp_lon
     do i=1,nbp_lon
        inertia(i,j,1:nsoilmx_PEM)=inertiafi_glo(ig0+i,1:nsoilmx_PEM)
        area(i,j)=areafi_glo(ig0+i)
     enddo
     ! handle redundant point in longitude
     inertia(nbp_lon+1,j,1:nsoilmx_PEM)=inertia(1,j,1:nsoilmx_PEM)
     area(nbp_lon+1,j)=area(1,j)
    enddo
   endif
   
   ! write "header" of file (longitudes, latitudes, geopotential, ...)
   if (klon_glo>1) then ! general 3D case
     call iniwritesoil(nid,ngrid,inertia,area,nbp_lon+1,nbp_lat)
   else ! 1D model
     call iniwritesoil(nid,ngrid,inertiafi_glo(1,:),areafi_glo(1),1,1)
   endif

  endif ! of if (is_master)
  
  ! set zitau to -1 to be compatible with zitau incrementation step below
  zitau=-1
  
else
  ! If not an initialization call, simply open the NetCDF file
  if (is_master) then
   ierr=NF_OPEN(filename,NF_WRITE,nid)
  endif
endif ! of if (firstname.eq."1234567890")

! 2. Increment local time counter, if necessary
if (name.eq.firstname) then
  ! if we run across 'firstname', then it is a new time step
  zitau = zitau + dt_pem
endif

! 3. Write data, if the time index matches the sample rate
if (mod(zitau+1,isample).eq.0) then

! 3.1 If first call at this date, update 'time' variable
  if (name.eq.firstname) then
    ntime=ntime+1
    date = float(zitau + 1)
    
    if (is_master) then
     ! Get NetCDF ID for "time"
     ierr=NF_INQ_VARID(nid,"time",varid)
     ! Add the current value of date to the "time" array
!#ifdef NC_DOUBLE
!     ierr=NF_PUT_VARA_DOUBLE(nid,varid,ntime,1,date)
!#else
     ierr=NF_PUT_VARA_REAL(nid,varid,[ntime],[1],[date])
!#endif
     if (ierr.ne.NF_NOERR) then
      write(*,*)"writediagsoil: Failed writing date to time variable"
      call abort_physic("writediagsoil","failed writing time",1) 
     endif
    endif ! of if (is_master)
  endif ! of if (name.eq.firstname)

! 3.2 Write the variable to the NetCDF file
if (dimpx.eq.3) then ! Case of a 3D variable
  ! A. Recast data along 'dynamics' grid
#ifdef CPP_PARA
  ! gather field on a "global" (without redundant longitude) array
  call Gather(px,dx3_glop)
!$OMP MASTER
  if (is_mpi_root) then
    call Grid1Dto2D_glo(dx3_glop,dx3_glo)
    ! copy dx3_glo() to dx3(:) and add redundant longitude
    data3(1:nbp_lon,:,:)=dx3_glo(1:nbp_lon,:,:)
    data3(nbp_lon+1,:,:)=data3(1,:,:)
  endif
!$OMP END MASTER
!$OMP BARRIER
#else
  if (klon_glo>1) then ! General case
   do l=1,nsoilmx_PEM
    ! handle the poles
    do i=1,nbp_lon+1
      data3(i,1,l)=px(1,l)
      data3(i,nbp_lat,l)=px(ngrid,l)
    enddo
    ! rest of the grid
    do j=2,nbp_lat-1
      ig0=1+(j-2)*nbp_lon
      do i=1,nbp_lon
        data3(i,j,l)=px(ig0+i,l)
      enddo
      data3(nbp_lon+1,j,l)=data3(1,j,l) ! extra (modulo) longitude
    enddo
   enddo
  else ! 1D model case
   data3_1d(1,1:nsoilmx_PEM)=px(1,1:nsoilmx_PEM)
  endif
#endif
  
  ! B. Write (append) the variable to the NetCDF file
  if (is_master) then
  ! B.1. Get the ID of the variable
  ierr=NF_INQ_VARID(nid,name,varid)
  if (ierr.ne.NF_NOERR) then
    ! If we failed geting the variable's ID, we assume it is because
    ! the variable doesn't exist yet and must be created.
    ! Start by obtaining corresponding dimensions IDs
    ierr=NF_INQ_DIMID(nid,"longitude",id(1))
    ierr=NF_INQ_DIMID(nid,"latitude",id(2))
    ierr=NF_INQ_DIMID(nid,"depth",id(3))
    ierr=NF_INQ_DIMID(nid,"time",id(4))
    ! Tell the world about it
    write(*,*) "====================="
    write(*,*) "writediagsoil: creating variable "//trim(name)
    call def_var(nid,name,title,units,4,id,varid,ierr)
  endif ! of if (ierr.ne.NF_NOERR)
  
  ! B.2. Prepare things to be able to write/append the variable
  corners(1)=1
  corners(2)=1
  corners(3)=1
  corners(4)=ntime
  
  if (klon_glo==1) then
    edges(1)=1
  else
    edges(1)=nbp_lon+1
  endif
  edges(2)=nbp_lat
  edges(3)=nsoilmx_PEM
  edges(4)=1
  
  ! B.3. Write the slab of data
!#ifdef NC_DOUBLE
!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data3)
!#else
  if (klon_glo>1) then
    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3)
  else
    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data3_1d)
  endif
!#endif
  if (ierr.ne.NF_NOERR) then
    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
               " to file "//trim(filename)//" at time",date
  endif
  endif ! of if (is_master)

elseif (dimpx.eq.2) then ! Case of a 2D variable

  ! A. Recast data along 'dynamics' grid
#ifdef CPP_PARA
  ! gather field on a "global" (without redundant longitude) array
  px2(:)=px(:,1)
  call Gather(px2,dx2_glop)
!$OMP MASTER
  if (is_mpi_root) then
    call Grid1Dto2D_glo(dx2_glop,dx2_glo)
    ! copy dx3_glo() to dx3(:) and add redundant longitude
    data2(1:nbp_lon,:)=dx2_glo(1:nbp_lon,:)
    data2(nbp_lon+1,:)=data2(1,:)
  endif
!$OMP END MASTER
!$OMP BARRIER
#else
  if (klon_glo>1) then ! general case
    ! handle the poles
    do i=1,nbp_lon+1
      data2(i,1)=px(1,1)
      data2(i,nbp_lat)=px(ngrid,1)
    enddo
    ! rest of the grid
    do j=2,nbp_lat-1
      ig0=1+(j-2)*nbp_lon
      do i=1,nbp_lon
        data2(i,j)=px(ig0+i,1)
      enddo
      data2(nbp_lon+1,j)=data2(1,j) ! extra (modulo) longitude
    enddo
  else ! 1D model case
    data2_1d=px(1,1)
  endif
#endif

  ! B. Write (append) the variable to the NetCDF file
  if (is_master) then
  ! B.1. Get the ID of the variable
  ierr=NF_INQ_VARID(nid,name,varid)
  if (ierr.ne.NF_NOERR) then
    ! If we failed geting the variable's ID, we assume it is because
    ! the variable doesn't exist yet and must be created.
    ! Start by obtaining corresponding dimensions IDs
    ierr=NF_INQ_DIMID(nid,"longitude",id(1))
    ierr=NF_INQ_DIMID(nid,"latitude",id(2))
    ierr=NF_INQ_DIMID(nid,"time",id(3))
    ! Tell the world about it
    write(*,*) "====================="
    write(*,*) "writediagsoil: creating variable "//trim(name)
    call def_var(nid,name,title,units,3,id,varid,ierr)
  endif ! of if (ierr.ne.NF_NOERR)

  ! B.2. Prepare things to be able to write/append the variable
  corners(1)=1
  corners(2)=1
  corners(3)=ntime
  
  if (klon_glo==1) then
    edges(1)=1
  else
    edges(1)=nbp_lon+1
  endif
  edges(2)=nbp_lat
  edges(3)=1
  
  ! B.3. Write the slab of data
!#ifdef NC_DOUBLE
!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data2)
!#else
  if (klon_glo>1) then ! General case
    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,data2)
  else
    ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,[data2_1d])
  endif
!#endif
  if (ierr.ne.NF_NOERR) then
    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
               " to file "//trim(filename)//" at time",date
  endif
  endif ! of if (is_master)

elseif (dimpx.eq.0) then ! Case of a 0D variable
#ifdef CPP_PARA
  write(*,*) "writediagsoil: dimps==0 case not implemented in // mode!!"
  call abort_physic("writediagsoil","dimps==0 not implemented",1)
#endif
  ! A. Copy data value
  data0=px(1,1)

  ! B. Write (append) the variable to the NetCDF file
  ! B.1. Get the ID of the variable
  ierr=NF_INQ_VARID(nid,name,varid)
  if (ierr.ne.NF_NOERR) then
    ! If we failed geting the variable's ID, we assume it is because
    ! the variable doesn't exist yet and must be created.
    ! Start by obtaining corresponding dimensions IDs
    ierr=NF_INQ_DIMID(nid,"time",id(1))
    ! Tell the world about it
    write(*,*) "====================="
    write(*,*) "writediagsoil: creating variable "//trim(name)
    call def_var(nid,name,title,units,1,id,varid,ierr)
  endif ! of if (ierr.ne.NF_NOERR)

  ! B.2. Prepare things to be able to write/append the variable
  corners(1)=ntime
  
  edges(1)=1

  ! B.3. Write the data
!#ifdef NC_DOUBLE
!  ierr=NF_PUT_VARA_DOUBLE(nid,varid,corners,edges,data0)
!#else
  ierr=NF_PUT_VARA_REAL(nid,varid,corners,edges,[data0])
!#endif
  if (ierr.ne.NF_NOERR) then
    write(*,*) "writediagsoil: Error: Failed writing "//trim(name)//&
               " to file "//trim(filename)//" at time",date
  endif

endif ! of if (dimpx.eq.3) elseif (dimpx.eq.2) ...
endif ! of if (mod(zitau+1,isample).eq.0)

! 4. Close the NetCDF file
if (is_master) then
  ierr=NF_CLOSE(nid)
endif

END SUBROUTINE writediagsoilpem


END MODULE writediagpem_mod