write_diag_misc_isban.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
      SUBROUTINE write_diag_misc_isba_n (DTCO, HSELECT, OSNOWDIMNC, U, OPATCH_BUDGET, D,  &
                                         ND, DM, NDM, IO, S, K, NP, TPSNOW, HPROGRAM)
!     #################################
!
!!****  *WRITE_DIAG_MISC_ISBA* - writes the ISBA diagnostic fields
!!
!!    PURPOSE
!!    -------
!!
!!
!!**  METHOD
!!    ------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      P. Le Moigne   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    10/2004
!!      B. Decharme    2008  Total Albedo, Total SWI and Floodplains
!!      B. Decharme 06/2009  key to write (or not) patch result
!!      A.L. Gibelin 04/09 : Add respiration diagnostics
!!      A.L. Gibelin 05/09 : Add carbon spinup
!!      A.L. Gibelin 07/09 : Suppress RDK and transform GPP as a diagnostic
!!      D. Carrer    04/11 : Add FAPAR and effective LAI
!!      B. Decharme  09/2012 : suppress NWG_LAYER (parallelization problems)
!!      B. Decharme  09/12 : Carbon fluxes in diag_evap
!!      B. Decharme  09/12   New diag for DIF:
!!                           F2 stress
!!                           Root zone swi, wg and wgi
!!                           swi, wg and wgi comparable to ISBA-FR-DG2 and DG3 layers
!!                           active layer thickness over permafrost
!!                           frozen layer thickness over non-permafrost
!!      B. Decharme  06/13   All snow outputs noted SN
!!                           XTSRAD_NAT instead of XAVG_TSRAD
!!                           delete NWG_SIZE
!!                           water table depth
!!
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_type_snow
!
USE modd_data_cover_n, ONLY : data_cover_t
USE modd_surf_atm_n, ONLY : surf_atm_t
USE modd_diag_n, ONLY : diag_t, diag_np_t
USE modd_diag_misc_isba_n, ONLY : diag_misc_isba_t, diag_misc_isba_np_t
USE modd_isba_options_n, ONLY : isba_options_t
USE modd_isba_n, ONLY : isba_s_t, isba_k_t, isba_np_t
!
USE modd_xios, ONLY : lallow_add_dim
!
USE modd_surf_par,        ONLY :   nundef, xundef
!
USE modd_assim, ONLY : lassim, cassim_isba, nvar, cvar, nobstype, nboutput, cobs
!                                 
USE modd_agri,            ONLY :   lagrip
!
USE modi_init_io_surf_n
USE modi_write_surf
USE modi_write_field_2d_patch
USE modi_write_field_1d_patch
USE modi_end_io_surf_n
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
 CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: HSELECT
LOGICAL, INTENT(IN) :: OSNOWDIMNC
TYPE(surf_atm_t), INTENT(INOUT) :: U
LOGICAL, INTENT(IN) :: OPATCH_BUDGET
TYPE(diag_t), INTENT(INOUT) :: D
TYPE(diag_np_t), INTENT(INOUT) :: ND
TYPE(diag_misc_isba_t), INTENT(INOUT) :: DM
TYPE(diag_misc_isba_np_t), INTENT(INOUT) :: NDM
TYPE(isba_options_t), INTENT(INOUT) :: IO
TYPE(isba_s_t), INTENT(INOUT) :: S
TYPE(isba_k_t), INTENT(INOUT) :: K
TYPE(isba_np_t), INTENT(INOUT) :: NP
TYPE(surf_snow), INTENT(IN) :: TPSNOW
!
 CHARACTER(LEN=6),  INTENT(IN)  :: HPROGRAM ! program calling
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER           :: IRESP          ! IRESP  : return-code if a problem appears
 CHARACTER(LEN=1) :: YVAR, YOBS, YTIM
 CHARACTER(LEN=12) :: YRECFM         ! Name of the article to be read
 CHARACTER(LEN=100):: YCOMMENT       ! Comment string
 CHARACTER(LEN=2)  :: YLVL
 CHARACTER(LEN=20) :: YFORM
!
REAL, DIMENSION(SIZE(DM%XSWI,1)) :: ZMAX
INTEGER           :: JL, JJ, JVAR, JOBS, JP, JI, JT, JK, ISIZE
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!         Initialisation for IO
!
IF (lhook) CALL dr_hook('WRITE_DIAG_MISC_ISBA_N',0,zhook_handle)
!
IF ( dm%LPROSNOW ) THEN
  CALL init_io_surf_n(dtco, u, hprogram,'NATURE','ISBA  ','WRITE','ISBA_PROGNOSTIC.OUT.nc')
ELSE
  CALL init_io_surf_n(dtco, u, hprogram,'NATURE','ISBA  ','WRITE','ISBA_DIAGNOSTICS.OUT.nc')
ENDIF
!
!-------------------------------------------------------------------------------
!
IF (dm%LSURF_MISC_BUDGET) THEN
  !
  !*       2.     Miscellaneous fields :
  !
  !-------------------------------------------------------------------------------
  !
  !        2.1    Halstead coefficient
  !               --------------------
  !
  yrecfm='HV_ISBA'
  ycomment='Halstead coefficient averaged over tile nature (-)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XHV(:),iresp,hcomment=ycomment)
  !
  !        2.2    Snow fractions
  !               --------------
  !
  yrecfm='PSNG_ISBA'
  ycomment='snow fraction over ground averaged over tile nature (-)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XPSNG(:),iresp,hcomment=ycomment)
  !
  yrecfm='PSNV_ISBA'
  ycomment='snow fraction over vegetation averaged over tile nature (-)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XPSNV(:),iresp,hcomment=ycomment)
  !
  yrecfm='PSN_ISBA'
  ycomment='total snow fraction averaged over tile nature (-)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XPSN(:),iresp,hcomment=ycomment)
  !
  !        2.3    Total Albedo and surface temperature
  !               ------------------------------------
  !
  IF (tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN
    !        
    yrecfm='TS_ISBA'
    ycomment='total surface temperature (isba+snow) over tile nature'
    CALL write_surf(hselect, hprogram,yrecfm,d%XTS(:),iresp,hcomment=ycomment)
    !
    yrecfm='TSRAD_ISBA'
    ycomment='total radiative surface temperature (isba+snow) over tile nature'
    CALL write_surf(hselect, hprogram,yrecfm,s%XTSRAD_NAT(:),iresp,hcomment=ycomment)
    !
  END IF
  !
  !        2.4    Soil Wetness Index, Water content and active layer depth
  !               --------------------------------------------------------
  !  
  IF(io%CISBA=='DIF')THEN
    zmax(:) = 0.
    !
    DO jp = 1,io%NPATCH
      DO ji = 1,np%AL(jp)%NSIZE_P
        jj = np%AL(jp)%NR_P(ji)
        !
        IF (np%AL(jp)%NWG_LAYER(ji)/=nundef.AND.np%AL(jp)%NWG_LAYER(ji)>zmax(jj)) THEN
          zmax(jj) = np%AL(jp)%NWG_LAYER(ji)
        ENDIF
      ENDDO
    ENDDO
    !
    DO jj=1,SIZE(dm%XSWI,1)

      DO jl = 1,io%NGROUND_LAYER
        IF(jl>zmax(jj))THEN  
          dm%XSWI (jj,jl) = xundef
          dm%XTSWI(jj,jl) = xundef
        ENDIF
      ENDDO 

    ENDDO
  ENDIF         
  !
  DO jl=1,io%NGROUND_LAYER
    !
    WRITE(ylvl,'(I2)') jl
    !
    yrecfm='SWI'//adjustl(ylvl(:len_trim(ylvl)))
    yrecfm=yrecfm(:len_trim(yrecfm))//'_ISBA'
    yform='(A29,I1.1,A4)'
    IF (jl >= 10)  yform='(A29,I2.2,A4)'
    WRITE(ycomment,fmt=yform) 'soil wetness index for layer ',jl,' (-)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XSWI(:,jl),iresp,hcomment=ycomment)
    !
    yrecfm='TSWI'//adjustl(ylvl(:len_trim(ylvl)))
    yrecfm=yrecfm(:len_trim(yrecfm))//'_ISBA'
    yform='(A29,I1.1,A4)'
    IF (jl >= 10)  yform='(A29,I2.2,A4)'
    WRITE(ycomment,fmt=yform) 'total swi (liquid+solid) for layer ',jl,' (-)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XTSWI(:,jl),iresp,hcomment=ycomment)
    !
  END DO
  !
  yrecfm='SWI_T_ISBA'
  ycomment='soil wetness index over the soil column (-)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XSOIL_SWI(:),iresp,hcomment=ycomment)
  !
  yrecfm='TSWI_T_ISBA'
  ycomment='total soil wetness index over the soil column (-)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XSOIL_TSWI(:),iresp,hcomment=ycomment)
  !
  yrecfm='WGTOT_T_ISBA'
  ycomment='total water content (liquid+solid) over the soil column (kg/m2)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XSOIL_TWG(:),iresp,hcomment=ycomment)
  !
  yrecfm='WGI_T_ISBA'
  ycomment='total ice content (solid) over the soil column (kg/m2)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XSOIL_TWGI(:),iresp,hcomment=ycomment)
  !
  yrecfm='WGTOT_ISBA'
  ycomment='total volumetric water content (liquid+solid) over the soil column (m3/m3)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XSOIL_WG(:),iresp,hcomment=ycomment)
  !
  IF (.NOT.lallow_add_dim) THEN
    yrecfm='WGI_ISBA'
    ycomment='total volumetric ice content (solid) over the soil column (m3/m3)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XSOIL_WGI(:),iresp,hcomment=ycomment)
  ENDIF
  !
  IF(io%CISBA=='DIF') THEN
    !
    IF (dm%LSURF_MISC_DIF)THEN
      !
      yrecfm='TSWI_D2_ISBA'
      ycomment='total soil wetness index over comparable FR-DG2 reservoir (-)'
      CALL write_surf(hselect, hprogram,yrecfm,dm%XFRD2_TSWI(:),iresp,hcomment=ycomment)
      !
      yrecfm='WG_D2_ISBA'
      ycomment='liquid water content over comparable FR-DG2 reservoir (m3/m3)'
      CALL write_surf(hselect, hprogram,yrecfm,dm%XFRD2_TWG(:),iresp,hcomment=ycomment)
      !
      yrecfm='WGI_D2_ISBA'
      ycomment='ice content over comparable FR-DG2 reservoir (m3/m3)'
      CALL write_surf(hselect, hprogram,yrecfm,dm%XFRD2_TWGI(:),iresp,hcomment=ycomment)  
      !
      yrecfm='TSWI_D3_ISBA'
      ycomment='total soil wetness index over comparable FR-DG3 reservoir (-)'
      CALL write_surf(hselect, hprogram,yrecfm,dm%XFRD3_TSWI(:),iresp,hcomment=ycomment)
      !
      yrecfm='WG_D3_ISBA'
      ycomment='liquid water content over comparable FR-DG3 reservoir (m3/m3)'
      CALL write_surf(hselect, hprogram,yrecfm,dm%XFRD3_TWG(:),iresp,hcomment=ycomment)
      !
      yrecfm='WGI_D3_ISBA'
      ycomment='ice content over comparable FR-DG3 reservoir (m3/m3)'
      CALL write_surf(hselect, hprogram,yrecfm,dm%XFRD3_TWGI(:),iresp,hcomment=ycomment)  
      !
    ENDIF
    !
    yrecfm='ALT_ISBA'
    ycomment='active layer thickness over permafrost (m)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XALT(:),iresp,hcomment=ycomment)
    !
    yrecfm='FLT_ISBA'
    ycomment='frozen layer thickness over non-permafrost (m)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XFLT(:),iresp,hcomment=ycomment)
    !
  ENDIF
  !
  !        2.5    Snow outputs
  !               -------------
  !
  yrecfm='WSN_T_ISBA'
  ycomment='Total_snow_reservoir (kg/m2)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XTWSNOW(:),iresp,hcomment=ycomment)
  !
  yrecfm='DSN_T_ISBA'
  ycomment='Total_snow_depth (m)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XTDSNOW(:),iresp,hcomment=ycomment)
  !
  yrecfm='TSN_T_ISBA'
  ycomment='Total_snow_temperature (K)'
  CALL write_surf(hselect, hprogram,yrecfm,dm%XTTSNOW(:),iresp,hcomment=ycomment)
  !
  IF (tpsnow%SCHEME=='CRO' .AND. dm%LPROSNOW) THEN
    !
    ycomment='accumulated snow thickness for past 1 days'
    CALL write_surf(hselect,hprogram,'SD_1DY_ISBA',dm%XSNDPT_1DY(:),iresp,hcomment=ycomment) 
    !      
    ycomment='accumulated snow thickness for past 3 days'        
    CALL write_surf(hselect,hprogram,'SD_3DY_ISBA',dm%XSNDPT_3DY(:),iresp,hcomment=ycomment)                 
    !
    ycomment='accumulated snow thickness for past 5 days'
    CALL write_surf(hselect,hprogram,'SD_5DY_ISBA',dm%XSNDPT_5DY(:),iresp,hcomment=ycomment)                  
    !
    ycomment='accumulated snow thickness for past 7 days'  
    CALL write_surf(hselect,hprogram,'SD_7DY_ISBA',dm%XSNDPT_7DY(:),iresp,hcomment=ycomment)                  
    !
    ycomment='accumulated snow water equivalent for past 1 days'  
    CALL write_surf(hselect,hprogram,'SWE_1DY_ISBA',dm%XSNSWE_1DY(:),iresp,hcomment=ycomment) 
    !      
    ycomment='accumulated snow water equivalent for past 3 days'      
    CALL write_surf(hselect,hprogram,'SWE_3DY_ISBA',dm%XSNSWE_3DY(:),iresp,hcomment=ycomment)    
    !      
    ycomment='accumulated snow water equivalent for past 5 days'  
    CALL write_surf(hselect,hprogram,'SWE_5DY_ISBA',dm%XSNSWE_5DY(:),iresp,hcomment=ycomment)                  
    !
    ycomment='accumulated snow water equivalent for past 7 days'     
    CALL write_surf(hselect,hprogram,'SWE_7DY_ISBA',dm%XSNSWE_7DY(:),iresp,hcomment=ycomment)          
    !      
    ycomment='Penetration of ram resistance sensor (2 daN)'              
    CALL write_surf(hselect,hprogram,'RAMSOND_ISBA',dm%XSNRAM_SONDE(:),iresp,hcomment=ycomment)        
    !      
    ycomment='Thickness of wet snow at the top of the snowpack' 
    CALL write_surf(hselect,hprogram,'WET_TH_ISBA',dm%XSN_WETTHCKN(:),iresp,hcomment=ycomment)          
    !      
    ycomment='Thickness of refrozen snow at the top of the snowpack'
    CALL write_surf(hselect,hprogram,'REFRZTH_ISBA',dm%XSN_REFRZNTHCKN(:),iresp,hcomment=ycomment)     
    !      
  ENDIF
  !  
  !        2.6    SGH scheme
  !               ----------
  !
  IF(io%CRUNOFF=='SGH '.OR.io%CRUNOFF=='DT92')THEN     
    yrecfm='FSAT_ISBA'
    ycomment='Soil saturated fraction (-)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XFSAT(:),iresp,hcomment=ycomment)
  ENDIF
  !
  IF(io%CRAIN=='SGH ')THEN
    yrecfm='MUF_ISBA'
    ycomment='fraction of the grid cell reached by the rainfall (-)'
    CALL write_surf(hselect, hprogram,yrecfm,k%XMUF(:),iresp,hcomment=ycomment)
  ENDIF
  !
  !        2.7    Flooding scheme
  !               ---------------
  !
  IF(io%LFLOOD)THEN
    !
    yrecfm='FFG_ISBA'
    ycomment='flood fraction over ground averaged over tile nature (-)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XFFG(:),iresp,hcomment=ycomment)
    !
    yrecfm='FFV_ISBA'
    ycomment='flood fraction over vegetation averaged over tile nature (-)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XFFV(:),iresp,hcomment=ycomment)
    !
    yrecfm='FF_ISBA'
    ycomment='total flood fraction averaged over tile nature (-)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XFF(:),iresp,hcomment=ycomment)
    !
    yrecfm='FFLOOD_ISBA'
    ycomment='Grid-cell potential flood fraction (-)'
    CALL write_surf(hselect, hprogram,yrecfm,k%XFFLOOD(:),iresp,hcomment=ycomment)
    !
    yrecfm='PIFLOOD_ISBA'
    ycomment='Grid-cell Potential_floodplain_infiltration (kg/m2/s)'
    CALL write_surf(hselect, hprogram,yrecfm,k%XPIFLOOD(:),iresp,hcomment=ycomment)
    !
  ENDIF
  !
  !        2.8    Total LAI
  !               ---------
  !
  IF(io%CPHOTO/='NON'.OR.io%NPATCH>1)THEN        
    yrecfm='LAI_ISBA'
    ycomment='leaf area index (m2/m2)'
    CALL write_surf(hselect, hprogram,yrecfm,dm%XLAI(:),iresp,hcomment=ycomment)
  ENDIF
  !
  !        2.9    Water table depth
  !               -----------------
  !
  IF(io%LWTD)THEN
    !
    yrecfm='FWTD_ISBA'
    ycomment='grid-cell fraction of water table to rise'
    CALL write_surf(hselect, hprogram,yrecfm,k%XFWTD(:),iresp,hcomment=ycomment)
    !
    yrecfm='WTD_ISBA'
    ycomment='water table depth from RRM model or observation (m)'
    CALL write_surf(hselect, hprogram,yrecfm,k%XWTD(:),iresp,hcomment=ycomment)
    !
  ENDIF
  !*       3.     Miscellaneous fields for each patch :
  !               -------------------------------------
  !
  isize = u%NSIZE_NATURE
  !
  !----------------------------------------------------------------------------
  !User wants (or not) patch output
  IF(opatch_budget .AND. io%NPATCH>1)THEN
    !----------------------------------------------------------------------------
    !
    !        3.1    Soil Wetness Index and active layer depth
    !               -----------------------------------------   
    !
    DO jl=1,io%NGROUND_LAYER
      !
      WRITE(ylvl,'(I2.0)') jl
      !
      yrecfm='SWI'//trim(adjustl(ylvl(:)))//'_'
      yform='(A39,I1.1,A4)'
      IF (jl >= 10)  yform='(A39,I2.2,A4)'    
      WRITE(ycomment,fmt=yform) 'soil wetness index per patch for layer ',jl,' (-)'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XSWI(:,jl),isize,s%XWORK_WR)
      ENDDO      
      !
      yrecfm='TSWI'//trim(adjustl(ylvl(:)))//'_'
      yform='(A39,I1.1,A4)'
      IF (jl >= 10)  yform='(A39,I2.2,A4)'
      WRITE(ycomment,fmt=yform) 'total swi (liquid+solid) per patch for layer ',jl,' (-)'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XTSWI(:,jl),isize,s%XWORK_WR)
      ENDDO        
      !
    END DO
    !
    IF(io%CISBA=='DIF')THEN
      !
   
      yrecfm='ALT_'
      ycomment='active layer thickness over permafrost per patch (m)'
      WRITE(ycomment,fmt=yform) 'total swi (liquid+solid) per patch for layer ',jl,' (-)'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XALT(:),isize,s%XWORK_WR)
      ENDDO          
      !
      yrecfm='FLT_'
      ycomment='frozen layer thickness over non-permafrost per patch (m)'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XFLT(:),isize,s%XWORK_WR)
      ENDDO        
      !
    ENDIF
    !    
    !        3.2    Snow fractions
    !               --------------
    !  
    yrecfm='PSNG_'
    ycomment='snow fraction per patch over ground '
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XPSNG(:),isize,s%XWORK_WR)
    ENDDO  
    !   
    yrecfm='PSNV_'
    ycomment='snow fraction per patch over vegetation'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XPSNV(:),isize,s%XWORK_WR)
    ENDDO      
    !
    yrecfm='PSN_'
    ycomment='total snow fraction per patch'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XPSN(:),isize,s%XWORK_WR)
    ENDDO      
    !
    !        3.3    SGH scheme
    !               ----------
    !
    IF(io%CRUNOFF=='SGH '.OR.io%CRUNOFF=='DT92')THEN    
      yrecfm='FSAT_'
      ycomment='Soil saturated fraction per patch (-)'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XFSAT(:),isize,s%XWORK_WR)
      ENDDO        
    ENDIF
    !
    !        3.3    Flood fractions
    !               --------------
    !
    IF(io%LFLOOD)THEN
      !        
      yrecfm='FFG_'
      ycomment='flood fraction per patch over ground '
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XFFG(:),isize,s%XWORK_WR)
      ENDDO       
      !
      yrecfm='FFV_'
      ycomment='flood fraction per patch over vegetation'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XFFV(:),isize,s%XWORK_WR)
      ENDDO        
      !
      yrecfm='FF_'
      ycomment='total flood fraction per patch'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,ndm%AL(jp)%XFF(:),isize,s%XWORK_WR)
      ENDDO 
      !
    ENDIF
    !
    !        3.4    Total Albedo
    !               ------------
    !
    !
    IF (tpsnow%SCHEME=='3-L' .OR. tpsnow%SCHEME=='CRO') THEN
      !
      yrecfm='TS_'
      ycomment='total surface temperature (isba+snow) per patch'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,nd%AL(jp)%XTS(:),isize,s%XWORK_WR)
      ENDDO       
      !
      yrecfm='TSRAD_'
      ycomment='total radiative surface temperature (isba+snow) per patch'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
            np%AL(jp)%NR_P,nd%AL(jp)%XTSRAD(:),isize,s%XWORK_WR)
      ENDDO  
      !      
    ENDIF
    !
    !        3.5    Halstead coefficient
    !               --------------------
    !
    yrecfm='HV_'
    ycomment='Halstead coefficient per patch'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XHV(:),isize,s%XWORK_WR)
    ENDDO      
    !
    !        3.6  Snow outputs 
    !        -----------------
    !
    yrecfm='WSN_T_'
    ycomment='X_Y_WSNOW_TOT (kg/m2) per patch'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XTWSNOW(:),isize,s%XWORK_WR)
    ENDDO      
    !
    yrecfm='DSN_T_'
    ycomment='X_Y_DSNOW_TOT (m) per patch'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XTDSNOW(:),isize,s%XWORK_WR)
    ENDDO      
    !
    yrecfm='TSN_T_'
    ycomment='X_Y_TSNOW_TOT (k) per patch'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XTTSNOW(:),isize,s%XWORK_WR)
    ENDDO      
    !
    IF (tpsnow%SCHEME=='CRO' .AND. dm%LPROSNOW) THEN
      !
      ycomment='accumulated snow thickness for past 1 days per patch'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SD_1DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNDPT_1DY(:),isize,s%XWORK_WR)
      ENDDO        
      !
      ycomment=   'accumulated snow thickness for past 3 days per patch'   
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SD_3DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNDPT_3DY(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment=  'accumulated snow thickness for past 5 days per patch'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SD_5DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNDPT_5DY(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='accumulated snow thickness for past 7 days per patch'      
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SD_7DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNDPT_7DY(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='accumulated snow water equivalent for past 1 days per patch'   
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SWE_1DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNSWE_1DY(:),isize,s%XWORK_WR)
      ENDDO     
      !
      ycomment='accumulated snow water equivalent for past 3 days per patch'      
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SWE_3DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNSWE_3DY(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='accumulated snow water equivalent for past 5 days per patch'  
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SWE_5DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNSWE_5DY(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='accumulated snow water equivalent for past 7 days per patch'      
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'SWE_7DY_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNSWE_7DY(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='Penetration of ram resistance sensor (2 daN) per patch'               
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'RAMSOND_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNRAM_SONDE(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='Thickness of wet snow at the top of the snowpack per patch'    
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'WET_TH_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSN_WETTHCKN(:),isize,s%XWORK_WR)
      ENDDO      
      !
      ycomment='Thickness of refrozen snow at the top of the snowpack per patch'
      DO jp=1,io%NPATCH
        CALL write_field_1d_patch(hselect,hprogram,'REFRZTH_',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSN_REFRZNTHCKN(:),isize,s%XWORK_WR)
      ENDDO      
      !   
    ENDIF
    !
  ENDIF
  !
  IF((opatch_budget.AND. io%NPATCH>1).OR.io%NPATCH==1)THEN
    !
    IF ( osnowdimnc ) THEN
      !
      IF ( dm%LVOLUMETRIC_SNOWLIQ ) THEN
        ycomment='snow liquid water (kg m-3)'
      ELSE
        ycomment='snow liquid water (m)'
      ENDIF
      DO jp=1,io%NPATCH
        CALL write_field_2d_patch(hselect,hprogram,'SNOWLIQ',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWLIQ(:,:),isize,'snow_layer',s%XWSN_WR)
      ENDDO
      !
      ycomment='snow temperature (K)'
      DO jp=1,io%NPATCH
        CALL write_field_2d_patch(hselect,hprogram,'SNOWTEMP',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWTEMP(:,:),isize,'snow_layer',s%XWSN_WR)
      ENDDO
      !        
      ycomment=  'snow layer thickness'
      DO jp=1,io%NPATCH
        CALL write_field_2d_patch(hselect,hprogram,'SNOWDZ',ycomment,jp,&
              np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWDZ(:,:),isize,'snow_layer',s%XWSN_WR)
      ENDDO
      !      
      IF (tpsnow%SCHEME=='CRO' .AND. dm%LPROSNOW) THEN
        !
        ycomment=  'snow layer dendricity'
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWDEND',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWDEND(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO        
        !
        ycomment='snow layer sphericity'  
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWSPHER',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWSPHER(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO        
        !
        ycomment='snow layer grain size'
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWSIZE',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWSIZE(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO        
        !
        ycomment='snow layer specific surface area'
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWSSA',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWSSA(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO        
        !
        ycomment='snow layer grain type'
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWTYPE',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWTYPEMEPRA(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO        
        !
        ycomment='snow layer ram resistance'
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWRAM',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWRAM(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO        
        !
        ycomment='snow layer shear resistance'
        DO jp=1,io%NPATCH
          CALL write_field_2d_patch(hselect,hprogram,'SNOWSHEAR',ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWSHEAR(:,:),isize,'snow_layer',s%XWSN_WR)
        ENDDO
        !
      ENDIF
      !
    ELSE
      !
      DO jl=1,tpsnow%NLAYER
        !
       WRITE(ylvl,'(I2)') jl
        !
        yrecfm='SNOWLIQ'//adjustl(ylvl(:len_trim(ylvl)))
        yform='(A18,I1.1,A9)'
        IF (jl >= 10)  yform='(A18,I2.2,A9)'
        IF ( dm%LVOLUMETRIC_SNOWLIQ ) THEN
          WRITE(ycomment,yform) 'snow liquid water ',jl,' (kg m-3)'
        ELSE
          WRITE(ycomment,yform) 'snow liquid water ',jl,' (m)     '
        ENDIF        
        WRITE(ycomment,fmt=yform) 'snow liquid water',jl,' (m)'
        DO jp=1,io%NPATCH
          CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWLIQ(:,jl),isize,s%XWORK_WR)
        ENDDO          
        !
        yrecfm='SNOWTEMP'//adjustl(ylvl(:len_trim(ylvl)))
        yform='(A16,I1.1,A4)'
        IF (jl >= 10)  yform='(A16,I2.2,A4)'
        WRITE(ycomment,fmt=yform) 'snow temperature',jl,' (K)'
        DO jp=1,io%NPATCH
          CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
                np%AL(jp)%NR_P,ndm%AL(jp)%XSNOWTEMP(:,jl),isize,s%XWORK_WR)
        ENDDO
        !
      END DO
      !
    ENDIF
    !
  ENDIF
  !
  IF (lagrip) THEN
    !
    !        2.8    Irrigation threshold
    !               --------------------
    !
    yrecfm='IRRISEUIL'
    ycomment='irrigation threshold per patch'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XSEUIL(:),isize,s%XWORK_WR)
    ENDDO    
    !
  ENDIF
  !
  IF (io%LTR_ML) THEN
    !
    yrecfm='FAPAR'
    ycomment='FAPAR (-)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XFAPAR(:),isize,s%XWORK_WR)
    ENDDO
    !
    yrecfm='FAPIR'
    ycomment='FAPIR (-)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XFAPIR(:),isize,s%XWORK_WR)
    ENDDO
    !
    yrecfm='FAPAR_BS'
    ycomment='FAPAR_BS (-)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XFAPAR_BS(:),isize,s%XWORK_WR)
    ENDDO
    !
    yrecfm='FAPIR_BS'
    ycomment='FAPIR_BS (-)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XFAPIR_BS(:),isize,s%XWORK_WR)
    ENDDO
    !
    yrecfm='DFAPARC'
    ycomment='DFAPARC (-)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XDFAPARC(:),isize,s%XWORK_WR)
    ENDDO
    !
    yrecfm='DFAPIRC'
    ycomment='DFAPIRC (-)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XDFAPIRC(:),isize,s%XWORK_WR)
    ENDDO
    !
    yrecfm='DLAI_EFFC'
    ycomment='DLAI_EFFC (m2/m2)'
    DO jp=1,io%NPATCH
      CALL write_field_1d_patch(hselect,hprogram,yrecfm,ycomment,jp,&
          np%AL(jp)%NR_P,ndm%AL(jp)%XDLAI_EFFC(:),isize,s%XWORK_WR)
    ENDDO
    !
  ENDIF
  ! 
  IF (lassim .AND. cassim_isba=="EKF  ") THEN
    !
    CALL end_io_surf_n(hprogram)
    CALL init_io_surf_n(dtco, u, hprogram,'NATURE','ISBA  ','WRITE','ISBA_ANALYSIS.OUT.nc')
    !
    DO jvar = 1,nvar
      WRITE(yvar,fmt='(I1.1)') jvar
      yrecfm="ANA_INCR"//yvar
      !YCOMMENT="by patch"
      ycomment="Analysis increment for control variable "//trim(cvar(jvar))
      DO jp = 1,io%NPATCH
        CALL write_field_1d_patch(hselect, hprogram, yrecfm, ycomment, jp,&
                                np%AL(jp)%NR_P,np%AL(jp)%XINCR(1:np%AL(jp)%NSIZE_P,jvar),isize,s%XWORK_WR)
      ENDDO
      !
      WRITE(yvar,fmt='(I1.1)') jvar
      DO jt = 1,nboutput
        WRITE(ytim,fmt='(I1.1)') jt
        DO jobs = 1,nobstype
          WRITE(yobs,fmt='(I1.1)') jobs
          yrecfm="HO"//yobs//"_"//yvar//"_"//ytim
          !YCOMMENT="by patch"
          ycomment="Jacobian matrix for observation "//trim(cobs(jobs))//" and control variable "//trim(cvar(jvar))  
          jk = (jt-1)*nobstype + jobs
          DO jp = 1,io%NPATCH
            CALL write_field_1d_patch(hselect, hprogram, yrecfm, ycomment, jp,&
                                      np%AL(jp)%NR_P,np%AL(jp)%XHO(1:np%AL(jp)%NSIZE_P,jk,jvar),isize,s%XWORK_WR)
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    !
    DO jt = 1,nboutput        
      WRITE(ytim,fmt='(I1.1)') jt
      DO jobs = 1,nobstype
        WRITE(yobs,fmt='(I1.1)') jobs
        yrecfm="INNOV"//yobs//"_"//ytim
        !YCOMMENT="not by patch"
        ycomment="Innovation for observation "//trim(cobs(jobs))
        jk = (jt-1)*nobstype + jobs
        CALL write_surf(hselect,hprogram,yrecfm,s%XINNOV(:,jk),iresp,hcomment=ycomment)
      ENDDO
    ENDDO
    !
    DO jt = 1,nboutput
      WRITE(ytim,fmt='(I1.1)') jt
      DO jobs = 1,nobstype
        WRITE(yobs,fmt='(I1.1)') jobs
        yrecfm="RESID"//yobs//"_"//ytim
        !YCOMMENT="not by patch"
        ycomment="Residuals for observation "//trim(cobs(jobs))
        jk = (jt-1)*nobstype + jobs
        CALL write_surf(hselect,hprogram,yrecfm,s%XRESID(:,jk),iresp,hcomment=ycomment)
      ENDDO
    ENDDO
    !
  ENDIF
  !
ENDIF
!         End of IO
!
!
 CALL end_io_surf_n(hprogram)
!
IF (lhook) CALL dr_hook('WRITE_DIAG_MISC_ISBA_N',1,zhook_handle)
!
END SUBROUTINE write_diag_misc_isba_n