write_lfifm1_for_diag.f90

!MNH_LIC Copyright 1994-2014 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.
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source: /home/cvsroot/MNH-VX-Y-Z/src/MNH/write_lfifm1_for_diag.f90,v $ $Revision: 1.3.2.5.2.4.2.3.2.3.2.9.2.2 $
! masdev4_7 BUG1 2007/06/15 17:47:18
!-----------------------------------------------------------------
!################################
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
!-------------------------------------------------------------------------------
!
!*       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 : CLUOUT0
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_CH_AEROSOL
USE MODD_CH_AERO_n
USE MODD_CH_MNHC_n
USE MODE_DUST_PSD
USE MODE_SALT_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_GATHER_ll
!
USE MODE_GRIDPROJ
USE MODE_FIELD
USE MODE_FMWRIT
USE MODE_ll
USE MODE_IO_ll
USE MODE_THERMO
USE MODE_MODELN_HANDLER
USE MODI_LIDAR
!
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 
INTEGER           :: IGRID          ! IGRID : grid indicator
INTEGER           :: ILENCH         ! ILENCH : length of comment string 
!
CHARACTER(LEN=28) :: YFMFILE        ! Temporary variable to store FM-file name
CHARACTER(LEN=28) :: YFMFILE2       ! Temporary variable to store FM-file name
CHARACTER(LEN=16) :: YRECFM         ! Name of the article to be written
CHARACTER(LEN=100):: YCOMMENT       ! Comment string
!
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  :: ZRHOT, ZTMP ! work array

!ECRITURE DANS UN FICHIER ASCII DE RESULTATS 
!INITIALISATION DU NOM DE FICHIER CREE EN PARALLELE AVEC CELUI LFI
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
INTEGER           :: I,IID
TYPE(TFIELDDATA)              :: TZFIELD
TYPE(TFIELDDATA),DIMENSION(2) :: TZFIELD2
!
! LIMA LIDAR
REAL,DIMENSION(:,:,:,:), ALLOCATABLE :: ZTMP1, ZTMP2, ZTMP3, ZTMP4
!
!-------------------------------------------------------------------------------
!
!*       0.     ARRAYS BOUNDS INITIALIZATION
!
YFMFILE2 = TPFILE%CNAME
!
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()
CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT0,IRESP)
!-------------------------------------------------------------------------------
!
!*       1.     WRITES IN THE LFI FILE
!               ---------------------- 
!
!*       1.0    YFMFILE2 and HDADFILE :
!
CALL IO_WRITE_FIELD(TPFILE,'MASDEV',   CLUOUT,IRESP,NMASDEV)
CALL IO_WRITE_FIELD(TPFILE,'BUGFIX',   CLUOUT,IRESP,NBUGFIX)
CALL IO_WRITE_FIELD(TPFILE,'BIBUSER',  CLUOUT,IRESP,CBIBUSER)
CALL IO_WRITE_FIELD(TPFILE,'PROGRAM',  CLUOUT,IRESP,CPROGRAM)
!
CALL IO_WRITE_FIELD(TPFILE,'L1D',      CLUOUT,IRESP,L1D)
CALL IO_WRITE_FIELD(TPFILE,'L2D',      CLUOUT,IRESP,L2D)
CALL IO_WRITE_FIELD(TPFILE,'PACK',     CLUOUT,IRESP,LPACK)
!
CALL IO_WRITE_FIELD(TPFILE,'MY_NAME',  CLUOUT,IRESP,TPFILE%CNAME)
CALL IO_WRITE_FIELD(TPFILE,'DAD_NAME', CLUOUT,IRESP,HDADFILE)
!
IF (LEN_TRIM(HDADFILE)>0) THEN
  CALL IO_WRITE_FIELD(TPFILE,'DXRATIO',CLUOUT,IRESP,NDXRATIO_ALL(1))
  CALL IO_WRITE_FIELD(TPFILE,'DYRATIO',CLUOUT,IRESP,NDYRATIO_ALL(1))
  CALL IO_WRITE_FIELD(TPFILE,'XOR',    CLUOUT,IRESP,NXOR_ALL(1))
  CALL IO_WRITE_FIELD(TPFILE,'YOR',    CLUOUT,IRESP,NYOR_ALL(1))
END IF
!
CALL IO_WRITE_FIELD(TPFILE,'SURF',     CLUOUT,IRESP,CSURF)
!
!*       1.1    Type and Dimensions :
!
CALL IO_WRITE_FIELD(TPFILE,'STORAGE_TYPE',CLUOUT,IRESP,'DI')
!
CALL IO_WRITE_FIELD(TPFILE,'IMAX',CLUOUT,IRESP,NIMAX_ll)
CALL IO_WRITE_FIELD(TPFILE,'JMAX',CLUOUT,IRESP,NJMAX_ll)
CALL IO_WRITE_FIELD(TPFILE,'KMAX',CLUOUT,IRESP,NKMAX)
!
CALL IO_WRITE_FIELD(TPFILE,'JPHEXT',CLUOUT,IRESP,JPHEXT)
!
!*       1.2    Grid variables :
!
IF (.NOT.LCARTESIAN) THEN
  CALL IO_WRITE_FIELD(TPFILE,'RPK',   CLUOUT,IRESP,XRPK)
  CALL IO_WRITE_FIELD(TPFILE,'LONORI',CLUOUT,IRESP,XLONORI)
  CALL IO_WRITE_FIELD(TPFILE,'LATORI',CLUOUT,IRESP,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)
!
  !LONOR and LATOR not in TFIELDLIST because local variables
  TZFIELD%CMNHNAME   = 'LONOR'
  TZFIELD%CSTDNAME   = ''
  TZFIELD%CLONGNAME  = 'MesoNH: LONOR'
  TZFIELD%CUNITS     = 'degree'
  TZFIELD%CDIR       = '--'
  TZFIELD%CCOMMENT   = 'Longitude of 1st mass point'
  TZFIELD%NGRID      = 0
  TZFIELD%NTYPE      = TYPEREAL
  TZFIELD%NDIMS      = 0
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,ZLONOR)
!
  TZFIELD%CMNHNAME   = 'LATOR'
  TZFIELD%CLONGNAME  = 'MesoNH: LATOR'
  TZFIELD%CCOMMENT   = 'Latitude of 1st mass point'
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,ZLATOR)
!
END IF 
!
CALL IO_WRITE_FIELD(TPFILE,'THINSHELL',CLUOUT,IRESP,LTHINSHELL)
CALL IO_WRITE_FIELD(TPFILE,'LAT0',CLUOUT,IRESP,XLAT0)
CALL IO_WRITE_FIELD(TPFILE,'LON0',CLUOUT,IRESP,XLON0)
CALL IO_WRITE_FIELD(TPFILE,'BETA',CLUOUT,IRESP,XBETA)
!
CALL IO_WRITE_FIELD(TPFILE,'XHAT',CLUOUT,IRESP,XXHAT)
CALL IO_WRITE_FIELD(TPFILE,'YHAT',CLUOUT,IRESP,XYHAT)
CALL IO_WRITE_FIELD(TPFILE,'ZHAT',CLUOUT,IRESP,XZHAT)
!
CALL IO_WRITE_FIELD(TPFILE,'ZS',   CLUOUT,IRESP,XZS)
CALL IO_WRITE_FIELD(TPFILE,'ZSMT', CLUOUT,IRESP,XZSMT)
CALL IO_WRITE_FIELD(TPFILE,'SLEVE',CLUOUT,IRESP,LSLEVE)
!
IF (LSLEVE) THEN
  CALL IO_WRITE_FIELD(TPFILE,'LEN1',CLUOUT,IRESP,XLEN1)
  CALL IO_WRITE_FIELD(TPFILE,'LEN2',CLUOUT,IRESP,XLEN2)
END IF
!
!
CALL IO_WRITE_FIELD(TPFILE,'DTMOD',CLUOUT,IRESP,TDTMOD)
CALL IO_WRITE_FIELD(TPFILE,'DTCUR',CLUOUT,IRESP,TDTCUR)
CALL IO_WRITE_FIELD(TPFILE,'DTEXP',CLUOUT,IRESP,TDTEXP)
CALL IO_WRITE_FIELD(TPFILE,'DTSEG',CLUOUT,IRESP,TDTSEG)
!
!*       1.3    Configuration  variables :
!
CALL IO_WRITE_FIELD(TPFILE,'CARTESIAN',CLUOUT,IRESP,LCARTESIAN)
CALL IO_WRITE_FIELD(TPFILE,'LBOUSS',   CLUOUT,IRESP,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',CLUOUT,IRESP,XRHODREFZ)
CALL IO_WRITE_FIELD(TPFILE,'THVREFZ',CLUOUT,IRESP,XTHVREFZ)
CALL IO_WRITE_FIELD(TPFILE,'EXNTOP', CLUOUT,IRESP,XEXNTOP)
!
CALL IO_WRITE_FIELD(TPFILE,'RHODREF',CLUOUT,IRESP,XRHODREF)
CALL IO_WRITE_FIELD(TPFILE,'THVREF', CLUOUT,IRESP,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',CLUOUT,IRESP,XPABST)
END IF
!
IF (INDEX(CISO,'TK') /= 0) THEN
  CALL IO_WRITE_FIELD(TPFILE,'THT',CLUOUT,IRESP,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  = 'MesoNH: 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
  CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,ZPOVO)
END IF
!
!
IF (LVAR_RS) THEN
  CALL IO_WRITE_FIELD(TPFILE,'UT',CLUOUT,IRESP,XUT)
  CALL IO_WRITE_FIELD(TPFILE,'VT',CLUOUT,IRESP,XVT)
  !
  IF (LWIND_ZM) THEN
    TZFIELD2(1)%CMNHNAME   = 'UM_ZM'
    TZFIELD2(1)%CSTDNAME   = ''
    TZFIELD2(1)%CLONGNAME  = 'MesoNH: 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(2)%CMNHNAME   = 'VM_ZM'
    TZFIELD2(2)%CSTDNAME   = ''
    TZFIELD2(2)%CLONGNAME  = 'MesoNH: 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
    !
    CALL UV_TO_ZONAL_AND_MERID(XUT,XVT,23,TPFILE=TPFILE,TZFIELDS=TZFIELD2)
  END IF
  !
  CALL IO_WRITE_FIELD(TPFILE,'WT',CLUOUT,IRESP,XWT)
  !
  !   write mixing ratio for water vapor required to plot radio-soundings
  !
  IF (LUSERV) THEN
    CALL IO_WRITE_FIELD(TPFILE,'RVT',CLUOUT,IRESP,XRT(:,:,:,IDX_RVT))
  END IF
END IF
!
!*   Latitude and Longitude arrays
!
IF (.NOT.LCARTESIAN) THEN
  CALL IO_WRITE_FIELD(TPFILE,'LAT',CLUOUT,IRESP,XLAT)
  CALL IO_WRITE_FIELD(TPFILE,'LON',CLUOUT,IRESP,XLON)
END IF
!
!
!-------------------------------------------------------------------------------
!
!*       1.6    Other pronostic variables
!
ZTEMP(:,:,:)=XTHT(:,:,:)*(XPABST(:,:,:)/ XP00) **(XRD/XCPD)
!
IF (LVAR_TURB) THEN
  IF (CTURB /= 'NONE') THEN
    YRECFM='TKET'
    YCOMMENT='X_Y_Z_Turbulent Kinetic Energy (M**2/S**2)'
    IGRID=1
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',XTKET,IGRID,ILENCH,YCOMMENT,IRESP)
    !
    IF( NRR > 1 ) THEN
      YRECFM='SRCT'
      YCOMMENT='X_Y_Z_normalized 2nd_order moment s_r_c/2Sigma_s2 (KG/KG**2)'
      IGRID=1
      ILENCH=LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',XSRCT,IGRID,ILENCH,YCOMMENT,IRESP)
      !
      YRECFM='SIGS'
      YCOMMENT='X_Y_Z_Sigma_s from turbulence scheme (KG/KG**2)'
      IGRID=1
      ILENCH=LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',XSIGS,IGRID,ILENCH,YCOMMENT,IRESP)
    END IF
    ! 
    IF(CTOM=='TM06') THEN
      YRECFM='BL_DEPTH'
      YCOMMENT='X_Y_BL_DEPTH (M)'
      IGRID=1
      ILENCH=LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',XBL_DEPTH,IGRID,ILENCH,YCOMMENT,IRESP)
    END IF
  END IF
END IF
!
!* Rains
!
IF (LVAR_PR .AND. LUSERR .AND. SIZE(XINPRR)>0 ) THEN
  !
  ! explicit species
  !
  ZWORK21(:,:) = XINPRR(:,:)*3.6E6
  YRECFM      ='INPRR'
  YCOMMENT    ='X_Y_INstantaneous PRecipitation Rate (MM/H)'
  IGRID       =1
  ILENCH      =LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
  ZWORK31(:,:,:)  = XINPRR3D(:,:,:)
  YRECFM      = 'INPRR3D'
  YCOMMENT    = 'X_Y_INstantaneous 3D Rain Precipitation flux (M/S)'
  IGRID       = 1
  ILENCH      = LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK31,IGRID,ILENCH, &
                                             YCOMMENT,IRESP) ! unit conversion
!
  ZWORK31(:,:,:)  = XEVAP3D(:,:,:)
  YRECFM      = 'EVAP3D'
  YCOMMENT    = 'X_Y_INstantaneous 3D Rain Evaporation flux (KG/KG/S)'
  IGRID       = 1
  ILENCH      = LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK31,IGRID,ILENCH, &
                                             YCOMMENT,IRESP) ! unit conversion
!
  ZWORK21(:,:) = XACPRR(:,:)*1.E3
  YRECFM      ='ACPRR'
  YCOMMENT    ='X_Y_ACcumulated PRecipitation Rate (MM)'
  IGRID       =1
  ILENCH      =LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
  IF (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'C2R2' .OR. CCLOUD == 'C3R5' .OR.&
      CCLOUD == 'KHKO' .OR. CCLOUD == 'LIMA') THEN 
    IF (SIZE(XINPRC) /= 0 ) THEN
      ZWORK21(:,:) = XINPRC(:,:)*3.6E6
      YRECFM      ='INPRC'
      YCOMMENT    ='X_Y_INstantaneous Cloud PRecipitation Rate (MM/H)'
      IGRID       =1
      ILENCH      =LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
      ZWORK21(:,:) = XACPRC(:,:)*1.E3
      YRECFM      ='ACPRC'
      YCOMMENT    ='X_Y_ACcumulated Cloud PRecipitation Rate (MM)'
      IGRID       =1
      ILENCH      =LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
    END IF 
    IF (SIZE(XINDEP) /= 0 ) THEN
      ZWORK21(:,:) = XINDEP(:,:)*3.6E6
      YRECFM      ='INDEP'
      YCOMMENT    ='X_Y_INstantaneous Cloud Deposition Rate (MM/H)'
      IGRID       =1
      ILENCH      =LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
      ZWORK21(:,:) = XACDEP(:,:)*1.E3
      YRECFM      ='ACDEP'
      YCOMMENT    ='X_Y_ACcumulated Cloud Deposition Rate (MM)'
      IGRID       =1
      ILENCH      =LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
    END IF 
  END IF 
  IF (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'C3R5' .OR. CCLOUD == 'LIMA') THEN
    ZWORK21(:,:)  = XINPRS(:,:)*3.6E6
    YRECFM      = 'INPRS'
    YCOMMENT    = 'X_Y_INstantaneaous PRecipitation rate for Snow (MM/H)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    ZWORK21(:,:)  = XACPRS(:,:)*1.0E3
    YRECFM      = 'ACPRS'
    YCOMMENT    = 'X_Y_ACccumuated PRecipitation rate for Snow (MM)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    ZWORK21(:,:)  = XINPRG(:,:)*3.6E6
    YRECFM      = 'INPRG'
    YCOMMENT    = 'X_Y_INstantaneaous PRecipitation rate for Graupel (MM/H)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    ZWORK21(:,:)  = XACPRG(:,:)*1.0E3
    YRECFM      = 'ACPRG'
    YCOMMENT    = 'X_Y_ACccumuated PRecipitation rate for Graupel (MM)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    IF (SIZE(XINPRH) /= 0 ) THEN
      ZWORK21(:,:)  = XINPRH(:,:)*3.6E6
      YRECFM      = 'INPRH'
      YCOMMENT    = 'X_Y_INstantaneaous PRecipitation rate for Hail (MM/H)'
      IGRID       = 1
      ILENCH      = LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
    !
      ZWORK21(:,:)  = XACPRH(:,:)*1.0E3
      YRECFM      = 'ACPRH'
      YCOMMENT    = 'X_Y_ACccumuated PRecipitation rate for Hail (MM)'
      IGRID       = 1
      ILENCH      = LEN(YCOMMENT)
      CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
    ENDIF
  !
    ZWORK21(:,:) = (XINPRR(:,:) + XINPRS(:,:) + XINPRG(:,:))*3.6E6
    IF (SIZE(XINPRC) /= 0 ) &     
      ZWORK21(:,:) = ZWORK21(:,:) + XINPRC(:,:)*3.6E6          
    IF (SIZE(XINPRH) /= 0 ) &       
      ZWORK21(:,:) = ZWORK21(:,:) + XINPRH(:,:)*3.6E6
    YRECFM      = 'INPRT'
    YCOMMENT    = 'X_Y_Total INstantaneaous PRecipitation rate (MM/H)'
    YCOMMENT    = 'X_Y_INPRT (MM/H)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    ZWORK21(:,:) = (XACPRR(:,:) + XACPRS(:,:) + XACPRG(:,:))*1.0E3
    IF (SIZE(XINPRC) /= 0 ) &      
      ZWORK21(:,:) = ZWORK21(:,:) + XACPRC(:,:)*1.0E3
    IF (SIZE(XINPRH) /= 0 ) &        
      ZWORK21(:,:) = ZWORK21(:,:) + XACPRH(:,:)*1.0E3
  !
    YRECFM      = 'ACPRT'
    YCOMMENT    = 'X_Y_Total ACcumulated PRecipitation rate (MM)'
    YCOMMENT    = 'X_Y_ACPRT (MM)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
  END IF
  !
  !* Convective rain
  !
  IF (CDCONV /= 'NONE') THEN
    ZWORK21(:,:) = XPRCONV(:,:)*3.6E6 
    YRECFM      ='PRCONV'
    YCOMMENT    ='X_Y_CONVective instantaneous Precipitation Rate (MM/H)'
    IGRID       =1
    ILENCH      =LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    ZWORK21(:,:) = XPACCONV(:,:)*1.E3 
    YRECFM      ='PACCONV'
    YCOMMENT    ='X_Y_CONVective ACcumulated Precipitation rate (MM)'
    IGRID       =1
    ILENCH      =LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
    ZWORK21(:,:) = XPRSCONV(:,:)*3.6E6 
    YRECFM      ='PRSCONV'
    YCOMMENT    ='X_Y_CONVective instantaneous Precipitation Rate for Snow (MM/H)'
    IGRID       = 1
    ILENCH      = LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  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
  YRECFM='PRECIP_WAT'
  YCOMMENT='(kg/m)'
  IGRID=1
  ILENCH=LEN(YCOMMENT)  
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK22,IGRID,ILENCH,YCOMMENT,IRESP)
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  
  YRECFM='UM90'
  YCOMMENT='(kg / s / m)'
  IGRID=2
  ILENCH=LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK31,IGRID,ILENCH,YCOMMENT,IRESP)
  !  
  YRECFM='VM90'
  YCOMMENT='(kg / s / m)'
  IGRID=3
  ILENCH=LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK32,IGRID,ILENCH,YCOMMENT,IRESP)
  !  composantes U et V du flux d'humidit intgr sur 3000 metres
  YRECFM='UM91'
  YCOMMENT='(kg / s / m)'
  IGRID=2
  ILENCH=LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK21,IGRID,ILENCH,YCOMMENT,IRESP)
  !
  YRECFM='VM91'
  YCOMMENT='(kg / s / m)'
  IGRID=3
  ILENCH=LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK22,IGRID,ILENCH,YCOMMENT,IRESP)
  !
  !   Convergence d'humidit
  YRECFM='HMCONV'
  YCOMMENT='X_Y_Horizontal CONVergence of moisture flux (kg / s / m)'
  YCOMMENT='(kg / s / m^3)'
  IGRID=1
  ILENCH=LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK35*(-1),IGRID,ILENCH,YCOMMENT,IRESP)
  !
  !   Convergence d'humidit intgr sur 3000 mtres
  YRECFM='HMCONV3000'
  YCOMMENT='X_Y_Horizontal CONVergence of moisture flux (kg / s / m)'
  YCOMMENT='(kg / s / m^3)'
  IGRID=1
  ILENCH=LEN(YCOMMENT)
  CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK25*(-1),IGRID,ILENCH,YCOMMENT,IRESP)
  !
  IF  (CCLOUD(1:3) == 'ICE' .OR. CCLOUD == 'LIMA') THEN
    !  composantes U et V du flux surfacique d'hydromtores  
    YRECFM='UM92'
    YCOMMENT='(kg / s / m)'
    IGRID=2
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK33,IGRID,ILENCH,YCOMMENT,IRESP)
    ! 
    YRECFM='VM92'
    YCOMMENT='(kg / s / m)'
    IGRID=3
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK34,IGRID,ILENCH,YCOMMENT,IRESP)
    !  composantes U et V du flux d'hydromtores intgr sur 3000 metres
    YRECFM='UM93'
    YCOMMENT='(kg / s / m)'
    IGRID=2
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK23,IGRID,ILENCH,YCOMMENT,IRESP)
    YRECFM='VM93'
    YCOMMENT='(kg / s / m)'
    IGRID=3
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK24,IGRID,ILENCH,YCOMMENT,IRESP)
    !   Convergence d'hydromtores
    YRECFM='HMCONV_TT'
    YCOMMENT='X_Y_Horizontal CONVergence of hydrometeor flux (kg / s / m)'
    YCOMMENT='(kg / s / m^3)'
    IGRID=1
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK36*(-1),IGRID,ILENCH,YCOMMENT,IRESP)
    !   Convergence d'hydromtores intgr sur 3000 mtres
    YRECFM='HMCONV3000_TT'
    YCOMMENT='X_Y_Horizontal CONVergence of hydrometeor flux (kg / s / m)'
    YCOMMENT='(kg / s / m^3)'
    IGRID=1
    ILENCH=LEN(YCOMMENT)
    CALL FMWRIT(YFMFILE2,YRECFM,CLUOUT,'XY',ZWORK26*(-1),IGRID,ILENCH,YCOMMENT,IRESP)
  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
    IF (LUSERV) THEN
      TZFIELD%CMNHNAME   = 'MRV'
      TZFIELD%CLONGNAME  = 'MesoNH: MRV'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRV'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RVT)*1.E3)
    END IF
    IF (LUSERC) THEN
      TZFIELD%CMNHNAME   = 'MRC'
      TZFIELD%CLONGNAME  = 'MesoNH: MRC'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRC'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RCT)*1.E3)
!
      TZFIELD%CMNHNAME   = 'VRC'
      TZFIELD%CLONGNAME  = 'MesoNH: VRC'
      TZFIELD%CUNITS     = '' !vol/vol
      TZFIELD%CCOMMENT   = 'X_Y_Z_VRC (vol/vol)'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RCT)*XRHODREF(:,:,:)/1.E3)
    END IF
    IF (LUSERR) THEN
      TZFIELD%CMNHNAME   = 'MRR'
      TZFIELD%CLONGNAME  = 'MesoNH: MRR'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRR'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RRT)*1.E3)
!
      TZFIELD%CMNHNAME   = 'VRR'
      TZFIELD%CLONGNAME  = 'MesoNH: VRR'
      TZFIELD%CUNITS     = '' !vol/vol
      TZFIELD%CCOMMENT   = 'X_Y_Z_VRR (vol/vol)'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RRT)*XRHODREF(:,:,:)/1.E3)
    END IF
    IF (LUSERI) THEN
      TZFIELD%CMNHNAME   = 'MRI'
      TZFIELD%CLONGNAME  = 'MesoNH: MRI'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRI'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RIT)*1.E3)
!
      IF (LUSECI) THEN
        CALL IO_WRITE_FIELD(TPFILE,'CIT',CLUOUT,IRESP,XCIT(:,:,:))
      END IF
    END IF
    IF (LUSERS) THEN
      TZFIELD%CMNHNAME   = 'MRS'
      TZFIELD%CLONGNAME  = 'MesoNH: MRS'
      TZFIELD%CUNITS     = 'g kg-1'
      TZFIELD%CCOMMENT   = 'X_Y_Z_MRS'
      CALL IO_WRITE_FIELD(TPFILE,TZFIELD,CLUOUT,IRESP,XRT(:,:,:,IDX_RST)*1.E3)
    END IF
    IF (LUSERG) THEN
      TZFIELD%CMNHNAME   = 'MRG'
      TZFIELD%CLONGNAME  = 'MesoNH: MRG'
      TZFIELD%CUNITS     = 'g kg-1'