init_surf_atmn.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 init_surf_atm_n (YSC, &
                            hprogram,hinit, oland_use,                   &
                             ki,ksv,ksw,                                &
                             hsv,pco2,prhoa,                            &
                             pzenith,pazim,psw_bands,pdir_alb,psca_alb, &
                             pemis,ptsrad,ptsurf,                       &
                             kyear, kmonth,kday, ptime,                 &
                             hatmfile,hatmfiletype,                     &
                             htest                                      )  
!#############################################################
!
!!****  *INIT_SURF_ATM_n* - routine to initialize GROUND
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2003
!      (P.Tulet )             01/11/03  initialisation of the surface chemistry!
!!     (D.Gazen)    01/12/03  change emissions handling for surf. externalization
!!     (P.LeMoigne) 18/07/05  get 1d mask only if associated tile exists    
!!     (B.Decharme)  03/2009  New keys read for arrange cover by user
!!     (B.Decharme)  04/2009  Read precipitation forcing from the restart file for ARPEGE/ALADIN run
!!     (A. Lemonsu)    2009   New key read for urban green areas
!!     (B.Decharme)  07/2011  Read pgd+prep
!!     (S. Queguiner)  2011   Modif chemistry (2.4)
!!     (B. Decharme)   2013   Read grid only once in AROME case
!!     (G. Tanguy)     2013   Add IF(ALLOCATED(NMASK_FULL))  before deallocate
!!      B. Decharme  04/2013  new coupling variables
!!                            Delete LPROVAR_TO_DIAG check
!!                            Delete NWG_LAYER_TOT
!!     (J.Escobar)      10/06/2013: replace DOUBLE PRECISION by REAL to handle problem for promotion of real on IBM SP
!!     (J.Durand)      2014   add activation of chemical deposition if LCH_EMIS=F
!!      R. Séférian 03/2014   Adding decoupling between CO2 seen by photosynthesis and radiative CO2
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
!
USE modd_surfex_n, ONLY : surfex_t
!
USE modd_surf_atm,       ONLY : xco2uncpl
!
USE modd_read_namelist,  ONLY : lnam_read
USE modd_surf_conf,      ONLY : cprogname
USE modd_dst_surf,       ONLY : ndstmde, ndst_mdebeg, lvarsig_dst, lrgfix_dst 
USE modd_slt_surf,       ONLY : nsltmde, nslt_mdebeg, lvarsig_slt, lrgfix_slt                                

USE modd_data_cover_par, ONLY : ntilesfc
USE modd_data_cover,     ONLY : lclim_lai, xdata_lai_all_years, xdata_lai, &
                                neco2_start_year, neco2_end_year  
!
USE modd_surf_par,       ONLY : xundef, nundef
USE modd_chs_aerosol,    ONLY : lvarsigi, lvarsigj
USE modd_write_surf_atm, ONLY : lnowrite_canopy, lnowrite_texfile  
!
USE modd_surfex_mpi, ONLY : xtime_init_sea, xtime_init_water, xtime_init_nature, xtime_init_town, &
                            nrank, npio, nsize
USE modd_surfex_omp, ONLY : nindx2sfx, nwork, nwork2, xwork, xwork2, xwork3, &
                            nwork_full, nwork2_full, xwork_full, xwork2_full, &
                            nblocktot
!
USE modd_mask, ONLY: nmask_full
!
USE modi_init_io_surf_n
USE modi_default_sso
USE modi_default_ch_surf_atm
USE modi_default_diag_surf_atm
USE modi_read_default_surf_atm_n
USE modi_read_surf_atm_conf_n
USE modi_read_surf_atm_date
USE modi_read_nam_prep_surf_n
USE modi_read_surf
USE modi_get_sizes_parallel
USE modi_sunpos
USE modi_get_size_full_n
USE modi_read_cover_n
USE modi_read_sso_n
USE modi_subscale_z0eff
USE modi_read_sso_canopy_n
USE modi_read_dummy_n
USE modi_read_grid
USE modi_read_gridtype
USE modi_end_io_surf_n
USE modi_prep_ctrl_surf_atm
USE modi_average_rad
USE modi_average_tsurf
USE modi_init_chemical_n
USE modi_ch_init_depconst
USE modi_ch_init_emission_n
USE modi_ch_init_snap_n
USE modi_open_namelist
USE modi_close_namelist
USE modi_abor1_sfx
USE modi_alloc_diag_surf_atm_n
USE modi_get_1d_mask
USE modi_ini_data_cover
USE modi_init_inland_water_n
USE modi_init_nature_n
USE modi_init_sea_n
USE modi_init_town_n
USE modi_read_arrange_cover
USE modi_read_cover_garden
USE modi_read_eco2_irrig
USE modi_read_lclim_lai
USE modi_read_lecoclimap
USE modi_surf_version
USE modi_get_luout
USE modi_set_surfex_filein
!
USE modi_init_cpl_gcm_n
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
#ifdef AIX64
!$ USE OMP_LIB
#endif
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include 'mpif.h'
#endif
!
#ifndef AIX64
!$ INCLUDE 'omp_lib.h'
#endif
!
!*       0.1   Declarations of arguments
!              -------------------------
!
!
TYPE(surfex_t), INTENT(INOUT) :: ysc
!
 CHARACTER(LEN=6),                 INTENT(IN)  :: hprogram  ! program calling surf. schemes
 CHARACTER(LEN=3),                 INTENT(IN)  :: hinit     ! choice of fields to initialize
LOGICAL,                          INTENT(IN)  :: oland_use ! 
INTEGER,                          INTENT(IN)  :: ki        ! number of points
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(KI),  INTENT(IN)  :: pazim     ! solar azimuthal angle (rad from N, clock)
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)
!
INTEGER,                          INTENT(IN)  :: kyear     ! current year (UTC)
INTEGER,                          INTENT(IN)  :: kmonth    ! current month (UTC)
INTEGER,                          INTENT(IN)  :: kday      ! current day (UTC)
REAL,                             INTENT(IN)  :: ptime     ! current time since
                                                          !  midnight (UTC, s)
!
 CHARACTER(LEN=28),                INTENT(IN)  :: hatmfile    ! atmospheric file name
 CHARACTER(LEN=6),                 INTENT(IN)  :: hatmfiletype! atmospheric file type
 CHARACTER(LEN=2),                 INTENT(IN)  :: htest       ! must be equal to 'OK'
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
 CHARACTER(LEN=3)  :: yread
!
INTEGER           :: iswb     ! number of shortwave bands
INTEGER           :: jtile    ! loop counter on tiles
INTEGER           :: iresp    ! error return code
INTEGER           :: iluout   ! unit of output listing file
INTEGER           :: ich      ! unit of input chemical file
INTEGER           :: iversion, ibugfix       ! surface version
!
INTEGER, DIMENSION(:), ALLOCATABLE :: isize_omp
!
LOGICAL           :: lzenith  ! is the PZENITH field initialized ?
!
REAL, DIMENSION(:,:), ALLOCATABLE   :: zfrac_tile     ! fraction of each surface type
REAL, DIMENSION(KI,KSW,NTILESFC)    :: zdir_alb_tile  ! direct albedo
REAL, DIMENSION(KI,KSW,NTILESFC)    :: zsca_alb_tile  ! diffuse albedo
REAL, DIMENSION(KI,NTILESFC)        :: zemis_tile     ! emissivity
REAL, DIMENSION(KI,NTILESFC)        :: ztsrad_tile    ! radiative temperature
REAL, DIMENSION(KI,NTILESFC)        :: ztsurf_tile    ! effective temperature
REAL, DIMENSION(KI)                 :: zzenith        ! zenith angle
REAL, DIMENSION(KI)                 :: zazim          ! azimuth angle
REAL, DIMENSION(KI)                 :: ztsun          ! solar time since midnight
!
REAL, DIMENSION(:),     ALLOCATABLE :: zp_zenith   ! zenithal angle
REAL, DIMENSION(:),     ALLOCATABLE :: zp_azim     ! azimuthal angle
REAL, DIMENSION(:),     ALLOCATABLE :: zp_co2      ! air CO2 concentration
REAL, DIMENSION(:),     ALLOCATABLE :: zp_rhoa     ! air density
REAL, DIMENSION(:,:),   ALLOCATABLE :: zp_dir_alb  ! direct albedo
REAL, DIMENSION(:,:),   ALLOCATABLE :: zp_sca_alb  ! diffuse albedo
REAL, DIMENSION(:),     ALLOCATABLE :: zp_emis     ! emissivity
REAL, DIMENSION(:),     ALLOCATABLE :: zp_tsrad    ! radiative temperature
REAL, DIMENSION(:),     ALLOCATABLE :: zp_tsurf    ! surface effective temperature
!
REAL :: xtime0
!
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('INIT_SURF_ATM_N',0,zhook_handle)
!
!
 cprogname=hprogram
!
IF (htest/='OK') THEN
   CALL abor1_sfx('INIT_SURF_ATMN: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!
!-------------------------------------------------------------------------------
!
 CALL surf_version
!
!-------------------------------------------------------------------------------
!
 CALL get_luout(hprogram,iluout)
!
IF (lnam_read) THEN
 !
 !*       0.     Defaults
 !               --------
 !
 !        0.1. Hard defaults
 !      
 CALL default_sso(ysc%USS%CROUGH,ysc%USS%XFRACZ0,ysc%USS%XCOEFBE)
 CALL default_ch_surf_atm(ysc%CHU%CCHEM_SURF_FILE,ysc%CHU%LCH_SURF_EMIS)
 CALL default_diag_surf_atm(ysc%DGU%N2M,ysc%DGU%LT2MMW,ysc%DGU%LSURF_BUDGET,ysc%DGU%L2M_MIN_ZS,&
                            ysc%DGU%LRAD_BUDGET, &
                            ysc%DGU%LCOEF,ysc%DGU%LSURF_VARS,ysc%DGU%LSURF_BUDGETC,          &
                            ysc%DGU%LRESET_BUDGETC,ysc%DGU%LSELECT, ysc%DGU%LPROVAR_TO_DIAG, &
                            ysc%DGU%LDIAG_GRID, ysc%DGU%LFRAC, ysc%DGU%XDIAG_TSTEP )                       
 !
ENDIF
!
!        0.2. Defaults from file header
!    
 CALL read_default_surf_atm_n(ysc%CHU, ysc%DGU, ysc%USS, &
                              hprogram)
!
!*       1.     Reading of configuration
!               ------------------------
!
!        1.1. general options (diagnostics, etc...)
!
 CALL read_surf_atm_conf_n(ysc%CHU, ysc%DGU, ysc%USS, &
                           hprogram)
!
IF(xco2uncpl/=xundef)THEN
  WRITE(iluout,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
  WRITE(iluout,*)'!!!                                           !!!'
  WRITE(iluout,*)'!!!          WARNING    WARNING               !!!'
  WRITE(iluout,*)'!!!                                           !!!'
  WRITE(iluout,*)'!!! Decoupling between CO2 for photosynthesis !!!' 
  WRITE(iluout,*)'!!! and atmospheric CO2 activated             !!!'
  WRITE(iluout,*)'!!! In NAM_SURF_ATM XCO2UNCPL =',xco2uncpl,'  !!!'
  WRITE(iluout,*)'!!!                                           !!!'
  WRITE(iluout,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'        
ENDIF
!
!        1.2. Date
!
SELECT CASE (hinit)
  CASE ('PGD')
    ysc%U%TTIME%TDATE%YEAR = nundef
    ysc%U%TTIME%TDATE%MONTH= nundef
    ysc%U%TTIME%TDATE%DAY  = nundef
    ysc%U%TTIME%TIME       = xundef
        
  CASE ('PRE')
    ! check that diagnostics are off if hinit=='pre'
    CALL prep_ctrl_surf_atm(ysc%DGU%N2M,ysc%DGU%LSURF_BUDGET,ysc%DGU%L2M_MIN_ZS,ysc%DGU%LRAD_BUDGET,&
                            ysc%DGU%LCOEF,ysc%DGU%LSURF_VARS,    &
                            ysc%DGU%LSURF_BUDGETC,ysc%DGU%LRESET_BUDGETC,lnowrite_texfile,ysc%DGU%LSELECT,iluout,&
                            ysc%DGU%LPROVAR_TO_DIAG)  
    ! preparation of fields  (date not present in PGD file)
    IF (lnam_read) CALL read_nam_prep_surf_n(hprogram)
    CALL read_surf_atm_date(&
                            hprogram,hinit,iluout,hatmfile,hatmfiletype,kyear,kmonth,kday,ptime,ysc%U%TTIME)

  CASE default
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
    CALL read_surf(&
                   hprogram,'DTCUR',ysc%U%TTIME,iresp)
    CALL end_io_surf_n(hprogram)

END SELECT
!
!-----------------------------------------------------------------------------------------------------
! READ PGD FILE
!-----------------------------------------------------------------------------------------------------
!
!        1.3. Schemes used
!
!         Initialisation for IO
!
 CALL set_surfex_filein(hprogram,'PGD ') ! change input file name to pgd name
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
!
 CALL read_surf(&
                   hprogram,'VERSION',iversion,iresp)
 CALL read_surf(&
                   hprogram,'BUG',ibugfix,iresp)
!
IF (iversion>7 .OR. iversion==7 .AND.ibugfix>=2) THEN
  CALL read_surf(&
                   hprogram,'STORAGETYPE',yread,iresp)
ENDIF
!         reading
!
 CALL read_surf(&
                   hprogram,'SEA   ',ysc%U%CSEA   ,iresp)
 CALL read_surf(&
                   hprogram,'WATER ',ysc%U%CWATER ,iresp)
 CALL read_surf(&
                   hprogram,'NATURE',ysc%U%CNATURE,iresp)
 CALL read_surf(&
                   hprogram,'TOWN  ',ysc%U%CTOWN  ,iresp)
!
 CALL read_surf(&
                   hprogram,'DIM_FULL  ',ysc%U%NDIM_FULL,  iresp)
IF (hinit=='PRE') THEN
  nindx2sfx = ysc%U%NDIM_FULL
  nsize = ysc%U%NDIM_FULL
  CALL end_io_surf_n(hprogram)
  !Initialize full mask with good dimension
  IF (ALLOCATED(nmask_full)) DEALLOCATE(nmask_full)
  CALL set_surfex_filein(hprogram,'PGD ')
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
  ALLOCATE(nwork(ysc%U%NDIM_FULL))
  ALLOCATE(xwork(ysc%U%NDIM_FULL))
  ALLOCATE(nwork2(ysc%U%NDIM_FULL,10))  
  ALLOCATE(xwork2(ysc%U%NDIM_FULL,10))
  ALLOCATE(xwork3(ysc%U%NDIM_FULL,10,10))
  IF (nrank==npio) THEN
    ALLOCATE(nwork_full(ysc%U%NDIM_FULL))
    ALLOCATE(xwork_full(ysc%U%NDIM_FULL))
    ALLOCATE(nwork2_full(ysc%U%NDIM_FULL,10))
    ALLOCATE(xwork2_full(ysc%U%NDIM_FULL,10))
  ELSE
    ALLOCATE(nwork_full(0))
    ALLOCATE(xwork_full(0))
    ALLOCATE(nwork2_full(0,0))
    ALLOCATE(xwork2_full(0,0))
  ENDIF
ENDIF  
 CALL read_surf(&
                   hprogram,'DIM_SEA   ',ysc%U%NDIM_SEA,   iresp)
 CALL read_surf(&
                   hprogram,'DIM_NATURE',ysc%U%NDIM_NATURE,iresp)
 CALL read_surf(&
                   hprogram,'DIM_WATER ',ysc%U%NDIM_WATER, iresp)
 CALL read_surf(&
                   hprogram,'DIM_TOWN  ',ysc%U%NDIM_TOWN,  iresp)
 CALL read_lecoclimap(&
                      hprogram,ysc%U%LECOCLIMAP)
 CALL read_arrange_cover(&
                         hprogram,ysc%U%LWATER_TO_NATURE,ysc%U%LTOWN_TO_ROCK)
 CALL read_cover_garden(&
                        hprogram,ysc%U%LGARDEN)
!
!* reads if climatological LAI is used or not for ecoclimap2. If not, looks for year to be used.
 CALL read_lclim_lai(&
                     hprogram,lclim_lai)
IF (.NOT. lclim_lai .AND. ysc%U%TTIME%TDATE%YEAR >= neco2_start_year &
                     .AND. ysc%U%TTIME%TDATE%YEAR <= neco2_end_year   ) ysc%DTCO%NYEAR=ysc%U%TTIME%TDATE%YEAR
 CALL ini_data_cover(ysc%DTCO, ysc%U)
 CALL read_eco2_irrig(&
                      ysc%DTCO, &
                      hprogram)
!
!*       2.     Cover fields and grid:
!               ---------------------
!
!        2.0. Get number of points on this proc
!
 CALL get_size_full_n(ysc%U, &
                      hprogram,ysc%U%NDIM_FULL,ysc%U%NSIZE_FULL)
!
!        2.1. Read cover
!
 CALL read_cover_n(ysc%DTCO, ysc%U, &
                   hprogram)
!
!        2.2. Read grid
!
ALLOCATE(ysc%UG%XLAT       (ysc%U%NSIZE_FULL))
ALLOCATE(ysc%UG%XLON       (ysc%U%NSIZE_FULL))
ALLOCATE(ysc%UG%XMESH_SIZE (ysc%U%NSIZE_FULL))
ALLOCATE(ysc%USS%XZ0EFFJPDIR(ysc%U%NSIZE_FULL))
 CALL read_grid(&
                hprogram,ysc%UG%CGRID,ysc%UG%XGRID_PAR,ysc%UG%XLAT,ysc%UG%XLON,ysc%UG%XMESH_SIZE,iresp,ysc%USS%XZ0EFFJPDIR)
ysc%UG%NGRID_PAR=SIZE(ysc%UG%XGRID_PAR)
!
!        2.3. Initialize zenith and azimuth angles if not done yet
!
lzenith = all(pzenith /= xundef)
IF (.NOT. lzenith) THEN
!$OMP PARALLEL
!$ NBLOCKTOT = OMP_GET_NUM_THREADS()
!$OMP END PARALLEL
  ALLOCATE(isize_omp(0:nblocktot-1))
  CALL get_sizes_parallel(ysc%DTCO, ysc%DGU, ysc%UG, ysc%U, &
                          nblocktot,ki,0,isize_omp)
  CALL sunpos(isize_omp,kyear, kmonth, kday, ptime, ysc%UG%XLON, ysc%UG%XLAT, ztsun, zzenith, zazim)
  DEALLOCATE(isize_omp)
ENDIF

!
IF (hprogram/='AROME '.AND.nrank==npio) THEN
  !
  IF (.NOT.ASSOCIATED(ysc%UG%XGRID_FULL_PAR)) THEN
    CALL read_gridtype(&
                       hprogram,ysc%UG%CGRID,ysc%UG%NGRID_PAR,ysc%U%NDIM_FULL,.false.,hdir='A')
    ALLOCATE(ysc%UG%XGRID_FULL_PAR(ysc%UG%NGRID_PAR))
    CALL read_gridtype(&
                       hprogram,ysc%UG%CGRID,ysc%UG%NGRID_PAR,ysc%U%NDIM_FULL,.true.,&
                       ysc%UG%XGRID_FULL_PAR,iresp,hdir='A')
  ENDIF
  !
ENDIF
!
!*       2.4     Allocation of chemical species name, chemical index of HSV array 
!
 CALL init_chemical_n(iluout, ksv, hsv, ysc%SV,        &
                     ysc%CHU%CCH_NAMES, ysc%CHU%CAER_NAMES     )
!
!        2.4 Initialize Chemical Emissions
!
 CALL read_surf(&
                   hprogram,'CH_EMIS',ysc%CHU%LCH_EMIS,iresp)
!
IF (ysc%CHU%LCH_EMIS) THEN
  !
  IF ( iversion<7 .OR. iversion==7 .AND. ibugfix<3 ) THEN
    ysc%CHU%CCH_EMIS='AGGR'
  ELSE
    CALL read_surf(&
                   hprogram,'CH_EMIS_OPT',ysc%CHU%CCH_EMIS,iresp)
  END IF
  !
  IF (ysc%SV%NBEQ > 0) THEN
    !
    CALL open_namelist(hprogram,ich,hfile=ysc%CHU%CCHEM_SURF_FILE)
    !
    IF (ysc%CHU%LCH_SURF_EMIS) THEN
      IF (ysc%CHU%CCH_EMIS=='AGGR') THEN
        CALL ch_init_emission_n(&
                                ysc%CHE, ysc%CHU, ysc%SV, &
                                hprogram,ysc%U%NSIZE_FULL,ich,prhoa) 
      ELSE
        CALL ch_init_snap_n(&
                            ysc%CHN, ysc%SV, &
                            hprogram,ysc%U%NSIZE_FULL,hinit,ich,prhoa)
      END IF
    ENDIF
    CALL close_namelist(hprogram,ich)
  ENDIF
  !
END IF
    !
    !*       2.5 Initialization of dry deposition scheme (chemistry)
    !    
!
IF (ysc%SV%NBEQ .GT. 0) THEN
  CALL open_namelist(hprogram,ich,hfile=ysc%CHU%CCHEM_SURF_FILE)
!
  IF (hinit=='ALL') CALL ch_init_depconst(ich,iluout,ysc%SV%CSV(ysc%SV%NSV_CHSBEG:ysc%SV%NSV_CHSEND))
!
  CALL close_namelist(hprogram,ich)
!
END IF
!
!*       2.5 Subgrid orography
!
 CALL read_sso_n(&
                 ysc%U, ysc%USS, &
                 hprogram)
!
!*       2.6 Orographic roughness length
!
ALLOCATE(ysc%USS%XZ0EFFIP(ysc%U%NSIZE_FULL))
ALLOCATE(ysc%USS%XZ0EFFIM(ysc%U%NSIZE_FULL))
ALLOCATE(ysc%USS%XZ0EFFJP(ysc%U%NSIZE_FULL))
ALLOCATE(ysc%USS%XZ0EFFJM(ysc%U%NSIZE_FULL))
ALLOCATE(ysc%USS%XZ0REL  (ysc%U%NSIZE_FULL))
!
 CALL subscale_z0eff(ysc%USS%XAOSIP,ysc%USS%XAOSIM,ysc%USS%XAOSJP,ysc%USS%XAOSJM,         &
                    ysc%USS%XHO2IP,ysc%USS%XHO2IM,ysc%USS%XHO2JP,ysc%USS%XHO2JM,0.,      &
                    ysc%USS%XZ0EFFIP,ysc%USS%XZ0EFFIM,ysc%USS%XZ0EFFJP,ysc%USS%XZ0EFFJM, &
                    ysc%USS%XZ0REL                               )
!
!*       2.7 Dummy fields
!
 CALL read_dummy_n(&
                   ysc%DUU,ysc% U, &
                   hprogram)
!
!         End of IO
!
 CALL end_io_surf_n(hprogram)
 CALL set_surfex_filein(hprogram,'PREP') ! restore input file name
!
!-----------------------------------------------------------------------------------------------------
! END READ PGD FILE
!-----------------------------------------------------------------------------------------------------
!
!
!         Initialisation for IO
!
 CALL init_io_surf_n(ysc%DTCO, ysc%DGU, ysc%U, &
                        hprogram,'FULL  ','SURF  ','READ ')
!
!*       2.8 Allocations and Initialization of diagnostics
!
IF (hinit=='ALL') CALL alloc_diag_surf_atm_n(&
                                             ysc%DGU, ysc%U, &
                                             hprogram,ksw)
!
!
!*       Canopy fields if Beljaars et al 2004 parameterization is used
!
IF (ysc%USS%CROUGH=='BE04') CALL read_sso_canopy_n(ysc%DTCO, ysc%SSCP, ysc%U, &
                                               hprogram,hinit)
!
!*       Physical fields need for ARPEGE/ALADIN climate run
!
 CALL init_cpl_gcm_n(&
                     ysc%U, &
                     hprogram,hinit)
!
!         End of IO
!
 CALL end_io_surf_n(hprogram)
!
!-----------------------------------------------------------------------------------------------------
!
!*       4.     Initialization of masks for each surface
!               ----------------------------------------
!
!* number of geographical points
ysc%U%NSIZE_NATURE    = count(ysc%U%XNATURE(:) > 0.0)
ysc%U%NSIZE_TOWN      = count(ysc%U%XTOWN(:)   > 0.0)
ysc%U%NSIZE_WATER     = count(ysc%U%XWATER(:)  > 0.0)
ysc%U%NSIZE_SEA       = count(ysc%U%XSEA(:)    > 0.0)
!
ALLOCATE(ysc%U%NR_NATURE (ysc%U%NSIZE_NATURE))
ALLOCATE(ysc%U%NR_TOWN   (ysc%U%NSIZE_TOWN  ))
ALLOCATE(ysc%U%NR_WATER  (ysc%U%NSIZE_WATER ))
ALLOCATE(ysc%U%NR_SEA    (ysc%U%NSIZE_SEA   ))
!
IF (ysc%U%NSIZE_SEA   >0)CALL get_1d_mask( ysc%U%NSIZE_SEA,    ysc%U%NSIZE_FULL, ysc%U%XSEA   , ysc%U%NR_SEA   )
IF (ysc%U%NSIZE_WATER >0)CALL get_1d_mask( ysc%U%NSIZE_WATER,  ysc%U%NSIZE_FULL, ysc%U%XWATER , ysc%U%NR_WATER )
IF (ysc%U%NSIZE_TOWN  >0)CALL get_1d_mask( ysc%U%NSIZE_TOWN,   ysc%U%NSIZE_FULL, ysc%U%XTOWN  , ysc%U%NR_TOWN  )
IF (ysc%U%NSIZE_NATURE>0)CALL get_1d_mask( ysc%U%NSIZE_NATURE, ysc%U%NSIZE_FULL, ysc%U%XNATURE, ysc%U%NR_NATURE)
!
!* number of shortwave spectral bands
iswb=SIZE(psw_bands)
!
!* tile number
ALLOCATE(zfrac_tile(ysc%U%NSIZE_FULL,ntilesfc))
jtile = 0
!
!
!*       5.     Default values
!               --------------
!
zdir_alb_tile = xundef
zsca_alb_tile = xundef
zemis_tile    = xundef
ztsrad_tile   = xundef
ztsurf_tile   = xundef
!
#ifdef SFX_MPI
xtime0 = mpi_wtime()
#endif
!
!*       6.     Initialization of sea
!               ---------------------
!
jtile               = jtile + 1
zfrac_tile(:,jtile) = ysc%U%XSEA(:)
!
! pack variables which are arguments to this routine
 CALL pack_surf_init_arg(ysc%U%NSIZE_SEA,ysc%U%NR_SEA)
!
! initialization
IF (ysc%U%NDIM_SEA>0) &
  CALL init_sea_n(ysc%DTCO, ysc%DGU, ysc%UG, ysc%U, &
                  ysc%SM, ysc%DGL,  &
                  hprogram,hinit,ysc%U%NSIZE_SEA,ksv,ksw,            &
                  hsv,zp_co2,zp_rhoa,                                &
                  zp_zenith,zp_azim,psw_bands,zp_dir_alb,zp_sca_alb, &
                  zp_emis,zp_tsrad,zp_tsurf,                         &
                  kyear,kmonth,kday,ptime, hatmfile,hatmfiletype,    &
                  'OK'                                               )  
!
!
 CALL unpack_surf_init_arg(jtile,ysc%U%NSIZE_SEA,ysc%U%NR_SEA)  
!
#ifdef SFX_MPI
xtime_init_sea = xtime_init_sea + (mpi_wtime() - xtime0)*100./max(1,ysc%U%NSIZE_SEA)
xtime0 = mpi_wtime()
#endif
!
!*       7.     Initialization of lakes
!               -----------------------
!
!
jtile               = jtile + 1
zfrac_tile(:,jtile) = ysc%U%XWATER(:)
!
! pack variables which are arguments to this routine
 CALL pack_surf_init_arg(ysc%U%NSIZE_WATER,ysc%U%NR_WATER)
!
! initialization
IF (ysc%U%NDIM_WATER>0) &
  CALL init_inland_water_n(ysc%DTCO, ysc%DGU,ysc%UG, &
                           ysc%U, ysc%WM, ysc%FM, ysc%DGL,    &
                           hprogram,hinit,ysc%U%NSIZE_WATER,ksv,ksw,          &
                           hsv,zp_co2,zp_rhoa,                                &
                           zp_zenith,zp_azim,psw_bands,zp_dir_alb,zp_sca_alb, &
                           zp_emis,zp_tsrad,zp_tsurf,                         &
                           kyear,kmonth,kday,ptime, hatmfile,hatmfiletype,    &
                           'OK'                                               )
!
 CALL unpack_surf_init_arg(jtile,ysc%U%NSIZE_WATER,ysc%U%NR_WATER)
!
#ifdef SFX_MPI
xtime_init_water = xtime_init_water + (mpi_wtime() - xtime0)*100./max(1,ysc%U%NSIZE_WATER)
xtime0 = mpi_wtime()
#endif
!
!*       8.     Initialization of vegetation scheme
!               -----------------------------------
!
!
jtile               = jtile + 1
zfrac_tile(:,jtile) = ysc%U%XNATURE(:)
!
! pack variables which are arguments to this routine
 CALL pack_surf_init_arg(ysc%U%NSIZE_NATURE,ysc%U%NR_NATURE)
!
! initialization
IF (ysc%U%NDIM_NATURE>0) &
  CALL init_nature_n(ysc%DTCO, ysc%DGU, ysc%UG, ysc%U, ysc%IM, &
                     ysc%DTZ, ysc%DGL, ysc%DST, ysc%SLT, ysc%SV, &
                     hprogram,hinit,oland_use,ysc%U%NSIZE_NATURE,ksv,ksw,   &
                     hsv,zp_co2,zp_rhoa,                                &
                     zp_zenith,zp_azim,psw_bands,zp_dir_alb,zp_sca_alb, &
                     zp_emis,zp_tsrad,zp_tsurf,                         &
                     kyear,kmonth,kday,ptime, hatmfile,hatmfiletype,    &
                     'OK'                                               )
!
!
 CALL unpack_surf_init_arg(jtile,ysc%U%NSIZE_NATURE,ysc%U%NR_NATURE)  
!
#ifdef SFX_MPI
xtime_init_nature = xtime_init_nature + (mpi_wtime() - xtime0)*100./max(1,ysc%U%NSIZE_NATURE)
xtime0 = mpi_wtime()
#endif
!
!*       9.     Initialization of urban scheme
!               ------------------------------
!
jtile               = jtile + 1
zfrac_tile(:,jtile) = ysc%U%XTOWN(:)
!
! pack variables which are arguments to this routine
 CALL pack_surf_init_arg(ysc%U%NSIZE_TOWN,ysc%U%NR_TOWN)
!
! initialization
IF (ysc%U%NDIM_TOWN>0) &
  CALL init_town_n(ysc%DTCO, ysc%DGU, ysc%UG, ysc%U, &
                   ysc%IM%CHI, ysc%IM%DTI, ysc%IM%I, &
                   ysc%TM, ysc%GDM, ysc%GRM, ysc%DGL, ysc%DST, ysc%SLT, &
                   hprogram,hinit,ysc%U%NSIZE_TOWN,ksv,ksw,             &
                   hsv,zp_co2,zp_rhoa,                                &
                   zp_zenith,zp_azim,psw_bands,zp_dir_alb,zp_sca_alb, &
                   zp_emis,zp_tsrad,zp_tsurf,                         &
                   kyear,kmonth,kday,ptime, hatmfile,hatmfiletype,    &
                   'OK'                                               )  
!
!
 CALL unpack_surf_init_arg(jtile,ysc%U%NSIZE_TOWN,ysc%U%NR_TOWN)  
!
#ifdef SFX_MPI
xtime_init_town = xtime_init_town + (mpi_wtime() - xtime0)*100./max(1,ysc%U%NSIZE_TOWN)
#endif
!
!
!*      10.     Output radiative and physical fields
!               ------------------------------------
!
IF (SIZE(pdir_alb)>0)                                                   &
  CALL average_rad(zfrac_tile,                                            &
                   zdir_alb_tile, zsca_alb_tile, zemis_tile, ztsrad_tile, &
                   pdir_alb,      psca_alb,      pemis,      ptsrad       ) 
!
IF (SIZE(ptsurf)>0) &
  CALL average_tsurf(zfrac_tile, ztsurf_tile, ptsurf)
!                  
DEALLOCATE(zfrac_tile)
!
!-------------------------------------------------------------------------------
!==============================================================================
IF (lhook) CALL dr_hook('INIT_SURF_ATM_N',1,zhook_handle)
 CONTAINS
!==============================================================================
SUBROUTINE pack_surf_init_arg(KSIZE,KMASK)
!
INTEGER, INTENT(IN)               :: ksize
INTEGER, INTENT(IN), DIMENSION(:) :: kmask
INTEGER :: jj
REAL(KIND=JPRB) :: zhook_handle
!
! input arguments:
!
IF (lhook) CALL dr_hook('PACK_SURF_INIT_ARG',0,zhook_handle)
ALLOCATE(zp_co2(ksize))
ALLOCATE(zp_rhoa(ksize))
ALLOCATE(zp_zenith(ksize))
ALLOCATE(zp_azim(ksize))
!
!
! output arguments:
!
ALLOCATE(zp_dir_alb(ksize,iswb))
ALLOCATE(zp_sca_alb(ksize,iswb))
ALLOCATE(zp_emis(ksize))
ALLOCATE(zp_tsrad(ksize))
ALLOCATE(zp_tsurf(ksize))
!
IF (ksize>0) THEN
  zp_co2    = 6.e-4
  zp_rhoa   = 1.2
  zp_zenith = 0.
  zp_azim   = 0.
  zp_dir_alb = xundef
  zp_sca_alb = xundef
  zp_emis    = xundef
  zp_tsrad   = xundef
  zp_tsurf   = xundef
END IF
!
DO jj=1,ksize
IF (SIZE(pco2)>0) &
     zp_co2(jj)     = pco2(kmask(jj))  
IF (SIZE(prhoa)>0) &
     zp_rhoa(jj)     = prhoa(kmask(jj))  
IF (SIZE(pzenith)>0) THEN
    IF (lzenith) THEN
       zp_zenith(jj)     = pzenith(kmask(jj)) 
    ELSE
       zp_zenith(jj)     = zzenith(kmask(jj)) 
    ENDIF
ENDIF
IF (SIZE(pazim  )>0) THEN
    IF (lzenith) THEN
       zp_azim(jj)     = pazim(kmask(jj)) 
    ELSE
       zp_azim(jj)     = zazim(kmask(jj)) 
    ENDIF
ENDIF
ENDDO
IF (lhook) CALL dr_hook('PACK_SURF_INIT_ARG',1,zhook_handle)
!
END SUBROUTINE pack_surf_init_arg
!==============================================================================
SUBROUTINE unpack_surf_init_arg(KTILE,KSIZE,KMASK)
!
INTEGER, INTENT(IN) :: ktile, ksize
!
INTEGER, INTENT(IN), DIMENSION(:) :: kmask
!
INTEGER :: jj   ! loop counter
REAL(KIND=JPRB) :: zhook_handle
!
!
IF (lhook) CALL dr_hook('UNPACK_SURF_INIT_ARG',0,zhook_handle)
DO jj=1,ksize
IF (SIZE(ztsrad_tile)>0) &
     ztsrad_tile(kmask(jj),ktile)  = zp_tsrad(jj)  
IF (SIZE(zdir_alb_tile)>0) &
     zdir_alb_tile(kmask(jj),:,ktile)= zp_dir_alb(jj,:)  
IF (SIZE(zsca_alb_tile)>0) &
     zsca_alb_tile(kmask(jj),:,ktile)= zp_sca_alb(jj,:)  
IF (SIZE(zemis_tile)>0) &
     zemis_tile(kmask(jj),ktile)  = zp_emis(jj)
IF (SIZE(ztsurf_tile)>0) &
     ztsurf_tile(kmask(jj),ktile)  = zp_tsurf(jj)
ENDDO
!
DEALLOCATE(zp_co2    )
DEALLOCATE(zp_rhoa   )
DEALLOCATE(zp_zenith )
DEALLOCATE(zp_azim   )
DEALLOCATE(zp_dir_alb)
DEALLOCATE(zp_sca_alb)
DEALLOCATE(zp_emis   )
DEALLOCATE(zp_tsrad  )
DEALLOCATE(zp_tsurf  )
IF (lhook) CALL dr_hook('UNPACK_SURF_INIT_ARG',1,zhook_handle)
!
END SUBROUTINE unpack_surf_init_arg
!==============================================================================
!
END SUBROUTINE init_surf_atm_n