write_lfifm1_for_diag.f90

!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
!################################
MODULE MODI_WRITE_LFIFM1_FOR_DIAG
!################################
INTERFACE
      SUBROUTINE WRITE_LFIFM1_FOR_DIAG(TPFILE,HDADFILE)
!
USE MODD_IO_ll, ONLY : TFILEDATA
!
TYPE(TFILEDATA),   INTENT(IN) :: TPFILE       ! outpput data file
CHARACTER(LEN=28), INTENT(IN) :: HDADFILE     ! corresponding FM-file name of 
                                              ! its DAD model
!
END SUBROUTINE WRITE_LFIFM1_FOR_DIAG
END INTERFACE
END MODULE MODI_WRITE_LFIFM1_FOR_DIAG
!
!     ##################################################
      SUBROUTINE WRITE_LFIFM1_FOR_DIAG(TPFILE,HDADFILE)
!     ##################################################
!
!!****  *WRITE_LFIFM1* - routine to write a LFIFM file for model 1
!!
!!    PURPOSE
!!    -------
!        The purpose of this routine is to write an initial LFIFM File 
!     of name YFMFILE2//'.lfi' with the FM routines.  
!
!!**  METHOD
!!    ------
!!      The data are written in the LFIFM file :
!!        - dimensions
!!        - grid variables
!!        - configuration variables
!!        - prognostic variables at time t and t-dt
!!        - 1D anelastic reference state
!!
!!      The localization on the model grid is also indicated :
!!
!!        IGRID = 1 for mass grid point
!!        IGRID = 2 for U grid point
!!        IGRID = 3 for V grid point
!!        IGRID = 4 for w grid point
!!        IGRID = 0 for meaningless case
!!          
!!
!!    EXTERNAL
!!    --------
!!      FMWRIT : FM-routine to write a record
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!      Module MODD_DIM1   : contains dimensions
!!      Module MODD_TIME1   : contains time variables and uses MODD_TIME
!!      Module MODD_GRID    : contains spatial grid variables for all models
!!      Module MODD_GRID1 : contains spatial grid variables
!!      Module MODD_REF     : contains reference state variables
!!      Module MODD_LUNIT1: contains logical unit variables.
!!      Module MODD_CONF    : contains configuration variables for all models
!!      Module MODD_CONF1  : contains configuration variables
!!      Module MODD_FIELD1  : contains prognostic variables
!!      Module MODD_GR_FIELD1 : contains surface prognostic variables
!!      Module MODD_LSFIELD1  : contains Larger Scale variables
!!      Module MODD_PARAM1    : contains parameterization options
!!      Module MODD_TURB1    : contains turbulence options
!!      Module MODD_FRC    : contains forcing variables
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!  	V. Ducrocq   *Meteo France* 
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    06/05/94 
!!       V. Ducrocq    27/06/94
!!       J.Stein       20/10/94 (name of the FMFILE)
!!       J.Stein       06/12/94 add the LS fields
!!       J.P. Lafore   09/01/95 add the DRYMASST
!!       J.Stein       20/01/95 add TKE and change the ycomment for the water
!!                              variables
!!       J.Stein       23/01/95 add a TKE switch and MODD_PARAM1
!!       J.Stein       16/03/95 remove R from the historical variables
!!       J.Stein       20/03/95 add the EPS var.
!!       J.Stein       30/06/95 add the variables related to the subgrid condens
!!       S. Belair     01/09/95 add surface variables and ground parameters
!!       J.-P. Pinty   15/09/95 add the radiation parameters
!!       J.Stein       23/01/96 add the TSZ0 option for the surface scheme
!!       M.Georgelin   13/12/95 add the forcing variables
!!       J.-P. Pinty   15/02/96 add external control for the forcing
!!       J.Stein P.Bougeault  15/03/96 add the cloud fraction and change the
!!                                     surface parameters for TSZ0 option
!!       J.Stein P.Jabouille  30/04/96 add the storage type
!!       J.Stein P.Jabouille  20/05/96 switch for XSIGS and XSRC
!!       J.Stein              10/10/96 change Xsrc into XSRCM and XRCT
!!       J.P. Lafore          30/07/96 add YFMFILE2 and HDADFILE writing
!!                                     corresponding to MY_NAME and DAD_NAME (for nesting)
!!       V.Masson             08/10/96 add LTHINSHELL
!!       J.-P. Pinty   15/12/96 add the microphysics (ice)
!!       J.-P. Pinty   11/01/97 add the deep convection
!!       J.-P. Pinty   27/01/97 split the recording of the SV array
!!       J.-P. Pinty   29/01/97 set recording of PRCONV and PACCONV in mm/h and
!!                                                         mm respectively
!!       J. Viviand    04/02/97 convert precipitation rates in mm/h
!!       P. Hereil     04/12/97 add the calculation of cloud top and moist PV
!!       P.Hereil N Asencio 3/02/98 add the calculation of  precipitation on large scale grid mesh
!!       N Asencio 2/10/98 suppress flux calculation if start file
!!       V Masson 25/11/98 places dummy arguments in module MODD_DIAG_FLAG
!!       V Masson 04/01/00 removes TSZ0 option
!!       J.-P. Pinty   29/11/02 add C3R5, ICE2, ICE4, CELEC
!!       V Masson 01/2004  removes surface (externalization)
!!       P. Tulet 01/2005   add dust, orilam
!!       M. Leriche 04/2007 add aqueous concentration in M
!!       O. Caumont 03/2008 add simulation of radar observations
!!       O. Caumont 14/09/2009 modifications to allow for polar outputs (radar diagnostics)
!!       October 2009 (G. Tanguy) add ILENCH=LEN(YCOMMENT) after
!!                                              change of YCOMMENT
!!       G. Tanguy  10/2009 add possibility to run radar after 
!!                          PREP_REAL_CASE with AROME
!!       O. Caumont 01/2011 [radar diagnostics] add control check for NMAX; revise comments
!!       O. Caumont 05/2011 [radar diagnostics] change output format
!!       G.Tanguy/ JP Pinty/ JP Chabureau 18/05/2011 : add lidar simulator
!!       S.Bielli 12/2012 : add latitude and longitude
!!       F. Duffourg 02/2013 : add new fields
!!      J.Escobar 21/03/2013: for HALOK get correctly local array dim/bound
!!       J. escobar 27/03/2014 : write LAT/LON only in not CARTESIAN case
!!       G.Delautier    2014 : remplace MODD_RAIN_C2R2_PARAM par MODD_RAIN_C2R2_KHKO_PARAM
!!       C. Augros 2014 : new radar simulator (T matrice)
!!       D.Ricard 2015 : add THETAES + POVOES  (LMOIST_ES=T)
!!      Modification    01/2016  (JP Pinty) Add LIMA
!!       C.Lac  04/2016 : add visibility and droplet deposition
!! 10/2017      (G.Delautier) New boundary layer height : replace LBLTOP by CBLTOP 
!!       T.Dauhut      10/2017 : add parallel 3D clustering
!!  Philippe Wautelet: 05/2016-04/2018: new data structures and calls for I/O
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE MODD_DIM_n
USE MODD_CONF
USE MODD_CONF_n
USE MODD_GRID
USE MODD_GRID_n
USE MODD_IO_ll, ONLY : TFILEDATA
USE MODD_METRICS_n
USE MODD_TIME
USE MODD_TIME_n
USE MODD_DYN_n
USE MODD_FIELD_n
USE MODD_GR_FIELD_n
USE MODD_LSFIELD_n
USE MODD_PARAM_n
USE MODD_CURVCOR_n
USE MODD_REF
USE MODD_REF_n
USE MODD_LUNIT, ONLY : TLUOUT0
USE MODD_LUNIT_n
USE MODD_TURB_n
USE MODD_RADIATIONS_n
USE MODD_FRC
USE MODD_PRECIP_n
USE MODD_CST
USE MODD_CLOUDPAR
USE MODD_DEEP_CONVECTION_n
USE MODD_PARAM_KAFR_n
USE MODD_NESTING
USE MODD_PARAMETERS
USE MODD_DIAG_FLAG
USE MODD_NSV
USE MODD_CH_M9_n,         ONLY : CNAMES, NEQAQ
USE MODD_RAIN_C2R2_DESCR, ONLY : C2R2NAMES
USE MODD_ICE_C1R3_DESCR,  ONLY : C1R3NAMES
USE MODD_ELEC_DESCR, ONLY : CELECNAMES
USE MODD_RAIN_C2R2_KHKO_PARAM
USE MODD_ICE_C1R3_PARAM
USE MODD_PARAM_ICE,       ONLY : LSEDIC
USE MODD_PARAM_LIMA,      ONLY : NMOD_CCN, NMOD_IFN, NMOD_IMM, NINDICE_CCN_IMM,&
                                 LSCAV, LHHONI, LAERO_MASS,                    &
                                 LLIMA_DIAG,                                   &
                                 NSPECIE, XMDIAM_IFN, XSIGMA_IFN, ZFRAC=>XFRAC,&
                                 XR_MEAN_CCN, XLOGSIG_CCN 
USE MODD_PARAM_LIMA_WARM, ONLY : CLIMA_WARM_CONC, CAERO_MASS
USE MODD_PARAM_LIMA_COLD, ONLY : CLIMA_COLD_CONC
USE MODD_LG,              ONLY : CLGNAMES
USE MODD_PASPOL,          ONLY : LPASPOL
USE MODD_CONDSAMP,        ONLY : LCONDSAMP
!
USE MODD_DIAG_FLAG
USE MODD_RADAR, ONLY: XLAT_RAD,XELEV,&
     XSTEP_RAD,NBRAD,NBELEV,NBAZIM,NBSTEPMAX,&
     NCURV_INTERPOL,LATT,LCART_RAD,NPTS_H,NPTS_V,XGRID,&
     LREFR,LDNDZ,NMAX,CNAME_RAD,NDIFF,&
     XLON_RAD,XALT_RAD,XLAM_RAD,XDT_RAD,LWBSCS,LWREFL
!
USE MODI_RADAR_SIMULATOR
!
USE MODD_DUST
USE MODD_CSTS_DUST
USE MODD_SALT
USE MODD_BLOWSNOW
USE MODD_CH_AEROSOL
USE MODD_CH_AERO_n
USE MODD_CH_MNHC_n
USE MODE_DUST_PSD
USE MODE_SALT_PSD
USE MODE_BLOWSNOW_PSD
USE MODE_AERO_PSD
USE MODI_GRADIENT_M
USE MODI_GRADIENT_W
USE MODI_GRADIENT_U
USE MODI_GRADIENT_V
USE MODI_SHUMAN
USE MODI_RADAR_RAIN_ICE
USE MODI_INI_RADAR
USE MODI_COMPUTE_MEAN_PRECIP
USE MODI_UV_TO_ZONAL_AND_MERID
USE MODI_CALCSOUND
USE MODI_FREE_ATM_PROFILE
USE MODI_GPS_ZENITH
USE MODI_CONTRAV
!
USE MODE_FM,              ONLY: IO_FILE_CLOSE_ll,IO_FILE_OPEN_ll
USE MODE_GRIDPROJ
USE MODE_FIELD
USE MODE_FMWRIT
USE MODE_GATHER_ll
USE MODE_ll
USE MODE_IO_ll
USE MODE_IO_MANAGE_STRUCT,ONLY: IO_FILE_ADD2LIST
USE MODE_THERMO
USE MODE_MODELN_HANDLER
USE MODI_LIDAR
USE MODI_CLUSTERING
!
USE MODD_MPIF
USE MODD_VAR_ll
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!
TYPE(TFILEDATA),   INTENT(IN) :: TPFILE       ! outpput data file
CHARACTER(LEN=28), INTENT(IN) :: HDADFILE     ! corresponding FM-file name of 
                                              ! its DAD model
!
!*       0.2   Declarations of local variables
!
INTEGER           :: IRESP          ! return-code for the file routines 
!
CHARACTER(LEN=3)  :: YFRC           ! to mark the time of the forcing
CHARACTER(LEN=31) :: YFGRI          ! file name for GPS stations
!
INTEGER           :: IIU,IJU,IKU,IIB,IJB,IKB,IIE,IJE,IKE ! Arrays bounds
! 
INTEGER                :: JLOOP,JI,JJ,JK,JSV,JT,JH,JV,JEL    ! loop index
INTEGER :: IMI ! Current model index
! 
REAL :: ZRV_OV_RD !  XRV / XRD
REAL :: ZGAMREF   ! Standard atmosphere lapse rate (K/m)
REAL :: ZX0D      ! work real scalar  
REAL :: ZLATOR, ZLONOR ! geographical coordinates of 1st mass point
REAL :: ZXHATM, ZYHATM ! conformal    coordinates of 1st mass point
REAL, DIMENSION(:), ALLOCATABLE   :: ZXHAT_ll    !  Position x in the conformal
                                                 ! plane (array on the complete domain)
REAL, DIMENSION(:), ALLOCATABLE   :: ZYHAT_ll    !   Position y in the conformal
                                                 ! plane (array on the complete domain)
!
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3))  :: ZPOVO
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3))  :: ZTEMP
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3))  :: ZVOX,ZVOY,ZVOZ 
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3))  :: ZCORIOZ 
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3))  :: ZWORK31,ZWORK32
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3))  :: ZWORK33,ZWORK34
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2))               :: ZWORK21,ZWORK22
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2))               :: ZWORK23,ZWORK24
REAL,DIMENSION(:,:,:,:,:), ALLOCATABLE                  :: ZWORK42 ! reflectivity on a cartesian grid (PREFL_CART)
REAL,DIMENSION(:,:,:,:,:), ALLOCATABLE                  :: ZWORK42_BIS
REAL,DIMENSION(:,:,:), ALLOCATABLE                      :: ZWORK43 ! latlon coordinates of cartesian grid points (PLATLON)
REAL,DIMENSION(:,:,:), ALLOCATABLE                      :: ZPHI,ZTHETAE,ZTHETAV
REAL,DIMENSION(:,:,:), ALLOCATABLE                      :: ZVISIKUN,ZVISIGUL,ZVISIZHA 
REAL,DIMENSION(:,:,:), ALLOCATABLE                      :: ZTHETAES
INTEGER, DIMENSION(:,:), ALLOCATABLE                    :: IWORK1
integer :: ICURR,INBOUT,IERR
!
REAL,DIMENSION(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),NSP+NCARB+NSOA,JPMODE):: ZPTOTA
REAL,DIMENSION(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),NMODE_DST*2):: ZSDSTDEP
REAL,DIMENSION(SIZE(XSVT,1),SIZE(XSVT,2),SIZE(XSVT,3),NMODE_SLT*2):: ZSSLTDEP
REAL,DIMENSION(:,:,:,:), ALLOCATABLE  :: ZSIG_DST, ZRG_DST, ZN0_DST
REAL,DIMENSION(:,:,:,:), ALLOCATABLE  :: ZSIG_SLT, ZRG_SLT, ZN0_SLT
REAL,DIMENSION(:,:,:), ALLOCATABLE    :: ZBET_SNW, ZRG_SNW
REAL,DIMENSION(:,:,:,:), ALLOCATABLE  :: ZMA_SNW
REAL,DIMENSION(:,:,:), ALLOCATABLE  :: ZRHOT, ZTMP ! work array
!
! GBOTUP = True does clustering from bottom up to top, False top down to surface
LOGICAL                                                   :: GBOTUP ! clustering propagation
LOGICAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3)) :: GCLOUD ! mask
INTEGER,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3)) :: ICLUSTERID, ICLUSTERLV
REAL,   DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3)) :: ZCLDSIZE

!ECRITURE DANS UN FICHIER ASCII DE RESULTATS 
!INITIALISATION DU NOM DE FICHIER CREE EN PARALLELE AVEC CELUI LFI
TYPE(TFILEDATA),POINTER :: TZRSFILE
INTEGER :: ILURS
CHARACTER(LEN=32) :: YRS
CHARACTER(LEN=3),DIMENSION(:),ALLOCATABLE  :: YRAD
CHARACTER(LEN=2*INT(NBSTEPMAX*XSTEP_RAD/XGRID)*2*9+1), DIMENSION(:), ALLOCATABLE :: CLATLON
CHARACTER(LEN=2*9) :: CBUFFER
CHARACTER(LEN=4)  :: YELEV
CHARACTER(LEN=3)  :: YGRID_SIZE
INTEGER :: IEL,IIELV
CHARACTER(LEN=5)  :: YVIEW   ! Upward or Downward integration
INTEGER           :: IACCMODE
!
!-------------------------------------------------------------------------------
INTEGER :: IAUX ! work variable 
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2),SIZE(XTHT,3)) :: ZWORK35,ZWORK36
REAL,DIMENSION(SIZE(XTHT,1),SIZE(XTHT,2))              :: ZWORK25,ZWORK26
REAL    :: ZEAU ! Mean precipitable water
INTEGER, DIMENSION(SIZE(XZZ,1),SIZE(XZZ,2))          ::IKTOP ! level in which is the altitude 3000m
REAL, DIMENSION(SIZE(XZZ,1),SIZE(XZZ,2),SIZE(XZZ,3)) :: ZDELTAZ ! interval (m) between two levels K
INTEGER :: ILUOUT0 ! Logical unit number for output-listing
!
CHARACTER(LEN=2)  :: INDICE
CHARACTER(LEN=100) :: YMSG
INTEGER           :: I,IID
TYPE(TFIELDDATA)              :: TZFIELD
TYPE(TFIELDDATA),DIMENSION(2) :: TZFIELD2
!
! LIMA LIDAR
REAL,DIMENSION(:,:,:,:), ALLOCATABLE :: ZTMP1, ZTMP2, ZTMP3, ZTMP4
!
! hauteur couche limite
REAL,DIMENSION(:,:,:),ALLOCATABLE :: ZZZ_GRID1
REAL,DIMENSION(:,:),ALLOCATABLE :: ZTHVSOL,ZSHMIX
REAL,DIMENSION(:,:,:),ALLOCATABLE :: ZZONWIND,ZMERWIND,ZFFWIND2,ZRIB
!
!-------------------------------------------------------------------------------
!
!*       0.     ARRAYS BOUNDS INITIALIZATION
!
CALL GET_DIM_EXT_ll ('B',IIU,IJU)
CALL GET_INDICE_ll (IIB,IJB,IIE,IJE)
IKU=NKMAX+2*JPVEXT
IKB=1+JPVEXT
IKE=IKU-JPVEXT

IMI = GET_CURRENT_MODEL_INDEX()
ILUOUT0 = TLUOUT0%NLU
TZRSFILE => NULL()
!-------------------------------------------------------------------------------
!
!*       1.     WRITES IN THE LFI FILE
!               ---------------------- 
!
!*       1.0    TPFILE%CNAME and HDADFILE :
!
CALL IO_WRITE_FIELD(TPFILE,'MASDEV',   NMASDEV)
CALL IO_WRITE_FIELD(TPFILE,'BUGFIX',   NBUGFIX)
CALL IO_WRITE_FIELD(TPFILE,'BIBUSER',  CBIBUSER)
CALL IO_WRITE_FIELD(TPFILE,'PROGRAM',  CPROGRAM)
!
CALL IO_WRITE_FIELD(TPFILE,'L1D',      L1D)
CALL IO_WRITE_FIELD(TPFILE,'L2D',      L2D)
CALL IO_WRITE_FIELD(TPFILE,'PACK',     LPACK)
!
CALL IO_WRITE_FIELD(TPFILE,'MY_NAME',  TPFILE%CNAME)
CALL IO_WRITE_FIELD(TPFILE,'DAD_NAME', HDADFILE)
!
IF (LEN_TRIM(HDADFILE)>0) THEN
  CALL IO_WRITE_FIELD(TPFILE,'DXRATIO',NDXRATIO_ALL(1))
  CALL IO_WRITE_FIELD(TPFILE,'DYRATIO',NDYRATIO_ALL(1))
  CALL IO_WRITE_FIELD(TPFILE,'XOR',    NXOR_ALL(1))
  CALL IO_WRITE_FIELD(TPFILE,'YOR',    NYOR_ALL(1))
END IF
!
CALL IO_WRITE_FIELD(TPFILE,'SURF',     CSURF)
!
!*       1.1    Type and Dimensions :
!
CALL IO_WRITE_FIELD(TPFILE,'STORAGE_TYPE','DI')
!
CALL IO_WRITE_FIELD(TPFILE,'IMAX',NIMAX_ll)
CALL IO_WRITE_FIELD(TPFILE,'JMAX',NJMAX_ll)
CALL IO_WRITE_FIELD(TPFILE,'KMAX',NKMAX)
!
CALL IO_WRITE_FIELD(TPFILE,'JPHEXT',JPHEXT)
!
!*       1.2    Grid variables :
!
IF (.NOT.LCARTESIAN) THEN
  CALL IO_WRITE_FIELD(TPFILE,'RPK',   XRPK)
  CALL IO_WRITE_FIELD(TPFILE,'LONORI',XLONORI)
  CALL IO_WRITE_FIELD(TPFILE,'LATORI',XLATORI)
! 
!* diagnostic of 1st mass point
!
  ALLOCATE(ZXHAT_ll(NIMAX_ll+ 2 * JPHEXT),ZYHAT_ll(NJMAX_ll+2 * JPHEXT))
  CALL GATHERALL_FIELD_ll('XX',XXHAT,ZXHAT_ll,IRESP) !//
  CALL GATHERALL_FIELD_ll('YY',XYHAT,ZYHAT_ll,IRESP) !//
  ZXHATM = 0.5 * (ZXHAT_ll(1)+ZXHAT_ll(2))
  ZYHATM = 0.5 * (ZYHAT_ll(1)+ZYHAT_ll(2))
  CALL SM_LATLON(XLATORI,XLONORI,ZXHATM,ZYHATM,ZLATOR,ZLONOR)
  DEALLOCATE(ZXHAT_ll,ZYHAT_ll)
!
  CALL IO_WRITE_FIELD(TPFILE,'LONOR',ZLONOR)
  CALL IO_WRITE_FIELD(TPFILE,'LATOR',ZLATOR)
!
END IF 
!
CALL IO_WRITE_FIELD(TPFILE,'THINSHELL',LTHINSHELL)
CALL IO_WRITE_FIELD(TPFILE,'LAT0',XLAT0)
CALL IO_WRITE_FIELD(TPFILE,'LON0',XLON0)
CALL IO_WRITE_FIELD(TPFILE,'BETA',XBETA)
!
CALL IO_WRITE_FIELD(TPFILE,'XHAT',XXHAT)
CALL IO_WRITE_FIELD(TPFILE,'YHAT',XYHAT)
CALL IO_WRITE_FIELD(TPFILE,'ZHAT',XZHAT)
CALL IO_WRITE_FIELD(TPFILE,'ZTOP',XZTOP)
!
CALL IO_WRITE_FIELD(TPFILE,'ZS',   XZS)
CALL IO_WRITE_FIELD(TPFILE,'ZSMT', XZSMT)
CALL IO_WRITE_FIELD(TPFILE,'SLEVE',LSLEVE)
!
IF (LSLEVE) THEN
  CALL IO_WRITE_FIELD(TPFILE,'LEN1',XLEN1)
  CALL IO_WRITE_FIELD(TPFILE,'LEN2',XLEN2)
END IF
!
!
CALL IO_WRITE_FIELD(TPFILE,'DTMOD',TDTMOD)
CALL IO_WRITE_FIELD(TPFILE,'DTCUR',TDTCUR)
CALL IO_WRITE_FIELD(TPFILE,'DTEXP',TDTEXP)
CALL IO_WRITE_FIELD(TPFILE,'DTSEG',TDTSEG)
!
!*       1.3    Configuration  variables :
!
CALL IO_WRITE_FIELD(TPFILE,'CARTESIAN',LCARTESIAN)
CALL IO_WRITE_FIELD(TPFILE,'LBOUSS',   LBOUSS)
!
IF (LCARTESIAN .AND. LWIND_ZM) THEN
  LWIND_ZM=.FALSE.
  PRINT*,'YOU ARE IN CARTESIAN GEOMETRY SO LWIND_ZM IS FORCED TO FALSE'
END IF
!*       1.4    Reference state variables :
!
CALL IO_WRITE_FIELD(TPFILE,'RHOREFZ',XRHODREFZ)
CALL IO_WRITE_FIELD(TPFILE,'THVREFZ',XTHVREFZ)
CALL IO_WRITE_FIELD(TPFILE,'EXNTOP', XEXNTOP)
!
CALL IO_WRITE_FIELD(TPFILE,'RHODREF',XRHODREF)
CALL IO_WRITE_FIELD(TPFILE,'THVREF', XTHVREF)
!
!
!*       1.5    Variables necessary for plots
!
! PABST,THT,POVOM for cross sections at constant pressure 
! level or constant theta level or constant PV level
!
IF (INDEX(CISO,'PR') /= 0) THEN
  CALL IO_WRITE_FIELD(TPFILE,'PABST',XPABST)
END IF
!
IF (INDEX(CISO,'TK') /= 0) THEN
  CALL IO_WRITE_FIELD(TPFILE,'THT',XTHT)
END IF
!
ZCORIOZ(:,:,:)=SPREAD( XCORIOZ(:,:),DIM=3,NCOPIES=IKU )
ZVOX(:,:,:)=GY_W_VW(1,IKU,1,XWT,XDYY,XDZZ,XDZY)-GZ_V_VW(1,IKU,1,XVT,XDZZ)
ZVOX(:,:,2)=ZVOX(:,:,3)
ZVOY(:,:,:)=GZ_U_UW(1,IKU,1,XUT,XDZZ)-GX_W_UW(1,IKU,1,XWT,XDXX,XDZZ,XDZX)
ZVOY(:,:,2)=ZVOY(:,:,3)
ZVOZ(:,:,:)=GX_V_UV(1,IKU,1,XVT,XDXX,XDZZ,XDZX)-GY_U_UV(1,IKU,1,XUT,XDYY,XDZZ,XDZY)
ZVOZ(:,:,2)=ZVOZ(:,:,3)
ZVOZ(:,:,1)=ZVOZ(:,:,3)
ZWORK31(:,:,:)=GX_M_M(1,IKU,1,XTHT,XDXX,XDZZ,XDZX)
ZWORK32(:,:,:)=GY_M_M(1,IKU,1,XTHT,XDYY,XDZZ,XDZY)
ZWORK33(:,:,:)=GZ_M_M(1,IKU,1,XTHT,XDZZ)
ZPOVO(:,:,:)= ZWORK31(:,:,:)*MZF(1,IKU,1,MYF(ZVOX(:,:,:)))     &
             + ZWORK32(:,:,:)*MZF(1,IKU,1,MXF(ZVOY(:,:,:)))     &
             + ZWORK33(:,:,:)*(MYF(MXF(ZVOZ(:,:,:))) + ZCORIOZ(:,:,:))
ZPOVO(:,:,:)= ZPOVO(:,:,:)*1E6/XRHODREF(:,:,:)
ZPOVO(:,:,1)  =-1.E+11
ZPOVO(:,:,IKU)=-1.E+11
IF (INDEX(CISO,'EV') /= 0) THEN
  TZFIELD%CMNHNAME   = 'POVOT'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'POVOT'
  TZFIELD%CUNITS     = 'PVU' ! 1 PVU = 1e-6 m^2 s^-1 K kg^-1
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'X_Y_Z_POtential VOrticity'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZPOVO)
END IF
!
!
IF (LVAR_RS) THEN
  CALL IO_WRITE_FIELD(TPFILE,'UT',XUT)
  CALL IO_WRITE_FIELD(TPFILE,'VT',XVT)
  !
  IF (LWIND_ZM) THEN
    TZFIELD2(1)%CMNHNAME   = 'UM_ZM'
    TZFIELD2(1)%CSTDNAME   = ''
    TZFIELD2(1)%CLONGNAME  = 'UM_ZM'
    TZFIELD2(1)%CUNITS     = 'm s-1'
    TZFIELD2(1)%CDIR       = 'XY'
    TZFIELD2(1)%CCOMMENT   = 'Zonal component of horizontal wind'
    TZFIELD2(1)%NGRID      = 2
    TZFIELD2(1)%NTYPE      = TYPEREAL
    TZFIELD2(1)%NDIMS      = 3
    TZFIELD2(1)%LTIMEDEP   = .TRUE.
    !
    TZFIELD2(2)%CMNHNAME   = 'VM_ZM'
    TZFIELD2(2)%CSTDNAME   = ''
    TZFIELD2(2)%CLONGNAME  = 'VM_ZM'
    TZFIELD2(2)%CUNITS     = 'm s-1'
    TZFIELD2(2)%CDIR       = 'XY'
    TZFIELD2(2)%CCOMMENT   = 'Meridian component of horizontal wind'
    TZFIELD2(2)%NGRID      = 3
    TZFIELD2(2)%NTYPE      = TYPEREAL
    TZFIELD2(2)%NDIMS      = 3
    TZFIELD2(2)%LTIMEDEP   = .TRUE.
    !
    CALL UV_TO_ZONAL_AND_MERID(XUT,XVT,23,TPFILE=TPFILE,TZFIELDS=TZFIELD2)
  END IF
  !
  CALL IO_WRITE_FIELD(TPFILE,'WT',XWT)
  !
  !   write mixing ratio for water vapor required to plot radio-soundings
  !
  IF (LUSERV) THEN
    CALL IO_WRITE_FIELD(TPFILE,'RVT',XRT(:,:,:,IDX_RVT))
  END IF
END IF
!
!*   Latitude and Longitude arrays
!
IF (.NOT.LCARTESIAN) THEN
  CALL IO_WRITE_FIELD(TPFILE,'LAT',XLAT)
  CALL IO_WRITE_FIELD(TPFILE,'LON',XLON)
END IF
!
!
!-------------------------------------------------------------------------------
!
!*       1.6    Other pronostic variables
!
ZTEMP(:,:,:)=XTHT(:,:,:)*(XPABST(:,:,:)/ XP00) **(XRD/XCPD)
!
IF (LVAR_TURB) THEN
  IF (CTURB /= 'NONE') THEN
    CALL IO_WRITE_FIELD(TPFILE,'TKET',XTKET)
    !
    IF( NRR > 1 ) THEN
      CALL IO_WRITE_FIELD(TPFILE,'SRCT',XSRCT)
      CALL IO_WRITE_FIELD(TPFILE,'SIGS',XSIGS)
    END IF
    ! 
    IF(CTOM=='TM06') THEN
      CALL IO_WRITE_FIELD(TPFILE,'BL_DEPTH',XBL_DEPTH)
    END IF
  END IF
END IF
!
!* Rains
!
IF (LVAR_PR .AND. LUSERR .AND. SIZE(XINPRR)>0 ) THEN
  !
  ! explicit species
  !
  CALL FIND_FIELD_ID_FROM_MNHNAME('INPRR',IID,IRESP)
  TZFIELD = TFIELDLIST(IID)
  TZFIELD%CUNITS = 'mm hour-1'
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XINPRR*3.6E6)
  !
  CALL IO_WRITE_FIELD(TPFILE,'INPRR3D',XINPRR3D)
  CALL IO_WRITE_FIELD(TPFILE,'EVAP3D', XEVAP3D)
  !
  CALL FIND_FIELD_ID_FROM_MNHNAME('ACPRR',IID,IRESP)
  TZFIELD = TFIELDLIST(IID)
  TZFIELD%CUNITS = 'mm'
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XACPRR*1.0E3)
  !
  IF (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'C2R2' .OR. CCLOUD == 'C3R5' .OR.&
      CCLOUD == 'KHKO' .OR. CCLOUD == 'LIMA') THEN 
    IF (SIZE(XINPRC) /= 0 ) THEN
      CALL FIND_FIELD_ID_FROM_MNHNAME('INPRC',IID,IRESP)
      TZFIELD = TFIELDLIST(IID)
      TZFIELD%CUNITS = 'mm hour-1'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XINPRC*3.6E6)
      !
      CALL FIND_FIELD_ID_FROM_MNHNAME('ACPRC',IID,IRESP)
      TZFIELD = TFIELDLIST(IID)
      TZFIELD%CUNITS = 'mm'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XACPRC*1.0E3)
    END IF 
    IF (SIZE(XINDEP) /= 0 ) THEN
      CALL FIND_FIELD_ID_FROM_MNHNAME('INDEP',IID,IRESP)
      TZFIELD = TFIELDLIST(IID)
      TZFIELD%CUNITS = 'mm hour-1'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XINDEP*3.6E6)
      !
      CALL FIND_FIELD_ID_FROM_MNHNAME('ACDEP',IID,IRESP)
      TZFIELD = TFIELDLIST(IID)
      TZFIELD%CUNITS = 'mm'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XACDEP*1.0E3)
    END IF 
  END IF 
  IF (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'C3R5' .OR. CCLOUD == 'LIMA') THEN
    CALL FIND_FIELD_ID_FROM_MNHNAME('INPRS',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm hour-1'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XINPRS*3.6E6)
    !
    CALL FIND_FIELD_ID_FROM_MNHNAME('ACPRS',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XACPRS*1.0E3)
    !
    CALL FIND_FIELD_ID_FROM_MNHNAME('INPRG',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm hour-1'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XINPRG*3.6E6)
    !
    CALL FIND_FIELD_ID_FROM_MNHNAME('ACPRG',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XACPRG*1.0E3)
  !
    IF (SIZE(XINPRH) /= 0 ) THEN
      CALL FIND_FIELD_ID_FROM_MNHNAME('INPRH',IID,IRESP)
      TZFIELD = TFIELDLIST(IID)
      TZFIELD%CUNITS = 'mm hour-1'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XINPRH*3.6E6)
      !
      CALL FIND_FIELD_ID_FROM_MNHNAME('ACPRH',IID,IRESP)
      TZFIELD = TFIELDLIST(IID)
      TZFIELD%CUNITS = 'mm'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XACPRH*1.0E3)
    ENDIF
  !
    ZWORK21(:,:) = XINPRR(:,:) + XINPRS(:,:) + XINPRG(:,:)
    IF (SIZE(XINPRC) /= 0 ) &     
      ZWORK21(:,:) = ZWORK21(:,:) + XINPRC(:,:)
    IF (SIZE(XINPRH) /= 0 ) &       
      ZWORK21(:,:) = ZWORK21(:,:) + XINPRH(:,:)
    CALL FIND_FIELD_ID_FROM_MNHNAME('INPRT',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm hour-1'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK21*3.6E6)
  !
    ZWORK21(:,:) = (XACPRR(:,:) + XACPRS(:,:) + XACPRG(:,:))*1.0E3
    IF (SIZE(XINPRC) /= 0 ) &      
      ZWORK21(:,:) = ZWORK21(:,:) + XACPRC(:,:)*1.0E3
    IF (SIZE(XINPRH) /= 0 ) &        
      ZWORK21(:,:) = ZWORK21(:,:) + XACPRH(:,:)*1.0E3
  !
    CALL FIND_FIELD_ID_FROM_MNHNAME('ACPRT',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK21*1.0E3)
  !
  END IF
  !
  !* Convective rain
  !
  IF (CDCONV /= 'NONE') THEN
    CALL FIND_FIELD_ID_FROM_MNHNAME('PRCONV',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm hour-1'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XPRCONV*3.6E6)
    !
    CALL FIND_FIELD_ID_FROM_MNHNAME('PACCONV',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XPACCONV*1.0E3)
    !
    CALL FIND_FIELD_ID_FROM_MNHNAME('PRSCONV',IID,IRESP)
    TZFIELD = TFIELDLIST(IID)
    TZFIELD%CUNITS = 'mm hour-1'
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XPRSCONV*3.6E6)
  END IF
END IF
IF (LVAR_PR ) THEN
  !Precipitable water in kg/m**2 
  ZWORK21(:,:) = 0.
  ZWORK22(:,:) = 0.
  ZWORK23(:,:) = 0.
  ZWORK31(:,:,:) = DZF(1,IKU,1,XZZ(:,:,:))
  DO JK = IKB,IKE
    !* Calcul de qtot
    IF  (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'LIMA') THEN
      ZWORK23(IIB:IIE,IJB:IJE) = XRT(IIB:IIE,IJB:IJE,JK,1) + &
      XRT(IIB:IIE,IJB:IJE,JK,2) + XRT(IIB:IIE,IJB:IJE,JK,3) + &
      XRT(IIB:IIE,IJB:IJE,JK,4) + XRT(IIB:IIE,IJB:IJE,JK,5) + &
      XRT(IIB:IIE,IJB:IJE,JK,6)
    ELSE
      ZWORK23(IIB:IIE,IJB:IJE) = XRT(IIB:IIE,IJB:IJE,JK,1)
    ENDIF
    !* Calcul de l'eau precipitable
    ZWORK21(IIB:IIE,IJB:IJE)=XRHODREF(IIB:IIE,IJB:IJE,JK)* &
    ZWORK23(IIB:IIE,IJB:IJE)* ZWORK31(IIB:IIE,IJB:IJE,JK)
    !* Sum 
    ZWORK22(IIB:IIE,IJB:IJE) = ZWORK22(IIB:IIE,IJB:IJE)+ZWORK21(IIB:IIE,IJB:IJE)
    ZWORK21(:,:) = 0.
    ZWORK23(:,:) = 0.
  END DO
  !* Precipitable water in kg/m**2
  TZFIELD%CMNHNAME   = 'PRECIP_WAT'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'PRECIP_WAT'
  TZFIELD%CUNITS     = 'kg m-2'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = ''
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 2
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK22)
ENDIF
!
!
!* Flux d'humidit et d'hydromtores
IF (LHU_FLX) THEN
  ZWORK35(:,:,:) = XRHODREF(:,:,:) * XRT(:,:,:,1)
  ZWORK31(:,:,:) = MXM(ZWORK35(:,:,:)) * XUT(:,:,:)
  ZWORK32(:,:,:) = MYM(ZWORK35(:,:,:)) * XVT(:,:,:)
  IF  (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'LIMA') THEN
    ZWORK36(:,:,:) = ZWORK35(:,:,:) + XRHODREF(:,:,:) * (XRT(:,:,:,2) + &
    XRT(:,:,:,3) + XRT(:,:,:,4) + XRT(:,:,:,5) + XRT(:,:,:,6))
    ZWORK33(:,:,:) = MXM(ZWORK36(:,:,:)) * XUT(:,:,:)
    ZWORK34(:,:,:) = MYM(ZWORK36(:,:,:)) * XVT(:,:,:)
  ENDIF
  ZWORK35(:,:,:) = GX_U_M(1,IKU,1,ZWORK31,XDXX,XDZZ,XDZX) + GY_V_M(1,IKU,1,ZWORK32,XDYY,XDZZ,XDZY)
  ZWORK36(:,:,:) = GX_U_M(1,IKU,1,ZWORK33,XDXX,XDZZ,XDZX) + GY_V_M(1,IKU,1,ZWORK34,XDYY,XDZZ,XDZY)
  !
  ! Integration sur 3000 m
  !
  IKTOP(:,:)=0
  DO JK=1,IKU-1
    WHERE (((XZZ(:,:,JK) -XZS(:,:))<= 3000.0) .AND. ((XZZ(:,:,JK+1) -XZS(:,:))> 3000.0))
      IKTOP(:,:)=JK
    END WHERE
  END DO
  ZDELTAZ(:,:,:)=DZF(1,IKU,1,XZZ) 
  ZWORK21(:,:) = 0.
  ZWORK22(:,:) = 0.
  ZWORK25(:,:) = 0.  
  DO JJ=1,IJU
    DO JI=1,IIU
      IAUX=IKTOP(JI,JJ)
      DO JK=IKB,IAUX-1 
        ZWORK21(JI,JJ) = ZWORK21(JI,JJ) + ZWORK31(JI,JJ,JK) * ZDELTAZ(JI,JJ,JK)
        ZWORK22(JI,JJ) = ZWORK22(JI,JJ) + ZWORK32(JI,JJ,JK) * ZDELTAZ(JI,JJ,JK)
        ZWORK25(JI,JJ) = ZWORK25(JI,JJ) + ZWORK35(JI,JJ,JK) * ZDELTAZ(JI,JJ,JK)
      ENDDO
      IF (IAUX >= IKB) THEN
        ZDELTAZ(JI,JJ,IAUX)= 3000. - (XZZ(JI,JJ,IAUX) -XZS(JI,JJ))
        ZWORK21(JI,JJ) = ZWORK21(JI,JJ) + ZWORK31(JI,JJ,IAUX) * ZDELTAZ(JI,JJ,IAUX) 
        ZWORK22(JI,JJ) = ZWORK22(JI,JJ) + ZWORK32(JI,JJ,IAUX) * ZDELTAZ(JI,JJ,IAUX)
        ZWORK25(JI,JJ) = ZWORK25(JI,JJ) + ZWORK35(JI,JJ,IAUX) * ZDELTAZ(JI,JJ,IAUX)
      ENDIF
    ENDDO
  ENDDO
  IF  (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'LIMA') THEN
    ZWORK23(:,:) = 0.
    ZWORK24(:,:) = 0.
    ZWORK26(:,:) = 0.
    DO JJ=1,IJU
      DO JI=1,IIU
        IAUX=IKTOP(JI,JJ)
        DO JK=IKB,IAUX-1 
          ZWORK23(JI,JJ) = ZWORK23(JI,JJ) + ZWORK33(JI,JJ,JK) * ZDELTAZ(JI,JJ,JK)
          ZWORK24(JI,JJ) = ZWORK24(JI,JJ) + ZWORK34(JI,JJ,JK) * ZDELTAZ(JI,JJ,JK)
          ZWORK26(JI,JJ) = ZWORK26(JI,JJ) + ZWORK36(JI,JJ,JK) * ZDELTAZ(JI,JJ,JK)
        ENDDO
        IF (IAUX >= IKB) THEN
          ZDELTAZ(JI,JJ,IAUX)= 3000. - (XZZ(JI,JJ,IAUX) -XZS(JI,JJ))
          ZWORK23(JI,JJ) = ZWORK23(JI,JJ) + ZWORK33(JI,JJ,IAUX) * ZDELTAZ(JI,JJ,IAUX) 
          ZWORK24(JI,JJ) = ZWORK24(JI,JJ) + ZWORK34(JI,JJ,IAUX) * ZDELTAZ(JI,JJ,IAUX)
          ZWORK26(JI,JJ) = ZWORK26(JI,JJ) + ZWORK36(JI,JJ,IAUX) * ZDELTAZ(JI,JJ,IAUX)
        ENDIF
      ENDDO
    ENDDO
  ENDIF
  ! Ecriture
  !  composantes U et V du flux surfacique d'humidit  
  TZFIELD%CMNHNAME   = 'UM90'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'UM90'
  TZFIELD%CUNITS     = 'kg s-1 m-2'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = ''
  TZFIELD%NGRID      = 2
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK31)
  !  
  TZFIELD%CMNHNAME   = 'VM90'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'VM90'
  TZFIELD%CUNITS     = 'kg s-1 m-2'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = ''
  TZFIELD%NGRID      = 3
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK32)
  !  composantes U et V du flux d'humidit intgr sur 3000 metres
  TZFIELD%CMNHNAME   = 'UM91'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'UM91'
  TZFIELD%CUNITS     = 'kg s-1 m-1'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = ''
  TZFIELD%NGRID      = 2
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 2
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK21)
  !
  TZFIELD%CMNHNAME   = 'VM91'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'VM91'
  TZFIELD%CUNITS     = 'kg s-1 m-1'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = ''
  TZFIELD%NGRID      = 3
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 2
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK22)
  !
  !   Convergence d'humidit
  TZFIELD%CMNHNAME   = 'HMCONV'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'HMCONV'
  TZFIELD%CUNITS     = 'kg s-1 m-3'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'X_Y_Horizontal CONVergence of moisture flux'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 3
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,-ZWORK35)
  !
  !   Convergence d'humidit intgr sur 3000 mtres
  TZFIELD%CMNHNAME   = 'HMCONV3000'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'HMCONV3000'
  TZFIELD%CUNITS     = 'kg s-1 m-3'
  TZFIELD%CDIR       = 'XY'
  TZFIELD%CCOMMENT   = 'X_Y_Horizontal CONVergence of moisture flux'
  TZFIELD%NGRID      = 1
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 2
  TZFIELD%LTIMEDEP   = .TRUE.
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,-ZWORK25)
  !
  IF  (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'LIMA') THEN
    !  composantes U et V du flux surfacique d'hydromtores  
    TZFIELD%CMNHNAME   = 'UM92'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'UM92'
    TZFIELD%CUNITS     = 'kg s-1 m-2'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = ''
    TZFIELD%NGRID      = 2
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 3
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK33)
    ! 
    TZFIELD%CMNHNAME   = 'VM92'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'VM92'
    TZFIELD%CUNITS     = 'kg s-1 m-2'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = ''
    TZFIELD%NGRID      = 3
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 3
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK34)
    !  composantes U et V du flux d'hydromtores intgr sur 3000 metres
    TZFIELD%CMNHNAME   = 'UM93'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'UM93'
    TZFIELD%CUNITS     = 'kg s-1 m-1'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = ''
    TZFIELD%NGRID      = 2
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 2
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK23)
    ! 
    TZFIELD%CMNHNAME   = 'VM93'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'VM93'
    TZFIELD%CUNITS     = 'kg s-1 m-1'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = ''
    TZFIELD%NGRID      = 3
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 2
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,ZWORK24)
    !   Convergence d'hydromtores
    TZFIELD%CMNHNAME   = 'HMCONV_TT'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'HMCONV_TT'
    TZFIELD%CUNITS     = 'kg s-1 m-3'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = 'X_Y_Horizontal CONVergence of hydrometeor flux'
    TZFIELD%NGRID      = 1
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 3
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,-ZWORK36)
    !   Convergence d'hydromtores intgr sur 3000 mtres
    TZFIELD%CMNHNAME   = 'HMCONV3000_TT'
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CLONGNAME  = 'HMCONV3000_TT'
    TZFIELD%CUNITS     = 'kg s-1 m-3'
    TZFIELD%CDIR       = 'XY'
    TZFIELD%CCOMMENT   = 'X_Y_Horizontal CONVergence of hydrometeor flux'
    TZFIELD%NGRID      = 1
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 2
    TZFIELD%LTIMEDEP   = .TRUE.
    CALL IO_WRITE_FIELD(TPFILE,TZFIELD,-ZWORK26)
  ENDIF
ENDIF
!
!* Moist variables
!
IF (LVAR_MRW .OR. LLIMA_DIAG) THEN
  IF (NRR >=1) THEN
    ! Moist variables are written individually in file
    TZFIELD%CSTDNAME   = ''
    TZFIELD%CDIR       = 'XY'
    TZFIELD%NGRID      = 1
    TZFIELD%NTYPE      = TYPEREAL
    TZFIELD%NDIMS      = 3
    TZFIELD%LTIMEDEP   = .TRUE.
    IF (LUSERV) THEN
      TZFIELD%CMNHNAME   = 'MRV'
      TZFIELD%CLONGNAME  = 'MRV'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRV'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XRT(:,:,:,IDX_RVT)*1.E3)
    END IF
    IF (LUSERC) THEN
      TZFIELD%CMNHNAME   = 'MRC'
      TZFIELD%CLONGNAME  = 'MRC'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRC'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XRT(:,:,:,IDX_RCT)*1.E3)
!
      TZFIELD%CMNHNAME   = 'VRC'
      TZFIELD%CLONGNAME  = 'VRC'
      TZFIELD%CUNITS     = '' !vol/vol
      TZFIELD%CCOMMENT   = 'X_Y_Z_VRC (vol/vol)'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XRT(:,:,:,IDX_RCT)*XRHODREF(:,:,:)/1.E3)
    END IF
    IF (LUSERR) THEN
      TZFIELD%CMNHNAME   = 'MRR'
      TZFIELD%CLONGNAME  = 'MRR'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRR'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XRT(:,:,:,IDX_RRT)*1.E3)
!
      TZFIELD%CMNHNAME   = 'VRR'
      TZFIELD%CLONGNAME  = 'VRR'
      TZFIELD%CUNITS     = '' !vol/vol
      TZFIELD%CCOMMENT   = 'X_Y_Z_VRR (vol/vol)'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XRT(:,:,:,IDX_RRT)*XRHODREF(:,:,:)/1.E3)
    END IF
    IF (LUSERI) THEN
      TZFIELD%CMNHNAME   = 'MRI'
      TZFIELD%CLONGNAME  = 'MRI'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRI'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,XRT(:,:,:,IDX_RIT)*1.E3)
!
      IF (LUSECI) THEN