pgd_orography.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 pgd_orography (DTCO, UG, U, USS, HPROGRAM, HFILE, HFILETYPE, OZS)
!     ##############################################################
!
!!**** *PGD_OROGRAPHY* monitor for averaging and interpolations of cover fractions
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!    V. Masson        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    10/12/97
!!    12/2008 E. Martin : add case 'MAX' for choice of orography
!!    11/2012 M. Lafaysse : read ZS from a NETCDF file at the same resolution
!!    07/2013 M. Lafaysse : explicit slope from resolved orography
!!
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
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_sso_n, ONLY : sso_t
!
!
USE modd_surfex_mpi, ONLY : nrank, npio
!
USE modd_pgd_grid,       ONLY : nl, cgrid, xgrid_par
USE modd_pgdwork,        ONLY : xall, nsize_all, xssqo, lssqo, nsso, &
                                xext_all, xsumval, nsize
USE modd_surf_par,       ONLY : xundef, nundef
!
USE modi_get_luout
USE modi_open_aux_io_surf
USE modi_close_aux_io_surf
USE modi_read_nam_pgd_orography
USE modi_read_surf
USE modi_treat_field
USE modi_read_pgd_netcdf
USE modi_interpol_field
USE modi_sso
USE modi_subscale_aos
USE modi_get_size_full_n
USE modi_test_nam_var_surf
!
USE modi_read_sso_n
USE modi_init_io_surf_n
USE modi_end_io_surf_n
#ifdef SFX_ASC
USE modd_io_surf_asc, ONLY : cfilein
#endif
#ifdef SFX_FA
USE modd_io_surf_fa,  ONLY : cfilein_fa
#endif
#ifdef SFX_LFI
USE modd_io_surf_lfi, ONLY : cfilein_lfi
#endif
!
USE modi_explicit_slope

USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!

USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(data_cover_t), INTENT(INOUT) :: DTCO
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
TYPE(surf_atm_t), INTENT(INOUT) :: U
TYPE(sso_t), INTENT(INOUT) :: USS
!
 CHARACTER(LEN=6),     INTENT(IN)  :: HPROGRAM ! program calling
 CHARACTER(LEN=28),    INTENT(IN)  :: HFILE    ! atmospheric file name
 CHARACTER(LEN=6),     INTENT(IN)  :: HFILETYPE! atmospheric file type
LOGICAL,              INTENT(IN)  :: OZS      ! .true. if orography is imposed by atm. model
!
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER :: ILUOUT    ! output listing logical unit
!

REAL,DIMENSION(:),POINTER :: ZSLOPE ! degrees
INTEGER::JJ
REAL,PARAMETER :: PP_DEG2RAD= 3.141592654/180.
LOGICAL:: LPRESENT

LOGICAL, DIMENSION(NL)   :: GSSO        ! mask where SSO are computed
LOGICAL, DIMENSION(NL)   :: GSSO_ANIS   ! mask where SSO anisotropy is computed
LOGICAL, DIMENSION(NL)   :: GZ0EFFI     ! mask where z0  is  computed in subgrid direction x
LOGICAL, DIMENSION(NL)   :: GZ0EFFJ     ! mask where z0  is  computed in subgrid direction y
INTEGER, DIMENSION(NL)   :: IFLAG       ! flag for SSO and z0 fields interpolations
INTEGER                  :: IRESP       ! error code
REAL                     :: ZEPS = 1.e-4! a small number
INTEGER                  :: IDIM_FULL   ! total size of orographic array in atmospheric file
INTEGER                  :: IZS         ! size of orographic array in atmospheric file
!
!*    0.3    Declaration of namelists
!            ------------------------
!
 CHARACTER(LEN=28)        :: YZS         ! file name for orography
 CHARACTER(LEN=6)         :: YFILETYPE   ! data file type
 CHARACTER(LEN=28)        :: YSLOPE         ! file name for orography
 CHARACTER(LEN=6)         :: YSLOPEFILETYPE   ! data file type
REAL                     :: XUNIF_ZS    ! uniform orography
 CHARACTER(LEN=3)         :: COROGTYPE   ! orogpraphy type 
!                                       ! 'AVG' : average orography
!                                       ! 'SIL' : silhouette orography
!                                       ! 'ENV' : enveloppe orography
REAL                     :: XENV        ! parameter for enveloppe orography:
!                                       ! zs = avg_zs + XENV * SSO_STEDV
LOGICAL                  :: LIMP_ZS     ! Imposed orography from another PGD file
LOGICAL                  :: LEXPLICIT_SLOPE ! Slope is computed from explicit ZS field and not subgrid orography

REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
!-------------------------------------------------------------------------------
!
!*    1.      Initializations
!             ---------------
!
IF (lhook) CALL dr_hook('PGD_OROGRAPHY',0,zhook_handle)
 CALL get_luout(hprogram,iluout)
!
!-------------------------------------------------------------------------------
!
!*    2.      Reading of namelist
!             -------------------
!
 CALL read_nam_pgd_orography(hprogram, yzs, yfiletype, xunif_zs, &
                              corogtype, xenv, limp_zs, &
                              yslope, yslopefiletype, lexplicit_slope)  
!
 CALL test_nam_var_surf(iluout,'YSLOPEFILETYPE',yslopefiletype,'      ','NETCDF')
!
!-------------------------------------------------------------------------------
!
!*    3.      Allocations of orographic arrays
!             --------------------------------
!
ALLOCATE(u%XZS        (nl))
!
ALLOCATE(uss%XAVG_ZS    (nl))
ALLOCATE(uss%XSIL_ZS    (nl))
ALLOCATE(uss%XSSO_STDEV (nl))
ALLOCATE(uss%XMIN_ZS    (nl))
ALLOCATE(uss%XMAX_ZS    (nl))
!
ALLOCATE(uss%XSSO_ANIS  (nl))
ALLOCATE(uss%XSSO_DIR   (nl))
ALLOCATE(uss%XSSO_SLOPE (nl))
!
ALLOCATE(uss%XAOSIP     (nl))
ALLOCATE(uss%XAOSIM     (nl))
ALLOCATE(uss%XAOSJP     (nl))
ALLOCATE(uss%XAOSJM     (nl))
ALLOCATE(uss%XHO2IP     (nl))
ALLOCATE(uss%XHO2IM     (nl))
ALLOCATE(uss%XHO2JP     (nl))
ALLOCATE(uss%XHO2JM     (nl))
!
u%XZS       (:) = xundef
uss%XAVG_ZS   (:) = xundef
uss%XSIL_ZS   (:) = xundef
uss%XSSO_STDEV(:) = xundef
uss%XMIN_ZS   (:) = 99999.
uss%XMAX_ZS   (:) =-99999. 
!
uss%XSSO_ANIS (:) = xundef
uss%XSSO_DIR  (:) = xundef
uss%XSSO_SLOPE(:) = xundef
!
uss%XAOSIP    (:) = xundef
uss%XAOSIM    (:) = xundef
uss%XAOSJP    (:) = xundef
uss%XAOSJM    (:) = xundef
uss%XHO2IP    (:) = xundef
uss%XHO2IM    (:) = xundef
uss%XHO2JP    (:) = xundef
uss%XHO2JM    (:) = xundef
!-------------------------------------------------------------------------------
!
!*    4.      Allocations of work arrays
!             --------------------------
!
!-------------------------------------------------------------------------------
!
!*    5.      Uniform field is prescribed
!             ---------------------------
!
IF (ozs) THEN
!
!*    5.1     Use of imposed field
!             --------------------
!
  CALL open_aux_io_surf(hfile,hfiletype,'FULL  ')
  CALL read_surf(hfiletype,'DIM_FULL  ',idim_full,iresp)
  CALL get_size_full_n(hprogram,idim_full,u%NSIZE_FULL,izs)
  IF (izs /= nl) THEN
    WRITE(iluout,*) ' '
    WRITE(iluout,*) '***********************************************************'
    WRITE(iluout,*) '* Error in orography preparation                          *'
    WRITE(iluout,*) '* Prescribed orography from atmospheric model does not    *'
    WRITE(iluout,*) '* have the correct number of points                       *'
    WRITE(iluout,*) '* number of points in atmospheric orography: ', izs
    WRITE(iluout,*) '* number of points in the surface          : ', nl
    WRITE(iluout,*) '***********************************************************'
    WRITE(iluout,*) ' '
    CALL abor1_sfx('PGD_OROGRAPHY: ATMOSPHERIC PRESCRIBED OROGRAPHY DOES NOT HAVE THE CORRECT NB OF POINTS')
  END IF
  CALL read_surf(hfiletype,'ZS',u%XZS(:),iresp)
  CALL close_aux_io_surf(hfile,hfiletype)
  !
  uss%XAVG_ZS(:)    = u%XZS(:)
  uss%XSIL_ZS(:)    = u%XZS(:)
  uss%XMIN_ZS(:)    = u%XZS(:)
  uss%XMAX_ZS(:)    = u%XZS(:)
  uss%XSSO_STDEV(:) = 0.
  uss%XHO2IP(:)     = 0.
  uss%XHO2IM(:)     = 0.
  uss%XHO2JP(:)     = 0.
  uss%XHO2JM(:)     = 0.
  uss%XAOSIP(:)     = 0.
  uss%XAOSIM(:)     = 0.
  uss%XAOSJP(:)     = 0.
  uss%XAOSJM(:)     = 0.
  uss%XSSO_ANIS(:)  = 0.
  uss%XSSO_DIR(:)   = 0.
  uss%XSSO_SLOPE(:) = 0.

  IF (lhook) CALL dr_hook('PGD_OROGRAPHY',1,zhook_handle)
  RETURN

!
ELSE IF (xunif_zs/=xundef) THEN
!
!*    5.2     Use of the presribed cover fractions
!             ------------------------------------
!
  u%XZS(:)        = xunif_zs
  uss%XAVG_ZS(:)    = u%XZS(:)
  uss%XSIL_ZS(:)    = u%XZS(:)
  uss%XMIN_ZS(:)    = u%XZS(:)
  uss%XMAX_ZS(:)    = u%XZS(:)
  uss%XSSO_STDEV(:) = 0.
  uss%XHO2IP(:)     = 0.
  uss%XHO2IM(:)     = 0.
  uss%XHO2JP(:)     = 0.
  uss%XHO2JM(:)     = 0.
  uss%XAOSIP(:)     = 0.
  uss%XAOSIM(:)     = 0.
  uss%XAOSJP(:)     = 0.
  uss%XAOSJM(:)     = 0.
  uss%XSSO_ANIS(:)  = 0.
  uss%XSSO_DIR(:)   = 0.
  uss%XSSO_SLOPE(:) = 0.

  IF (lhook) CALL dr_hook('PGD_OROGRAPHY',1,zhook_handle)
  RETURN
!
!*    5.3     No data
!             -------
!
ELSEIF (len_trim(yzs)==0) THEN

  WRITE(iluout,*) ' '
  WRITE(iluout,*) '***********************************************************'
  WRITE(iluout,*) '* Error in orography preparation                          *'
  WRITE(iluout,*) '* There is no prescribed orography and no input file      *'
  WRITE(iluout,*) '***********************************************************'
  WRITE(iluout,*) ' '
  CALL abor1_sfx('PGD_OROGRAPHY: NO PRESCRIBED OROGRAPHY NOR INPUT FILE')
!  
ELSEIF (limp_zs) THEN !LIMP_ZS (impose topo from input file at the same resolution)
!
  IF(yfiletype=='NETCDF')THEN
     
!      CALL ABOR1_SFX('Use another format than netcdf for topo input file with LIMP_ZS')
    CALL read_pgd_netcdf(ug, u, uss, &
                          hprogram,'SURF  ','      ',yzs,'ZS                  ',u%XZS)
     
    uss%XSIL_ZS(:)    = u%XZS(:)
    uss%XAVG_ZS(:)    = u%XZS(:)
    uss%XMIN_ZS(:)    = u%XZS(:)
    uss%XMAX_ZS(:)    = u%XZS(:)
    uss%XSSO_STDEV(:) = 0.
    uss%XHO2IP(:)     = 0.
    uss%XHO2IM(:)     = 0.
    uss%XHO2JP(:)     = 0.
    uss%XHO2JM(:)     = 0.
    uss%XAOSIP(:)     = 0.
    uss%XAOSIM(:)     = 0.
    uss%XAOSJP(:)     = 0.
    uss%XAOSJM(:)     = 0.
    uss%XSSO_ANIS(:)  = 0.
    uss%XSSO_DIR(:)   = 0.
     
     
    ! read slope in file
    IF (len_trim(yslope)/=0) THEN
      ALLOCATE(zslope(nl))

    ! Read field on the same grid as FORCING
      CALL read_pgd_netcdf(ug, u, uss, &
                          hprogram,'SURF  ','      ',yslope,'slope               ',zslope)

      DO jj=1,nl
        uss%XSSO_SLOPE(jj)=tan(zslope(jj)*pp_deg2rad)
      END DO
      DEALLOCATE(zslope)     
    ELSE
      uss%XSSO_SLOPE=0.
    ENDIF
      
  ELSE
#ifdef SFX_ASC
    cfilein     = adjustl(adjustr(yzs)//'.txt')
#endif
#ifdef SFX_FA
    cfilein_fa  = adjustl(adjustr(yzs)//'.fa')
#endif
#ifdef SFX_LFI
    cfilein_lfi = adjustl(yzs)
#endif
    CALL init_io_surf_n(dtco, u, yfiletype,'FULL  ','SURF  ','READ ')
  ENDIF     
!   
  CALL read_surf(yfiletype,'ZS',u%XZS(:),iresp) 
  CALL read_sso_n(u%NSIZE_FULL, u%XSEA, uss, yfiletype)
!
  CALL end_io_surf_n(yfiletype)
!
  IF (lhook) CALL dr_hook('PGD_OROGRAPHY',1,zhook_handle)
  RETURN
!

ELSE
  !
  ALLOCATE(nsize_all(u%NDIM_FULL,1))
  ALLOCATE(xext_all(u%NDIM_FULL,2))
  ALLOCATE(xall(u%NDIM_FULL,2,1))  
  nsize_all(:,1) = 0.
  xext_all(:,1) = -99999.
  xext_all(:,2) = 99999.
  xall(:,:,1) = 0.  
  !
  !-------------------------------------------------------------------------------
!
!*    6.      Averages the field
!             ------------------
!
  CALL treat_field(ug, u, uss, &
                  hprogram,'SURF  ',yfiletype,'A_OROG',yzs, 'ZS                  ' )  
!

  DEALLOCATE(xsumval )
  !
ENDIF  
!
IF (.NOT.ALLOCATED(nsize)) THEN
  ALLOCATE(nsize(nl,1))
  nsize(:,1) = 0
ENDIF
!
!-------------------------------------------------------------------------------
!
!*    7.      Coherence with land sea mask
!             ----------------------------
!
WHERE (u%XSEA(:)==1. .AND. nsize(:,1)==0) nsize(:,1) = -1
!
!-------------------------------------------------------------------------------
!
!*    8.      Interpolation if some points are not initialized (no data for these points)
!             ------------------------------------------------
!
! note that if no orography data exists near points that need to be defined,
! these points are probably small isolated islands, and a default value of 1m is assumed.
!
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize(:,1),uss%XAVG_ZS,   'average orography',pdef=1.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize(:,1),uss%XSIL_ZS,   'silhouette orography',pdef=1.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize(:,1),uss%XMIN_ZS,   'minimum orography',pdef=1.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize(:,1),uss%XMAX_ZS,   'maximum orography',pdef=1.)
!
iflag(:) = nsize(:,1)
WHERE (nsize(:,1)==1) iflag(:) = 0 ! only 1 data point was not enough for standard deviation
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XSSO_STDEV,'standard deviation of orography',pdef=0.)
!
!-------------------------------------------------------------------------------
!
!*    9.      Coherence with land sea mask
!             ----------------------------
!
uss%XAVG_ZS   (:) = uss%XAVG_ZS   (:) * (1. - u%XSEA(:))
uss%XSIL_ZS   (:) = uss%XSIL_ZS   (:) * (1. - u%XSEA(:))
!
WHERE (u%XSEA(:)==1.)
  uss%XSSO_STDEV(:) = xundef
END WHERE
!
WHERE (u%XWATER(:)==1.)
  uss%XSSO_STDEV(:) = 0.
END WHERE
!
WHERE(u%XSEA(:)>0.)
  uss%XMIN_ZS(:) = 0.
END WHERE
!
WHERE(u%XSEA(:)==1.)
  uss%XMAX_ZS(:) = 0.
END WHERE
!
!* slightly modifies the orography values when there are by coincidence equal to
!  default value.
!
WHERE (uss%XAVG_ZS==xundef) uss%XAVG_ZS = uss%XAVG_ZS + zeps
WHERE (uss%XSIL_ZS==xundef) uss%XSIL_ZS = uss%XSIL_ZS + zeps
WHERE (uss%XMIN_ZS==xundef) uss%XMIN_ZS = uss%XMIN_ZS + zeps
WHERE (uss%XMAX_ZS==xundef) uss%XMAX_ZS = uss%XMAX_ZS + zeps
!
!-------------------------------------------------------------------------------
!
!*   10.      Choice of orography
!             -------------------
!
SELECT CASE (corogtype)
  CASE ('AVG')
    u%XZS(:) = uss%XAVG_ZS(:)
  CASE ('ENV')
    u%XZS(:) = uss%XAVG_ZS(:)
    WHERE (u%XSEA(:)<1.) u%XZS(:) = uss%XAVG_ZS(:) + xenv * uss%XSSO_STDEV
  CASE ('SIL')
    u%XZS(:) = uss%XSIL_ZS(:)
  CASE ('MAX')
    u%XZS(:) = uss%XMAX_ZS(:)
  CASE DEFAULT
    CALL abor1_sfx('PGD_OROGRAPHY: OROGRAPHY TYPE NOT SUPPORTED '//corogtype)
END SELECT
!
!-------------------------------------------------------------------------------
!
!*   12.      Subgrid scale orography characteristics
!             ---------------------------------------
!
 CALL sso(u, ug, uss, gsso, gsso_anis)
!
iflag(:) = nsize(:,1)
WHERE(.NOT. gsso(:))                 iflag(:) = 0
WHERE(u%XSEA(:)==1. .AND. iflag(:)==0) iflag(:) = -1
!
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XSSO_DIR,  'subgrid orography direction',pdef=0.)
!
IF (lexplicit_slope) THEN
  CALL explicit_slope(ug, u%NDIM_FULL, u%XZS, uss%XSSO_SLOPE) 
ELSEIF (len_trim(yslope)==0) THEN
  CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XSSO_SLOPE,'subgrid orography slope',pdef=0.)  
END IF
!
iflag(:) = nsize(:,1)
WHERE(.NOT. gsso_anis(:))            iflag(:) = 0
WHERE(u%XSEA(:)==1. .AND. iflag(:)==0) iflag(:) = -1
!
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XSSO_ANIS, 'subgrid orography anisotropy',pdef=0.)
!
WHERE (u%XSEA(:)==1.)
  uss%XSSO_ANIS (:) = xundef
  uss%XSSO_DIR  (:) = xundef
  uss%XSSO_SLOPE(:) = xundef
END WHERE
!
WHERE (u%XWATER(:)==1.)
  uss%XSSO_ANIS (:) = 1.
  uss%XSSO_DIR  (:) = 0.
  uss%XSSO_SLOPE(:) = 0.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*   13.      Subgrid scale orography roughness
!             ---------------------------------
!
 CALL subscale_aos(u, ug, uss, gz0effi, gz0effj)
!
iflag(:) = nsize(:,1)
WHERE(.NOT. gz0effi(:))              iflag(:) = 0
WHERE(u%XSEA(:)==1. .AND. iflag(:)==0) iflag(:) = -1
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XAOSIP, 'subgrid orography A/S, direction i+',pdef=0.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XAOSIM, 'subgrid orography A/S, direction i-',pdef=0.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XHO2IP, 'subgrid orography h/2, direction i+',pdef=0.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XHO2IM, 'subgrid orography h/2, direction i-',pdef=0.)
!
iflag(:) = nsize(:,1)
WHERE(.NOT. gz0effj(:))              iflag(:) = 0
WHERE(u%XSEA(:)==1. .AND. iflag(:)==0) iflag(:) = -1
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XAOSJP, 'subgrid orography A/S, direction j+',pdef=0.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XAOSJM, 'subgrid orography A/S, direction j-',pdef=0.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XHO2JP, 'subgrid orography h/2, direction j+',pdef=0.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XHO2JM, 'subgrid orography h/2, direction j-',pdef=0.)
!
WHERE (u%XSEA(:)==1.)
  uss%XHO2IP(:) = xundef
  uss%XHO2IM(:) = xundef
  uss%XHO2JP(:) = xundef
  uss%XHO2JM(:) = xundef
  uss%XAOSIP(:) = xundef
  uss%XAOSIM(:) = xundef
  uss%XAOSJP(:) = xundef
  uss%XAOSJM(:) = xundef
END WHERE
!
WHERE (u%XWATER(:)==1.)
  uss%XHO2IP(:) = 0.
  uss%XHO2IM(:) = 0.
  uss%XHO2JP(:) = 0.
  uss%XHO2JM(:) = 0.
  uss%XAOSIP(:) = 0.
  uss%XAOSIM(:) = 0.
  uss%XAOSJP(:) = 0.
  uss%XAOSJM(:) = 0.
END WHERE
!-------------------------------------------------------------------------------
DEALLOCATE(nsize    )
IF (lhook) CALL dr_hook('PGD_OROGRAPHY',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE pgd_orography