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 (DGU, DTCO, UG, U, USS, &
                                hprogram,psea,pwater,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_diag_surf_atm_n, ONLY : diag_surf_atm_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_surf_atm_sso_n, ONLY : surf_atm_sso_t
!
!
USE modd_pgd_grid,       ONLY : nl, cgrid, xgrid_par
USE modd_pgdwork,        ONLY : xsumval, xsumval2, nsize, xssqo, lssqo, nsso
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_ini_ssowork
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(diag_surf_atm_t), INTENT(INOUT) :: dgu
TYPE(data_cover_t), INTENT(INOUT) :: dtco
TYPE(surf_atm_grid_t), INTENT(INOUT) :: ug
TYPE(surf_atm_t), INTENT(INOUT) :: u
TYPE(surf_atm_sso_t), INTENT(INOUT) :: uss
!
 CHARACTER(LEN=6),     INTENT(IN)  :: hprogram ! program calling
REAL, DIMENSION(:),   INTENT(IN)  :: psea     ! sea  fraction
REAL, DIMENSION(:),   INTENT(IN)  :: pwater   ! lake fraction
 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
!             --------------------------
!
ALLOCATE(nsize(nl))
ALLOCATE(xsumval(nl))
ALLOCATE(xsumval2(nl))
!
nsize(:) = 0.
xsumval(:) = 0.
xsumval2(:) = 0.
!
 CALL ini_ssowork
!
!-------------------------------------------------------------------------------
!
!*    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(u, &
                       hprogram,idim_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.

  DEALLOCATE(nsize    )
  DEALLOCATE(xsumval  )
  DEALLOCATE(xsumval2 )

  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.

  DEALLOCATE(nsize    )
  DEALLOCATE(xsumval  )
  DEALLOCATE(xsumval2 )

  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(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(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, dgu, u, &
                         yfiletype,'FULL  ','SURF  ','READ ')
  ENDIF     
!   
  CALL read_surf(&
                 yfiletype,'ZS',u%XZS(:),iresp) 
  CALL read_sso_n(&
                  u, uss, &
                  yfiletype)
!
  CALL end_io_surf_n(yfiletype)
!
  DEALLOCATE(nsize    )
  DEALLOCATE(xsumval  )
  DEALLOCATE(xsumval2 )
!
  IF (lhook) CALL dr_hook('PGD_OROGRAPHY',1,zhook_handle)
  RETURN
!

END IF
!
!-------------------------------------------------------------------------------

!
!*    6.      Averages the field
!             ------------------
!
 CALL treat_field(ug, u, uss, &
                  hprogram,'SURF  ',yfiletype,'A_OROG',yzs,  &
                   'ZS                  '                     )  
!
DEALLOCATE(xsumval  )
DEALLOCATE(xsumval2 )
!
!-------------------------------------------------------------------------------
!
!*    7.      Coherence with land sea mask
!             ----------------------------
!
WHERE (psea(:)==1. .AND. nsize(:)==0) nsize(:) = -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,uss%XAVG_ZS,   'average orography',pdef=1.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize,uss%XSIL_ZS,   'silhouette orography',pdef=1.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize,uss%XMIN_ZS,   'minimum orography',pdef=1.)
 CALL interpol_field(ug, u, &
                     hprogram,iluout,nsize,uss%XMAX_ZS,   'maximum orography',pdef=1.)
!
iflag(:) = nsize(:)
WHERE (nsize(:)==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. - psea(:))
uss%XSIL_ZS   (:) = uss%XSIL_ZS   (:) * (1. - psea(:))
!
WHERE (psea(:)==1.)
  uss%XSSO_STDEV(:) = xundef
END WHERE
!
WHERE (pwater(:)==1.)
  uss%XSSO_STDEV(:) = 0.
END WHERE
!
WHERE(psea(:)>0.)
  uss%XMIN_ZS(:) = 0.
END WHERE
!
WHERE(psea(:)==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 (psea(:)<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(ug, uss, &
          gsso,gsso_anis,psea)
!
iflag(:) = nsize(:)
WHERE(.NOT. gsso(:))                 iflag(:) = 0
WHERE(psea(:)==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%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(:)
WHERE(.NOT. gsso_anis(:))            iflag(:) = 0
WHERE(psea(:)==1. .AND. iflag(:)==0) iflag(:) = -1
!
 CALL interpol_field(ug, u, &
                     hprogram,iluout,iflag,uss%XSSO_ANIS, 'subgrid orography anisotropy',pdef=0.)
!
WHERE (psea(:)==1.)
  uss%XSSO_ANIS (:) = xundef
  uss%XSSO_DIR  (:) = xundef
  uss%XSSO_SLOPE(:) = xundef
END WHERE
!
WHERE (pwater(:)==1.)
  uss%XSSO_ANIS (:) = 1.
  uss%XSSO_DIR  (:) = 0.
  uss%XSSO_SLOPE(:) = 0.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*   13.      Subgrid scale orography roughness
!             ---------------------------------
!
 CALL subscale_aos(ug, uss, &
                   gz0effi,gz0effj,psea)
!
iflag(:) = nsize(:)
WHERE(.NOT. gz0effi(:))              iflag(:) = 0
WHERE(psea(:)==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(:)
WHERE(.NOT. gz0effj(:))              iflag(:) = 0
WHERE(psea(:)==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 (psea(:)==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 (pwater(:)==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