explicit_slope.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 explicit_slope (UG,KDIM_FULL, &
                           PZS,PSSO_SLOPE)
!     #########################################################################
!!    AUTHOR
!!    ------
!!      M. Lafaysse      * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    19/07/13
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
!
USE modi_read_and_send_mpi
USE modi_gather_and_write_mpi
!
USE modi_get_grid_dim
USE modi_get_mesh_dim

IMPLICIT NONE
!
!*       0.1   DECLARATIONS OF DUMMY ARGUMENTS :
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
!
INTEGER, INTENT(IN) :: KDIM_FULL
REAL,DIMENSION(:),INTENT(IN)::PZS ! resolved model orography
REAL,DIMENSION(:),INTENT(OUT)::PSSO_SLOPE ! resolved slope tangent



!
!
!*       0.2   DECLARATIONS OF LOCAL VARIABLES
!

INTEGER :: IX       ! number of points in X direction
INTEGER :: IY       ! number of points in Y direction

INTEGER :: INNX !  number of points in X direction for large domain
INTEGER :: INNY !  number of points in Y direction for large domain

LOGICAL::GRECT=.true.

INTEGER :: JX       ! loop counter
INTEGER :: JY       ! loop counter

REAL, DIMENSION(:,:), ALLOCATABLE :: ZMAP         ! map factor
REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS ! orography in a 2D array
REAL, DIMENSION(:,:), ALLOCATABLE :: ZSSO_SLOPE ! explicit slope in a 2D array
REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS_XY ! orography at southwest corner of the mesh
REAL, DIMENSION(:,:),ALLOCATABLE :: ZZSL   ! orography in a 2D array
REAL, DIMENSION(:),ALLOCATABLE :: ZXHAT  ! X coordinate
REAL, DIMENSION(:),ALLOCATABLE :: ZYHAT  ! Y coordinate


REAL,DIMENSION(:), ALLOCATABLE :: ZDX        ! grid mesh size in x direction
REAL,DIMENSION(:), ALLOCATABLE  :: ZDY       ! grid mesh size in y direction

REAL, DIMENSION(:), ALLOCATABLE :: ZZS0
REAL, DIMENSION(:), ALLOCATABLE :: ZSSO_SLOPE0
!
! parameters
REAL,    PARAMETER :: XPI=4.*atan(1.)  ! Pi
INTEGER, PARAMETER :: JPHEXT = 1 ! number of points around the physical domain
!
INTEGER :: IIB, IIE, IJB, IJE
INTEGER :: JI, JJ, JB
INTEGER :: JT
!
REAL                   :: ZDZSDX   ! slope in X and Y direction
REAL                   :: ZDZSDY   ! of a triangle surface
REAL                   :: ZSURF    ! surface of 4 triangles
!



!-------------------------------------------------------------------------------
!
!*    1.1     Gets the geometry of the grid
!            -----------------------------
!
 CALL get_grid_dim(ug%G%CGRID,SIZE(ug%XGRID_FULL_PAR),ug%XGRID_FULL_PAR,grect,ix,iy)
!
innx=ix+2
inny=iy+2

!

!*    1.2    Grid dimension (meters)
!            -----------------------
!
ALLOCATE(zdx(ix*iy))
ALLOCATE(zdy(ix*iy))

 CALL get_mesh_dim(ug%G%CGRID,SIZE(ug%XGRID_FULL_PAR),ix*iy,ug%XGRID_FULL_PAR,zdx,zdy,ug%G%XMESH_SIZE)

!
!*    2.     If grid is not rectangular, nothing is done
!            -------------------------------------------
!
 ALLOCATE(zzs0(kdim_full))
 CALL gather_and_write_mpi(pzs,zzs0)
 !
!IF (.NOT. GRECT) RETURN
!
IF (SIZE(zzs0) /= ix * iy) RETURN
!
!-------------------------------------------------------------------------------
!
!*    3.1     Grid rectangular: orography is put in a 2D array
!            ------------------------------------------------
!
ALLOCATE(zzs(ix,iy))
ALLOCATE(zzsl(innx,inny))

DO jy=1,iy
  DO jx=1,ix
    zzs(jx,jy) = zzs0( jx + (jy-1)*ix ) 
  END DO
END DO

DEALLOCATE(zzs0)
!
zzsl(2:innx-1,2:inny-1) = zzs(:,:)
zzsl(1,:) = zzsl(2,:)
zzsl(innx,:) = zzsl(innx-1,:)
zzsl(:,1) = zzsl(:,2)
zzsl(:,inny) = zzsl(:,inny-1)

!------------------------------------------------------------------------------------------
!
!*    3.2.    Orography of SW corner of grid meshes
!     -------------------------------------
!

ALLOCATE(zzs_xy(innx,inny))

zzs_xy(2:innx,2:inny) = 0.25*(  zzsl(2:innx,2:inny)   + zzsl(1:innx-1,2:inny)   &
                          + zzsl(2:innx,1:inny-1) + zzsl(1:innx-1,1:inny-1) )
!
zzs_xy(1,:) = zzs_xy(2,:)
zzs_xy(:,1) = zzs_xy(:,2)

!
!*    3.3     Initialize Grid meshes
!      -----------
!
ALLOCATE(zxhat(innx))
ALLOCATE(zyhat(inny))

DO jx=1,innx
  zxhat(jx) = zdx(1)*jx
END DO
DO jy=1,inny
  zyhat(jy) = zdy(1)*jy
END DO

DEALLOCATE(zdx,zdy)

!-------------------------------------------------------------------------------
!
iib= 1+jphext
iie=innx-jphext
ijb=1+jphext
ije=inny-jphext
!
ALLOCATE(zmap(innx,inny))
zmap(:,:)=1.0
ALLOCATE(zsso_slope(ix,iy))

!-------------------------------------------------------------------------------
!
!*       1.    LOOP ON GRID MESHES
!              -------------------
!
!* discretization of the grid mesh in four triangles
!
!
DO jj=ijb,ije
  DO ji=iib,iie
    zsurf=0.
    DO jt=1,4
!
!* slopes in x and y
!
      SELECT CASE (jt)
        CASE (1)
          zdzsdx=(    2.* zzsl(ji,jj)                   &
                   - (zzs_xy(ji,jj)+zzs_xy(ji,jj+1)) )    &
                 / (zxhat(ji+1)-zxhat(ji))  * zmap(ji,jj)
          zdzsdy=(  zzs_xy(ji,jj+1) - zzs_xy(ji,jj) )     &
                 / (zyhat(jj+1)-zyhat(jj))  * zmap(ji,jj)
        CASE (2)
           zdzsdx=(  zzs_xy(ji+1,jj+1) -zzs_xy(ji,jj+1))  &
                 / (zxhat(ji+1)-zxhat(ji))  * zmap(ji,jj)
           zdzsdy=(  (zzs_xy(ji+1,jj+1)+zzs_xy(ji,jj+1))  &
                     - 2.* zzsl(ji,jj) )                  &
                 / (zyhat(jj+1)-zyhat(jj))  * zmap(ji,jj)
        CASE (3)
          zdzsdx=(  (zzs_xy(ji+1,jj)+zzs_xy(ji+1,jj+1))   &
                   - 2.* zzsl(ji,jj)                    )  &
                 / (zxhat(ji+1)-zxhat(ji))  * zmap(ji,jj)
          zdzsdy=(  zzs_xy(ji+1,jj+1) - zzs_xy(ji+1,jj) ) &
                 / (zyhat(jj+1)-zyhat(jj))  * zmap(ji,jj)
        CASE (4)
           zdzsdx=(  zzs_xy(ji+1,jj) - zzs_xy(ji,jj) )    &
                 / (zxhat(ji+1)-zxhat(ji))  * zmap(ji,jj)
           zdzsdy=(  2.* zzsl(ji,jj)                       &
                   - (zzs_xy(ji+1,jj)+zzs_xy(ji,jj)) )    &
                 / (zyhat(jj+1)-zyhat(jj))  * zmap(ji,jj)
      END SELECT
!
!
      ! If slope is higher than 60 degrees : numerical problems
      zdzsdx=min(2.0,max(-2.0,zdzsdx))
      zdzsdy=min(2.0,max(-2.0,zdzsdy))
    
     ! total surface of 4 triangles
      zsurf=zsurf+0.25*sqrt(1. + zdzsdx**2 + zdzsdy**2)

    END DO
    
    !equivalent tangent slope of a homogeneous surface with the same area
    zsso_slope(ji-jphext,jj-jphext)=sqrt(zsurf**2-1)
    
  END DO
END DO
DEALLOCATE(zzsl)
DEALLOCATE(zzs)
DEALLOCATE(zzs_xy)
DEALLOCATE(zmap)
!
ALLOCATE(zsso_slope0(kdim_full))
!
DO jy=1,iy
  DO jx=1,ix
    zsso_slope0( jx + (jy-1)*ix )=zsso_slope(jx,jy)
  END DO
END DO
!
 CALL read_and_send_mpi(zsso_slope0,psso_slope)
!
DEALLOCATE(zsso_slope0)
!
END SUBROUTINE explicit_slope