Newer
Older
!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 COMPUTE_ISBA_PARAMETERS (DTCO, OREAD_BUDGETC, UG, U, &
IO, DTI, SB, S, IG, K, NK, NIG, NP, NPE, &
NAG, NISS, ISS, NCHI, CHI, MGN, MSF, ID, &
GB, NGB, NDST, DST, SLT,BLOWSNW, SV, &
HPROGRAM,HINIT, &
KI,KMONTH,KDAY,KSV,KSW,HSV,PCO2,PRHOA, &
PEMIS,PTSRAD,PTSURF, HTEST, PMEGAN_FIELDS )
!#############################################################
!
!!**** *COMPUTE_ISBA_PARAMETERS_n* - routine to initialize ISBA
!!
!! PURPOSE
!! -------
!!
!!** METHOD
!! ------
!!
!! EXTERNAL
!! --------
!!
!!
!! IMPLICIT ARGUMENTS
!! ------------------
!!
!! REFERENCE
!! ---------
!!
!!
!! AUTHOR
!! ------
!!
!! MODIFICATIONS
!! -------------
!! Original 01/2004
!! Modified by P. Le Moigne (11/2004): miscellaneous diagnostics
!! Modified by P. Le Moigne (06/2006): seeding and irrigation
!! Modified by B. Decharme (2008) : SGH and Flooding scheme
!! Modified by B. Decharme (01/2009): optional deep soil temperature as in Arpege
!! Modified by R. Hamdi (01/2009): Cp and L
!! Modified by B. Decharme (06/2009): read topographic index statistics
!! Modified by P. Le Moigne (01/2009): Beljaars sso
!! Modified by B. Decharme (08/2009): Active Trip coupling variable if Earth System Model
!! A.L. Gibelin 04/09 : change BSLAI_NITRO initialisation
!! A.L. Gibelin 04/09 : modifications for CENTURY model
!! A.L. Gibelin 06/09 : soil carbon initialisation
!! Modified by B. Decharme (09/2012): Bug in exponential profile calculation with DIF
!! F. Bouttier 08/13 : apply random perturbation patterns for ensembles
!! B. Vincendon 03/14 : bug correction for CISBA=3L and CKSAT=EXP (TOPD coupling)
!! Modified by B. Decharme (04/2013): Subsurface runoff if SGH (DIF option only)
!! Delete CTOPREG (never used)
!! Delete NWG_LAYER_TOT, NWG_SIZE
!! water table / Surface coupling
!! P. Samuelsson 02/14 : MEB
!! B. Decharme 01/16 : Bug when vegetation veg, z0 and emis are imposed whith interactive vegetation
!! B. Decharme 10/2016 bug surface/groundwater coupling
!! P. Tulet 06/2016 : call init_megan for coupling megan with surfex
!! R. Schoetter 01/2018 add UTCI
!! A. Druel 02/2019 Deep change to accept patch duplication (for irrigation or agricultural pratices), mainly for ECOCLIMAP-SG
!! Séférian/Decharme 08/16 : include fire , river-carbon cycle module + new land-use change implementation
!! B. Decharme 02/2022 : split init_diag_isba routine
!!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_ISBA_OPTIONS_n, ONLY : ISBA_OPTIONS_t
USE MODD_ISBA_n, ONLY : ISBA_S_t, ISBA_P_t, ISBA_PE_t, ISBA_K_t, ISBA_NK_t, &
ISBA_NP_t, ISBA_NPE_t
USE MODD_DATA_ISBA_n, ONLY : DATA_ISBA_t
USE MODD_SFX_GRID_n, ONLY : GRID_t, GRID_NP_t
USE MODD_AGRI_n, ONLY : AGRI_t, AGRI_NP_t
USE MODD_SSO_n, ONLY : SSO_t, SSO_NP_t
USE MODD_CH_ISBA_n, ONLY : CH_ISBA_t, CH_ISBA_NP_t
USE MODD_CANOPY_n, ONLY : CANOPY_t
USE MODD_GR_BIOG_n, ONLY : GR_BIOG_t, GR_BIOG_NP_t
USE MODD_SURFEX_n, ONLY : ISBA_DIAG_t
USE MODD_DATA_COVER_n, ONLY : DATA_COVER_t
USE MODD_SURF_ATM_GRID_n, ONLY : SURF_ATM_GRID_t
USE MODD_SURF_ATM_n, ONLY : SURF_ATM_t
USE MODD_SLT_n, ONLY : SLT_t
USE MODD_SV_n, ONLY : SV_t
USE MODD_BLOWSNW_n, ONLY : BLOWSNW_t
USE MODD_MEGAN_n, ONLY : MEGAN_t
USE MODD_MEGAN_SURF_FIELDS_n, ONLY : MEGAN_SURF_FIELDS_t
!
USE MODD_SFX_OASIS, ONLY : LCPL_LAND, LCPL_FLOOD, LCPL_GW, LCPL_CALVING, &
LCPL_RIVCARB
!
!
#ifdef TOPD
USE MODD_DUMMY_EXP_PROFILE,ONLY : XC_DEPTH_RATIO
USE MODD_ASSIM, ONLY : CASSIM_ISBA, LASSIM
!
USE MODD_DEEPSOIL, ONLY : LPHYSDOMC, LDEEPSOIL, XTDEEP_CLI, XGAMMAT_CLI
USE MODD_AGRI, ONLY : LAGRIP, LIRRIGMODE, NVEG_IRR, XTHRESHOLD
!
USE MODD_SGH_PAR, ONLY : NDIMTAB, XICE_DEPH_MAX, XF_DECAY
!
USE MODD_DATA_COVER_PAR, ONLY : NVEGTYPE
USE MODD_SURF_PAR, ONLY : XUNDEF, NUNDEF, LEN_HREC
USE MODE_BLOWSNW_SEDIM_LKT1D
USE MODI_GET_1D_MASK
USE MODI_GET_Z0REL
USE MODI_GET_LUOUT
USE MODI_ABOR1_SFX
USE MODI_INIT_IO_SURF_n
USE MODI_ALLOCATE_PHYSIO
USE MODI_INIT_ISBA_MIXPAR
USE MODI_CONVERT_PATCH_ISBA
USE MODI_INIT_VEG_PGD_n
USE MODI_INIT_TOP
USE MODI_EXP_DECAY_SOIL_FR
USE MODI_CARBON_INIT
USE MODI_SOILTEMP_ARP_PAR
USE MODI_END_IO_SURF_n
!
USE MODI_INIT_CHEMICAL_n
USE MODI_OPEN_NAMELIST
USE MODI_CH_INIT_DEP_ISBA_n
USE MODI_CLOSE_NAMELIST
USE MODI_INIT_DST
USE MODI_INIT_SLT
USE MODI_INIT_DST_VEG
USE MODI_AVERAGED_ALBEDO_EMIS_ISBA
USE MODI_INIT_SURF_TOPD
USE MODI_ISBA_SOC_PARAMETERS
USE MODI_DIAG_ISBA_INIT_n
USE MODI_DIAG_MISC_ISBA_INIT_n
!
USE MODI_READ_AND_SEND_MPI
USE MODI_ISBA_TO_TOPD
USE MODI_OPEN_FILE
USE MODI_CLOSE_FILE
USE MODI_AV_PGD
USE MODI_SURF_PATCH
USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
USE PARKIND1 ,ONLY : JPRB
!
IMPLICIT NONE
!
!* 0.1 Declarations of arguments
! -------------------------
TYPE(DATA_COVER_t), INTENT(INOUT) :: DTCO
LOGICAL, INTENT(IN) :: OREAD_BUDGETC
TYPE(SURF_ATM_GRID_t), INTENT(INOUT) :: UG
TYPE(SURF_ATM_t), INTENT(INOUT) :: U
!
TYPE(ISBA_OPTIONS_t), INTENT(INOUT) :: IO
TYPE(DATA_ISBA_t), INTENT(INOUT) :: DTI
TYPE(CANOPY_t), INTENT(INOUT) :: SB
TYPE(ISBA_S_t), INTENT(INOUT) :: S
TYPE(GRID_t), INTENT(INOUT) :: IG
TYPE(ISBA_K_t), INTENT(INOUT) :: K
TYPE(ISBA_NK_t), INTENT(INOUT) :: NK
TYPE(GRID_NP_t), INTENT(INOUT) :: NIG
TYPE(ISBA_NP_t), INTENT(INOUT) :: NP
TYPE(ISBA_NPE_t), INTENT(INOUT) :: NPE
TYPE(AGRI_NP_t), INTENT(INOUT) :: NAG
TYPE(SSO_NP_t), INTENT(INOUT) :: NISS
TYPE(SSO_t), INTENT(INOUT) :: ISS
TYPE(CH_ISBA_NP_t), INTENT(INOUT) :: NCHI
TYPE(CH_ISBA_t), INTENT(INOUT) :: CHI
TYPE(MEGAN_t), INTENT(INOUT) :: MGN
TYPE(MEGAN_SURF_FIELDS_t), INTENT(INOUT) :: MSF
TYPE(ISBA_DIAG_t), INTENT(INOUT) :: ID
TYPE(GR_BIOG_t), INTENT(INOUT) :: GB
TYPE(GR_BIOG_NP_t), INTENT(INOUT) :: NGB
!
TYPE(DST_NP_t), INTENT(INOUT) :: NDST
TYPE(DST_t), INTENT(INOUT) :: DST
TYPE(SLT_t), INTENT(INOUT) :: SLT
TYPE(SV_t), INTENT(INOUT) :: SV
TYPE(BLOWSNW_t), INTENT(INOUT) :: BLOWSNW
CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes
CHARACTER(LEN=3), INTENT(IN) :: HINIT ! choice of fields to initialize
INTEGER, INTENT(IN) :: KI ! number of points
INTEGER, INTENT(IN) :: KMONTH ! Month
INTEGER, INTENT(IN) :: KDAY ! Day
INTEGER, INTENT(IN) :: KSV ! number of scalars
INTEGER, INTENT(IN) :: KSW ! number of short-wave spectral bands
CHARACTER(LEN=6), DIMENSION(KSV), INTENT(IN) :: HSV ! name of all scalar variables
REAL, DIMENSION(KI), INTENT(IN) :: PCO2 ! CO2 concentration (kg/m3)
REAL, DIMENSION(KI), INTENT(IN) :: PRHOA ! air density
REAL, DIMENSION(KI), INTENT(IN) :: PZENITH ! solar zenithal angle
REAL, DIMENSION(KSW), INTENT(IN) :: PSW_BANDS ! middle wavelength of each band
REAL, DIMENSION(KI,KSW),INTENT(OUT) :: PDIR_ALB ! direct albedo for each band
REAL, DIMENSION(KI,KSW),INTENT(OUT) :: PSCA_ALB ! diffuse albedo for each band
REAL, DIMENSION(KI), INTENT(OUT) :: PEMIS ! emissivity
REAL, DIMENSION(KI), INTENT(OUT) :: PTSRAD ! radiative temperature
REAL, DIMENSION(KI), INTENT(OUT) :: PTSURF ! surface effective temperature (K)
REAL, DIMENSION(KI,MSF%NMEGAN_NBR),INTENT(IN), OPTIONAL :: PMEGAN_FIELDS
!
CHARACTER(LEN=2), INTENT(IN) :: HTEST ! must be equal to 'OK'
!
!
!* 0.2 Declarations of local variables
! -------------------------------
!
TYPE(GRID_t), POINTER :: GK
TYPE(ISBA_P_t), POINTER :: PK
TYPE(ISBA_K_t), POINTER :: KK
TYPE(AGRI_t), POINTER :: AGK
TYPE(SSO_t), POINTER :: ISSK
TYPE(DST_t), POINTER :: DSTK
REAL, DIMENSION(U%NDIM_FULL) :: ZF_PARAM, ZC_DEPTH_RATIO
!
REAL, DIMENSION(KI) :: ZTSRAD_NAT !radiative temperature
REAL, DIMENSION(KI) :: ZTSURF_NAT !effective temperature
REAL, DIMENSION(KI) :: ZWG1 ! work array for surface water content
REAL, DIMENSION(KI,IO%NPATCH) :: ZTG1 ! work array for surface temperature
REAL, DIMENSION(KI,IO%NPATCH) :: ZF
REAL, DIMENSION(:,:), ALLOCATABLE :: ZWORK
INTEGER :: ICH ! unit of input chemistry file
INTEGER :: ILUOUT ! unit of output listing file
INTEGER :: IRESP ! return code
INTEGER :: IDECADE, IDECADE2 ! decade of simulation
INTEGER :: JP ! loop counter on tiles
INTEGER :: ISIZE_LMEB_PATCH ! Number of patches with MEB=true
LOGICAL :: GFIX ! update fixed fields
LOGICAL :: GTIME ! update time varing fields
LOGICAL :: GMEB ! update meb fields
LOGICAL :: GIRR ! update irrigation fields
LOGICAL :: GALBSOIL ! update bare soil albedo fields
LOGICAL :: GALBVEG ! update only vegetation albedo
LOGICAL :: GURBTREE ! update urban tree fields
!
LOGICAL :: GDIM, GCAS1, GCAS2, GCAS3
INTEGER :: JVEG, IVERSION, IBUGFIX, IMASK, JMAXLOC
!
CHARACTER(LEN=4) :: YLVL
CHARACTER(LEN=LEN_HREC) :: YRECFM
!
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
! Initialisation for IO
!
IF (LHOOK) CALL DR_HOOK('COMPUTE_ISBA_PARAMETERS',0,ZHOOK_HANDLE)
!
IF (HTEST/='OK') THEN
CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
! PART 1 : Things depending only on options and / or needed first
! --------------------------------------------------------------
!* Physiographic data fields from land cover:
! -----------------------------------------
!
IF (S%TTIME%TDATE%MONTH /= NUNDEF) THEN
IDECADE = 3 * ( S%TTIME%TDATE%MONTH - 1 ) + MIN(S%TTIME%TDATE%DAY-1,29) / 10 + 1
ELSE
IDECADE = 1
END IF
!
IDECADE2 = IDECADE
!
! concern DATA_ISBA, so no dependence on patches
CALL INIT_ISBA_MIXPAR(DTCO, DTI, IG%NDIM, IO, IDECADE, IDECADE2, S%XCOVER, S%LCOVER, 'NAT', U%LECOSG)
!
ISIZE_LMEB_PATCH=COUNT(IO%LMEB_PATCH(:))
IF (ISIZE_LMEB_PATCH>0) THEN
CALL FIX_MEB_VEG(DTI, IG%NDIM, IO%LMEB_PATCH, IO%NPATCH)
ENDIF
!
!
!* Soil carbon
! -----------
!
IF(HINIT=='ALL'.AND.(IO%CRESPSL=='CNT'.OR.IO%CRESPSL=='DIF').AND.IO%CPHOTO=='NCB')THEN
CALL CARBON_INIT
ENDIF
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
! PART 2 : Arrays of vegtypes & patches
! -------------------------------------
!
! We need XVEGTYPE, XPATCH and XVEGTYPE_PATCH with dimension "PATCH" for some
! cases: initialized here
!
! Vegtypes first
ALLOCATE(S%XVEGTYPE(KI,NVEGTYPE))
IF (DTI%LDATA_VEGTYPE) THEN
S%XVEGTYPE = DTI%XPAR_VEGTYPE
ELSE
!classical ecoclimap case
DO JVEG=1,NVEGTYPE
CALL AV_PGD(DTCO, S%XVEGTYPE(:,JVEG),S%XCOVER ,DTCO%XDATA_VEGTYPE(:,JVEG),'NAT','ARI',S%LCOVER)
END DO
ENDIF
!
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
! Take care about irrigation fraction if ECOSG and IRRIGATION
ALLOCATE(S%XVEGTYPE2(KI,NVEGTYPE+NVEG_IRR))
S%XVEGTYPE2(:,:) = 0.
IF (DTI%LDATA_MIXPAR) THEN ! Probably sufficient case.
ALLOCATE(DTI%XPAR_VEGTYPE2(SIZE(DTI%XPAR_VEGTYPE,1),SIZE(DTI%XPAR_VEGTYPE,2)+NVEG_IRR)) ! because if .NOT.(LIRRIGMODE .OR. LAGRIP),
DTI%XPAR_VEGTYPE2(:,:) = 0. ! DTI%LDATA_MIXPAR automaticaly = .TRUE.
ELSE
ALLOCATE(DTI%XPAR_VEGTYPE2(1,1))
ENDIF
!
IF ( U%LECOSG .AND. ( LIRRIGMODE .OR. LAGRIP ) ) THEN
!
IF ( NVEG_IRR == 0 ) CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS - WITH LECOSG AND (LIRRIGMODE OR LAGRIP), NVEG_IRR HAVE TO BE >0')
!
IF ( NVEG_IRR /= SIZE(DTI%NPAR_VEG_IRR_USE)) &
CALL ABOR1_SFX('PCOMPUTE_ISBA_PARAMETER: WHEN LIRRIGMODE or LAGRIP+LECOSG, NVEG_IRR AND NPAR_VEG_IRR_USE HAVE TO BE EQUAL')
!
! IF ECOSG and IRRIGATION, compute fraction of non irrigated vegtype + irrigated vegtype
DO JVEG=1,NVEGTYPE+NVEG_IRR
IF ( JVEG <= NVEGTYPE ) THEN
IF ( ANY(DTI%NPAR_VEG_IRR_USE(:)==JVEG ) ) THEN
WHERE ( DTI%XPAR_IRRIGFRAC(:,JVEG) /= XUNDEF .AND. DTI%XPAR_IRRIGFRAC(:,JVEG) /= 0. )
S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,JVEG) * (1 - DTI%XPAR_IRRIGFRAC(:,JVEG) )
ELSEWHERE
S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,JVEG)
ENDWHERE
IF (DTI%LDATA_MIXPAR) THEN
WHERE ( DTI%XPAR_IRRIGFRAC(:,JVEG) /= XUNDEF .AND. DTI%XPAR_IRRIGFRAC(:,JVEG) /= 0. )
DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,JVEG) * (1 - DTI%XPAR_IRRIGFRAC(:,JVEG) )
ELSEWHERE
DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,JVEG)
ENDWHERE
ENDIF
ELSE
S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,JVEG)
IF ( DTI%LDATA_MIXPAR ) DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,JVEG)
ENDIF
ELSE
WHERE ( DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= XUNDEF .AND. &
DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= 0. )
S%XVEGTYPE2(:,JVEG) = S%XVEGTYPE(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) &
* DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE))
ELSEWHERE
S%XVEGTYPE2(:,JVEG) = 0.
ENDWHERE
IF (DTI%LDATA_MIXPAR) THEN
WHERE ( DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= XUNDEF .AND. &
DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) /= 0. )
DTI%XPAR_VEGTYPE2(:,JVEG) = DTI%XPAR_VEGTYPE(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE)) &
* DTI%XPAR_IRRIGFRAC(:,DTI%NPAR_VEG_IRR_USE(JVEG-NVEGTYPE))
ELSEWHERE
DTI%XPAR_VEGTYPE2(:,JVEG) = 0.
ENDWHERE
ENDIF
ENDIF
ENDDO
!
ELSE
!
IF ( NVEG_IRR /= 0 ) CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS - WITHOUT LECOSG AND (LIRRIGMODE OR LAGRIP), NVEG_IRR HAVE TO BE =0')
!
S%XVEGTYPE2(:,:) = S%XVEGTYPE(:,:)
IF (DTI%LDATA_MIXPAR) DTI%XPAR_VEGTYPE2(:,:) = DTI%XPAR_VEGTYPE(:,:)
!
ENDIF
!
!
!
ALLOCATE(S%XPATCH(KI, IO%NPATCH))
ALLOCATE(S%XVEGTYPE_PATCH(KI, NVEGTYPE+NVEG_IRR, IO%NPATCH))
CALL SURF_PATCH(IO%NPATCH,S%XVEGTYPE2,DTI%NPAR_VEG_IRR_USE,S%XPATCH,S%XVEGTYPE_PATCH)
!
!
! removing little fractions of patches must be done of the XPATCH with dimension
! "PATCH"
IF (IO%XRM_PATCH/=0.) THEN
!
WRITE(ILUOUT,*) " REMOVE PATCH below 5 % add to dominant patch "
! remove small fraction of PATCHES and add to MAIN PATCH
DO JI = 1,KI
!1) find most present patch maximum value
JMAXLOC = MAXVAL(MAXLOC(S%XPATCH(JI,:)))
!2) FIND small value of cover
DO JP = 1,IO%NPATCH
IF ( S%XPATCH(JI,JP)<IO%XRM_PATCH ) THEN
S%XPATCH(JI,JMAXLOC) = S%XPATCH(JI,JMAXLOC) + S%XPATCH(JI,JP)
S%XPATCH(JI,JP) = 0.0
ENDIF
ENDDO
ENDDO
!
ENDIF
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
! PART 3 : Loop on patches for general initialization
! --------------------------------------------------
! loop on patches
DO JP = 1, IO%NPATCH
!
KK => NK%AL(JP)
PK => NP%AL(JP)
PEK => NPE%AL(JP)
AGK => NAG%AL(JP)
ISSK => NISS%AL(JP)
! dimension of the patch
PK%NSIZE_P = COUNT(S%XPATCH(:,JP) > 0.0)
!
! mask of the patch in tile nature
ALLOCATE(PK%NR_P (PK%NSIZE_P))
CALL GET_1D_MASK(PK%NSIZE_P, KI, S%XPATCH(:,JP), PK%NR_P)
!
! the array of vegtypes, patches and vegtypes by patches reduced on this patches
ALLOCATE(KK%XVEGTYPE(PK%NSIZE_P,NVEGTYPE))
CALL PACK_SAME_RANK(PK%NR_P,S%XVEGTYPE,KK%XVEGTYPE)
!
ALLOCATE(PK%XPATCH(PK%NSIZE_P))
ALLOCATE(PK%XVEGTYPE_PATCH (PK%NSIZE_P,NVEGTYPE+NVEG_IRR))
CALL PACK_SAME_RANK(PK%NR_P,S%XPATCH(:,JP),PK%XPATCH)
CALL PACK_SAME_RANK(PK%NR_P,S%XVEGTYPE_PATCH(:,:,JP),PK%XVEGTYPE_PATCH)
!
!
! soon needed packed fields
!
IF (IO%LPERM) THEN
ALLOCATE(KK%XPERM(PK%NSIZE_P))
CALL PACK_SAME_RANK(PK%NR_P, K%XPERM, KK%XPERM)
ELSE
ALLOCATE(KK%XPERM(0))
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
!
!
ALLOCATE(KK%XSAND(PK%NSIZE_P,IO%NGROUND_LAYER))
ALLOCATE(KK%XCLAY(PK%NSIZE_P,IO%NGROUND_LAYER))
!
ALLOCATE(ISSK%XAOSIP(PK%NSIZE_P))
ALLOCATE(ISSK%XAOSIM(PK%NSIZE_P))
ALLOCATE(ISSK%XAOSJP(PK%NSIZE_P))
ALLOCATE(ISSK%XAOSJM(PK%NSIZE_P))
ALLOCATE(ISSK%XHO2IP(PK%NSIZE_P))
ALLOCATE(ISSK%XHO2IM(PK%NSIZE_P))
ALLOCATE(ISSK%XHO2JP(PK%NSIZE_P))
ALLOCATE(ISSK%XHO2JM(PK%NSIZE_P))
!
!
CALL PACK_SAME_RANK(PK%NR_P, K%XSAND, KK%XSAND)
CALL PACK_SAME_RANK(PK%NR_P, K%XCLAY, KK%XCLAY)
!
CALL PACK_SAME_RANK(PK%NR_P,ISS%XAOSIP,ISSK%XAOSIP)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XAOSIM,ISSK%XAOSIM)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XAOSJP,ISSK%XAOSJP)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XAOSJM,ISSK%XAOSJM)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XHO2IP,ISSK%XHO2IP)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XHO2IM,ISSK%XHO2IM)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XHO2JP,ISSK%XHO2JP)
CALL PACK_SAME_RANK(PK%NR_P,ISS%XHO2JM,ISSK%XHO2JM)
!
!-------------------------------------------------------------------------------
!
!* 2.5 Physiographic fields
! --------------------
!
CALL ALLOCATE_PHYSIO(IO, KK, PK, PEK, NVEGTYPE, U%LECOSG)
!-------------------------------------------------------------------------------
!
! Initialize general parameters
!
GFIX = .TRUE. ! update fixed fields
GTIME = .TRUE. ! update time varing fields
GMEB = .TRUE. ! update meb fields
GIRR = .TRUE. ! update irrigation fields
GURBTREE = .FALSE. ! update urban tree fields
!
CALL CONVERT_PATCH_ISBA(DTCO, DTI, IO, KMONTH, KDAY, IDECADE, IDECADE2, S%XCOVER, S%LCOVER, &
LAGRIP, U%LECOSG, LIRRIGMODE, 'NAT', JP, KK, PK, PEK, GFIX, &
GTIME, GMEB, GIRR, GURBTREE, PSOILGRID=IO%XSOILGRID, PPERM=KK%XPERM )
!
!-------------------------------------------------------------------------------
!
! Initialize vegetation albedo only
!
GALBSOIL = .FALSE. ! update bare soil albedo fields
GALBVEG = .TRUE. ! update only vegetation albedo
!
CALL CONVERT_PATCH_ALB_ISBA(DTCO, DTI, IO, IDECADE, IDECADE2, S%XCOVER, S%LCOVER, &
'NAT',JP, KK, PK, PEK, GALBSOIL, GALBVEG )
!
!-------------------------------------------------------------------------------
!
! in init_veg_pgd_n, things needed also by garden and greenroof
CALL INIT_VEG_PGD_n(ISSK, DTI, IO, S, K, KK, PK, PEK, AGK, KI, &
HPROGRAM, 'NATURE', ILUOUT, PK%NSIZE_P, S%TTIME%TDATE%MONTH, &
.TRUE., LDEEPSOIL, LPHYSDOMC, .TRUE., .TRUE., XTDEEP_CLI, XGAMMAT_CLI, &
LAGRIP,LIRRIGMODE,XTHRESHOLD, HINIT, PCO2, PRHOA)
!
!-------------------------------------------------------------------------------
!
! Other fields needed to be initialized for isba only
!
!Rainfall spatial distribution
!CRAIN used in HYDRO_VEG and HYDRO_SGH and VEG_SGH_UPDATE
IF(IO%CRAIN=='SGH')THEN
ALLOCATE(KK%XMUF(PK%NSIZE_P))
KK%XMUF(:)=0.0
!
ALLOCATE(KK%XFSAT(PK%NSIZE_P))
KK%XFSAT(:) = 0.0
!
! * Initialize flood scheme :
!
ALLOCATE(KK%XFFLOOD (PK%NSIZE_P))
ALLOCATE(KK%XPIFLOOD(PK%NSIZE_P))
ALLOCATE(KK%XFF (PK%NSIZE_P))
ALLOCATE(KK%XFFG (PK%NSIZE_P))
ALLOCATE(KK%XFFV (PK%NSIZE_P))
ALLOCATE(KK%XFFROZEN(PK%NSIZE_P))
ALLOCATE(KK%XALBF (PK%NSIZE_P))
ALLOCATE(KK%XEMISF (PK%NSIZE_P))
KK%XFFLOOD = 0.0
KK%XPIFLOOD = 0.0
KK%XFF = 0.0
KK%XFFG = 0.0
KK%XFFV = 0.0
KK%XFFROZEN = 0.0
KK%XALBF = 0.0
KK%XEMISF = 0.0
!
ENDDO
IF (DTI%LDATA_CONDSAT) DEALLOCATE(DTI%XPAR_CONDSAT)
!
!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------
!
! PART 4 : Initialization not depending on patches
! ------------------------------------------------
!
! Fields needed also unpacked
!
IF(IO%CRAIN=='SGH')THEN
ALLOCATE(K%XMUF(KI))
K%XMUF(:)=0.0
!-------------------------------------------------------------------------------
! PART 5: Initialize Chemical Deposition
! -----------------------------------
! 3.1 Chemical gazes
! --------------
!* for the time being, chemistry on vegetation works only for
! ISBA on nature tile (not for gardens), because subroutine INIT_CHEMICAL_n
! contains explicitely modules from ISBAn. It should be cleaned in a future
! version.
CALL INIT_CHEMICAL_n(ILUOUT, KSV, HSV, CHI%SVI, CHI%SLTI, CHI%DSTI, &
CHI%CCH_NAMES, CHI%CAER_NAMES, &
HDSTNAMES=CHI%CDSTNAMES, HSLTNAMES=CHI%CSLTNAMES, &
HSNWNAMES=CHI%CSNWNAMES )
IF (KSV /= 0) THEN
!
IF (CHI%SVI%NBEQ > 0) THEN
!* for the time being, chemistry deposition on vegetation works only for
! ISBA on nature tile (not for gardens), because subroutine CH_INIT_DEP_ISBA_n
! contains explicitely modules from ISBAn. It should be cleaned in a future
! version.
CALL OPEN_NAMELIST(HPROGRAM, ICH, HFILE=CHI%CCHEM_SURF_FILE)
CALL CH_INIT_DEP_ISBA_n(CHI, NCHI, NP, DTCO, IO%NPATCH, S%LCOVER, S%XCOVER, ICH, ILUOUT, KI, &
DTI%XPAR_IRRIGFRAC, DTI%NPAR_VEG_IRR_USE)
CALL CLOSE_NAMELIST(HPROGRAM, ICH)
END IF
!
DO JP = 1,IO%NPATCH
!
DSTK => NDST%AL(JP)
!
IF (CHI%SVI%NDSTEQ >=1) THEN
!

RODIER Quentin
committed
ALLOCATE (DSTK%XSFDST (PK%NSIZE_P, CHI%SVI%NDSTEQ)) !Output array
ALLOCATE (DSTK%XSFDSTM(PK%NSIZE_P, CHI%SVI%NDSTEQ)) !Output array
DSTK%XSFDSTM(:,:) = 0.
CALL INIT_DST_VEG(DSTK, U, HPROGRAM, PK%NSIZE_P, PK%NR_P, PK%XVEGTYPE_PATCH)
ELSE
ALLOCATE(DSTK%XSFDST (0,0))
ALLOCATE(DSTK%XSFDSTM(0,0))
END IF
!
ENDDO
!
!
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
IF (CHI%SVI%NSNWEQ >=1) THEN
ALLOCATE (BLOWSNW%XSNW_FSED(KI,CHI%SVI%NSNWEQ+1)) !Output array
ALLOCATE (BLOWSNW%XSNW_FTURB(KI,CHI%SVI%NSNWEQ+1)) !Output array
ALLOCATE (BLOWSNW%XSNW_FNET(KI,CHI%SVI%NSNWEQ+1)) !Output array
ALLOCATE (BLOWSNW%XSNW_FSALT(KI,CHI%SVI%NSNWEQ+1)) !Output array
ALLOCATE (BLOWSNW%XSFSNW(KI,CHI%SVI%NSNWEQ+1)) !Output array
ALLOCATE (BLOWSNW%XSNW_SUBL(KI,CHI%SVI%NSNWEQ+1)) !Output array
BLOWSNW%XSNW_FSED (:,:) = 0.
BLOWSNW%XSNW_FTURB(:,:) = 0.
BLOWSNW%XSNW_FNET (:,:) = 0.
BLOWSNW%XSNW_FSALT(:,:) = 0.
BLOWSNW%XSNW_SUBL (:,:) = 0.
BLOWSNW%XSFSNW (:,:) = 0.
!Read in look up tables of snow particles properties
!No arguments, all look up tables are defined in module
!mode_snowdrift_sedim_lkt
CALL BLOWSNW_SEDIM_LKT1D_SET
ELSE
ALLOCATE(BLOWSNW%XSNW_FSED(0,0))
ALLOCATE(BLOWSNW%XSNW_FTURB(0,0))
ALLOCATE(BLOWSNW%XSNW_FSALT(0,0))
ALLOCATE(BLOWSNW%XSNW_FNET(0,0))
ALLOCATE(BLOWSNW%XSNW_SUBL(0,0))
ALLOCATE(BLOWSNW%XSFSNW(0,0))
END IF
!-------------------------------------------------------------------------------
! PART 6: Specific options
! --------------------------
!6.A. DIF option :
!---------------
! Anisotropy coeficient for hydraulic conductivity for topmodel drainage (Fan et al. 2006)
! Soil organic matter effect and/or Exponential decay for DIF option
! Must be call before INIT_TOP
IF(IO%CISBA=='DIF' .AND. IO%CKSAT=='SGH') THEN
!
WRITE(ILUOUT,*)'THE KSAT EXP PROFILE WITH ISBA-DF IS NOT PHYSIC AND HAS BEEN REMOVED FOR NOW'
WRITE(ILUOUT,*)'A NEW PHYSICAL APPROACH WILL BE DEVELLOPED ACCOUNTING FOR COMPACTION IN ALL '
WRITE(ILUOUT,*)'HYDRODYNAMIC PARAMETERS (WSAT, PSISAT, KSAT, B) AND NOT ONLY IN KSAT '
CALL ABOR1_SFX('CKSAT=SGH is not physic with ISBA-DF and has been removed for now')
!
ENDIF
!
IF(IO%CISBA=='DIF' .AND. IO%LSOC)THEN
!
IF(.NOT.IO%LSOCP)THEN
WRITE(ILUOUT,*)'LSOC = T can be activated only if SOC data given in PGD fields'
CALL ABOR1_SFX('LSOC = T can be activated only if SOC data given in PGD fields')
!
ALLOCATE(S%XFRACSOC(KI,IO%NGROUND_LAYER))
CALL ISBA_SOC_PARAMETERS(IO%CRUNOFF, S%XSOC, K, NP, S%XFRACSOC, &
K%XWSAT, K%XWFC, K%XWWILT, IO%NPATCH )
!
ELSE
ALLOCATE(S%XFRACSOC(0,0))
ENDIF
!
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
!6.B. Topmodel
!--------------
!
ZF (:,:) = XUNDEF
ZM (:) = XUNDEF
!
!CRUNOFF used in hydro_sgh and isba_sgh_update
IF( IO%CRUNOFF=='SGH '.AND. HINIT/='PRE' .AND. .NOT.LASSIM ) THEN
!
! Subsurface flow by layer (m/s)
DO JP = 1,IO%NPATCH
PK => NP%AL(JP)
IF(IO%CISBA=='DIF') THEN
ALLOCATE(PK%XTOPQS(PK%NSIZE_P,IO%NGROUND_LAYER))
PK%XTOPQS(:,:) = 0.0
ELSE
ALLOCATE(PK%XTOPQS(0,0))
ENDIF
ENDDO
!
ALLOCATE(S%XTAB_FSAT(KI,NDIMTAB))
ALLOCATE(S%XTAB_WTOP(KI,NDIMTAB))
ALLOCATE(S%XTAB_QTOP(KI,NDIMTAB))
S%XTAB_FSAT(:,:) = 0.0
S%XTAB_WTOP(:,:) = 0.0
S%XTAB_QTOP(:,:) = 0.0
!
WHERE(K%XCLAY(:,1)==XUNDEF.AND.S%XTI_MEAN(:)/=XUNDEF) S%XTI_MEAN(:)=XUNDEF
CALL INIT_TOP(IO, S, K, NK, NP, ILUOUT, ZM )
!
!
DO JP = 1,IO%NPATCH
PK => NP%AL(JP)
ALLOCATE(PK%XTOPQS(0,0))
ENDDO
!
ALLOCATE(S%XTAB_FSAT(0,0))
ALLOCATE(S%XTAB_WTOP(0,0))
ALLOCATE(S%XTAB_QTOP(0,0))
!
!Exponential decay for ISBA-FR option
!CKSAT used in hydro_soil.F90 and soil.F90
IF ( IO%CISBA/='DIF' .AND. HINIT/='PRE' .AND. .NOT.LASSIM ) THEN
GCAS1 = (IO%CKSAT=='EXP' .AND. IO%CISBA=='3-L')
GCAS2 = (IO%CKSAT=='SGH')
GCAS3 = (HPROGRAM/='AROME ' .AND. HPROGRAM/='MESONH ')
!
IF ( GCAS1 .OR. GCAS2 ) THEN
ALLOCATE(S%XF_PARAM (KI))
S%XF_PARAM(:) = XUNDEF
CALL OPEN_FILE('ASCII ',NUNIT,HFILE='carte_f_dc.txt',HFORM='FORMATTED',HACTION='READ ')
DO JI = 1,U%NDIM_FULL
READ(NUNIT,*) ZF_PARAM(JI), ZC_DEPTH_RATIO(JI)
CALL READ_AND_SEND_MPI(ZF_PARAM,S%XF_PARAM,U%NR_NATURE)
IF (.NOT.ALLOCATED(XC_DEPTH_RATIO)) ALLOCATE(XC_DEPTH_RATIO (KI))
XC_DEPTH_RATIO(:) = XUNDEF
CALL READ_AND_SEND_MPI(ZC_DEPTH_RATIO,XC_DEPTH_RATIO,U%NR_NATURE)
#endif
WRITE(ILUOUT,*) "COMPUTE_ISBA_PARAMETERS: WITH CKSAT=EXP, IN NOT OFFLINE "//&
"MODE, TOPMODEL FILE FOR F_PARAM IS NOT READ "
ENDIF
!
! definition of ZF functions of options
!
! Exponential decay factor calculate using soil properties
! (eq. 11, Decharme et al., J. Hydrometeor, 2006)
DO JP = 1,IO%NPATCH
PK => NP%AL(JP)
!
DO JI = 1,PK%NSIZE_P
IMASK = PK%NR_P(JI)
IF ( GCAS2 .AND. IO%CRUNOFF=='SGH' .AND. ZM(IMASK)/=XUNDEF ) THEN
ZF(JI,JP) = (K%XWSAT(IMASK,1)-K%XWD0(IMASK,1)) / ZM(IMASK)
ELSEIF ( GCAS1 ) THEN
ZF(JI,JP) = S%XF_PARAM(IMASK)
ENDIF
ENDDO
ENDDO
!
DO JP = 1,IO%NPATCH
PK => NP%AL(JP)
!
WHERE ( ZF(1:PK%NSIZE_P,JP)==XUNDEF.AND.PK%XDG(:,2)/=XUNDEF )
ZF(1:PK%NSIZE_P,JP) = 4.0/PK%XDG(:,2)
ZF(1:PK%NSIZE_P,JP) = MIN(ZF(1:PK%NSIZE_P,JP),XF_DECAY)
!
ZC_DEPTH_RATIO(1:PK%NSIZE_P) = 1.
#ifdef TOPD
IF (ALLOCATED(XC_DEPTH_RATIO)) THEN
CALL PACK_SAME_RANK(PK%NR_P,XC_DEPTH_RATIO,ZC_DEPTH_RATIO(1:PK%NSIZE_P))
ENDIF
#endif
CALL EXP_DECAY_SOIL_FR(IO%CISBA, ZF(1:PK%NSIZE_P,JP), PK, ZC_DEPTH_RATIO(1:PK%NSIZE_P))
IF ( GCAS2 ) THEN
!
DO JI = 1,NP%AL(1)%NSIZE_P
IMASK = NP%AL(1)%NR_P(JI)
S%XF_PARAM(IMASK) = ZF(JI,1)
ENDDO
!
ENDIF
!
! 6.C. Initialize required coupling fields :
!-------------------------------------------
IO%LCPL_RRM = .FALSE.
IO%LFLOOD = .FALSE.
IO%LWTD = .FALSE.
! * River coupling
!
IO%LCPL_RRM = .TRUE.
ALLOCATE(S%XCPL_DRAIN (KI))
ALLOCATE(S%XCPL_RUNOFF(KI))
S%XCPL_DRAIN (:) = 0.0
S%XCPL_RUNOFF(:) = 0.0
S%XCPL_TWS (:) = 0.0
!
! * Calving coupling
!
IF(IO%LGLACIER)THEN
ALLOCATE(S%XCPL_ICEFLUX(KI))
S%XCPL_ICEFLUX(:) = 0.0
!
! * Groundwater coupling
!
IF(LCPL_GW)THEN
IO%LWTD = .TRUE.
ENDIF
!
! * Floodplains coupling
IO%LFLOOD = .TRUE.
ALLOCATE(S%XCPL_EFLOOD(KI))
ALLOCATE(S%XCPL_PFLOOD(KI))
ALLOCATE(S%XCPL_IFLOOD(KI))
S%XCPL_EFLOOD(:)= 0.0
S%XCPL_PFLOOD(:)= 0.0
S%XCPL_IFLOOD(:)= 0.0
ALLOCATE(S%XCPL_EFLOOD(0))
ALLOCATE(S%XCPL_PFLOOD(0))
ALLOCATE(S%XCPL_IFLOOD(0))
! * Riverine-carbon cycle coupling
!
IF(IO%LCLEACH)THEN
ALLOCATE(S%XCPL_DOCFLUX(KI))
S%XCPL_DOCFLUX(:) = 0.0
ELSE
ALLOCATE(S%XCPL_DOCFLUX(0))
ENDIF
!
ALLOCATE(S%XCPL_RUNOFF (0))
ALLOCATE(S%XCPL_DRAIN (0))
ALLOCATE(S%XCPL_ICEFLUX (0))
ALLOCATE(S%XCPL_EFLOOD (0))
ALLOCATE(S%XCPL_PFLOOD (0))
ALLOCATE(S%XCPL_IFLOOD (0))
ALLOCATE(S%XCPL_DOCFLUX (0))
ALLOCATE(S%XCPL_TWS (0))
ALLOCATE(K%XFWTD(KI))
ALLOCATE(K%XWTD (KI))
K%XFWTD(:) = 0.0
K%XWTD (:) = XUNDEF
IF(LCPL_FLOOD)THEN
ALLOCATE(K%XFFLOOD (KI))
ALLOCATE(K%XPIFLOOD(KI))
K%XFFLOOD (:) = 0.0
K%XPIFLOOD(:) = 0.0
!
ELSE
!
ALLOCATE(K%XFFLOOD (0))
ALLOCATE(K%XPIFLOOD(0))
!
ENDIF
!
!
ALLOCATE(K%XFWTD(0))
ALLOCATE(K%XWTD (0))
ALLOCATE(K%XFFLOOD (0))
ALLOCATE(K%XPIFLOOD(0))
!
ENDIF
!
! * Check some key :
!
IF(LCPL_CALVING)THEN
IF(.NOT.IO%LGLACIER)THEN
CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS: LGLACIER MUST BE ACTIVATED IF LCPL_CALVING')
ENDIF
IF(LCPL_RIVCARB)THEN
IF(.NOT.IO%LCLEACH)THEN
CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS: LCLEACH MUST BE ACTIVATED IF LCPL_RIVCARB')
ENDIF
ENDIF
!
IF(IO%LCLEACH.AND.IO%CPHOTO/='NCB')THEN
WRITE(ILUOUT,*)'COMPUTE_ISBA_PARAMETERS: Riverine - Carbon cycle coupling '
WRITE(ILUOUT,*)'COMPUTE_ISBA_PARAMETERS: Please check your prep namelist where CPHOTO must be NCB '
CALL ABOR1_SFX('COMPUTE_ISBA_PARAMETERS: Withe LCLEACH, NCB vegetation module is required')
ENDIF
!
!-------------------------------------------------------------------------------
!
!* 6.D. ISBA time-varying deep force-restore temperature initialization
! --------------------------------------------------------------------
!
!-------------------------------------------------------------------------------
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!-------------------------------------------------------------------------------
!
! PART 7: We packed needed fields and free unless ones
! -----------------------------------------------------
!
!
DO JP = 1,IO%NPATCH
!
KK => NK%AL(JP)
PK => NP%AL(JP)
ISSK => NISS%AL(JP)
GK => NIG%AL(JP)
!
ALLOCATE(KK%XMPOTSAT(PK%NSIZE_P,IO%NGROUND_LAYER))
ALLOCATE(KK%XBCOEF (PK%NSIZE_P,IO%NGROUND_LAYER))
! needed to be written as diagnostics, so not free
ALLOCATE(KK%XWWILT (PK%NSIZE_P,IO%NGROUND_LAYER))
ALLOCATE(KK%XWFC (PK%NSIZE_P,IO%NGROUND_LAYER))
ALLOCATE(KK%XWSAT (PK%NSIZE_P,IO%NGROUND_LAYER))
!
CALL PACK_SAME_RANK(PK%NR_P,K%XMPOTSAT,KK%XMPOTSAT)
CALL PACK_SAME_RANK(PK%NR_P,K%XBCOEF,KK%XBCOEF)
!
CALL PACK_SAME_RANK(PK%NR_P,K%XWWILT,KK%XWWILT)
CALL PACK_SAME_RANK(PK%NR_P,K%XWFC,KK%XWFC)
CALL PACK_SAME_RANK(PK%NR_P,K%XWSAT,KK%XWSAT)
!