init_slope_param.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_slope_param (UG, PGRID_FULL_PAR, PZS, KI, PLAT)
!
!!**** *INIT_SLOPE_PARAM
!!                         Compute physiographic field necessary to account for
!!                         slope and aspect effects on short-wave radiations
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!   
!
!!    EXTERNAL
!!    --------
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!   
!!
!!    AUTHOR
!!    ------
!!
!!    V. Vionnet        Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!
!!    Original    04/11
!!
!!
!!    Modifications Matthieu Lafaysse 04/14 :
!!       * Computations of XSLOPANG, XSLOPAZI and XSURF_TRIANGLE in this routine (optimization)
!!       * Call to GET_GRID_DIM moved in get_size_parallel to initialize NIX and NIY
!!       * Modifications linked to parallelization : gather full domain on all MPI threads
!!       * Solar azimuth is now counted from North clockwise : N --> 0 ; E --> pi/2 ; S --> pi ; W --> -pi/2
!!         for consistency with solar azimuth.
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
!
USE modd_sfx_grid_n, ONLY : grid_t
!
USE modi_get_mesh_dim
!
USE modd_surfex_mpi, ONLY : nrank, npio, nproc, ncomm, nsize, nindex, nsize_task

USE modd_slope_effect
!
IMPLICIT NONE
!
#ifdef SFX_MPI
include "mpif.h"
#endif

!*    0.1    Declaration of arguments
!            ------------------------
!
!
TYPE(grid_t), INTENT(INOUT) :: UG
REAL, DIMENSION(:), INTENT(IN) :: PGRID_FULL_PAR
!
INTEGER,              INTENT(IN)  :: KI        ! number of points
REAL, DIMENSION(KI),   INTENT(IN) :: PZS      ! orography of this MPI thread (or total domain if Open MP)
REAL,DIMENSION(:),INTENT(IN):: PLAT ! latitudes
!
!
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
LOGICAL :: GRECT    ! true when grid is rectangular

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

REAL, DIMENSION(:),ALLOCATABLE :: ZZS1D_FULL !orography of total domain
REAL, DIMENSION(:,:), ALLOCATABLE :: ZZS ! orography in a 2D array

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

REAL,    PARAMETER :: XPI=4.*atan(1.)  ! Pi
INTEGER, PARAMETER :: JPHEXT = 1 ! number of points around the physical domain

REAL                   :: ZDZSDX   ! slope in X and Y direction
REAL                   :: ZDZSDY   ! of a triangle surface

INTEGER :: IIB, IIE, IJB, IJE
INTEGER :: JI, JJ
INTEGER :: JT       ! loop counter (triangles)
INTEGER :: INFOMPI
INTEGER,DIMENSION(:),ALLOCATABLE::IDISPLS
INTEGER::JPOINT,JPROC
INTEGER::IRANK,INEXTRANK,IPOS
INTEGER::IINDY

!-------------------------------------------------------------------------------
!
!*    1.1     Gets the geometry of the grid
!            -----------------------------
!


! Toutes les threads MPI doivent connaître le champ d'altitude complet
ALLOCATE(zzs1d_full(nix*niy))

#ifdef SFX_MPI

IF (nproc>1) THEN

  ALLOCATE(idispls(0:nproc-1))
  idispls(0)=0

  !IRANK est le rang du thread associé au 1er point
  ipos=1
  irank=nindex(ipos)
  idispls(irank)=0

  ! Attention cette méthode ne marche que si NINDEX est continu (découpage linéaire)
  DO jproc=1,nproc-1
     !INEXTRANK est le rang du thread suivant : on décale IPOS du nb de points affectés au thread précedent
     !                                        et on ajoute à IDISPLS(INEXTRANK) l'espace correspondant
     ipos=ipos+nsize_task(irank)
     inextrank=nindex(ipos)
     idispls(inextrank)=idispls(irank)+nsize_task(irank)*kind(pzs)/4
     irank=inextrank
  END DO

  CALL mpi_allgatherv(pzs,SIZE(pzs)*kind(pzs)/4,mpi_real,zzs1d_full,          &
                      nsize_task(:)*kind(pzs)/4,idispls,mpi_real,ncomm,infompi)

ELSE
  zzs1d_full=pzs
ENDIF

#else
  zzs1d_full=pzs
#endif

nnx=nix+2
nny=niy+2

!

!*    1.2    Grid dimension (meters)
!            -----------------------
!
ALLOCATE(zdx(nix*niy))
ALLOCATE(zdy(nix*niy))

IF (nrank==npio) THEN

  CALL get_mesh_dim(ug%CGRID,SIZE(pgrid_full_par),nix*niy,pgrid_full_par,zdx,zdy,ug%XMESH_SIZE)
  
ENDIF

#ifdef SFX_MPI
IF (nproc>1) THEN
  CALL mpi_bcast(zdx,SIZE(zdx)*kind(zdx)/4,mpi_real,npio,ncomm,infompi)
  CALL mpi_bcast(zdy,SIZE(zdy)*kind(zdy)/4,mpi_real,npio,ncomm,infompi)  
END IF
#endif
!
!-------------------------------------------------------------------------------
!
!*    2.     If grid is not rectangular, nothing is done
!            -------------------------------------------
!
IF (SIZE(zzs1d_full) /= nix * niy) stop "BIG PROBLEM WITH SIZE"
!
!-------------------------------------------------------------------------------
!
!*    3.1     Grid rectangular: orography is put in a 2D array
!            ------------------------------------------------
!
ALLOCATE(zzs(nix,niy))
ALLOCATE(xzsl(nnx,nny))
!
!RJ: next one does not work with NPROC>4
!RJ 'Fortran runtime error: Index '13' of dimension 1 of array 'plat' above upper bound of 12'
! 2d grid should be increasing in latitude
lrevertgrid=(plat(1)>plat(1+nix))

IF (lrevertgrid) THEN
  DO jy=1,niy
    iindy=niy-jy+1
    DO jx=1,nix
      zzs(jx,iindy) = zzs1d_full( jx + (jy-1)*nix ) 
    END DO
  END DO
ELSE
  DO jy=1,niy
    DO jx=1,nix
      zzs(jx,jy) = zzs1d_full( jx + (jy-1)*nix ) 
    END DO
  END DO
ENDIF

xzsl(2:nnx-1,2:nny-1) = zzs(:,:)
xzsl(1,:) = xzsl(2,:)
xzsl(nnx,:) = xzsl(nnx-1,:)
xzsl(:,1) = xzsl(:,2)
xzsl(:,nny) = xzsl(:,nny-1)

!------------------------------------------------------------------------------------------
!
!*    3.2.    Orography of SW corner of grid meshes
!     -------------------------------------
!
ALLOCATE(xzs_xy(nnx,nny))
xzs_xy(2:nnx,2:nny) = 0.25*(xzsl(2:nnx,2:nny)  +xzsl(1:nnx-1,2:nny) +&
                            xzsl(2:nnx,1:nny-1)+xzsl(1:nnx-1,1:nny-1))
!
xzs_xy(1,:) = xzs_xy(2,:)
xzs_xy(:,1) = xzs_xy(:,2)
!
!*    3.3     Initialize Grid meshes
!      -----------
!
ALLOCATE(xxhat(nnx))
ALLOCATE(xyhat(nny))

xxhat(1)=zdx(1)
xyhat(1)=zdy(1)
DO jx=2,nnx
  xxhat(jx) = xxhat(jx-1)+zdx(jx)
END DO
DO jy=2,nny
  xyhat(jy) = xyhat(jy-1)+zdy(jy)
END DO

DEALLOCATE(zzs)

!-------------------------------------------------------------------------------
!
!*       4.    Compute slope angles
!              -------------------

iib= 1+jphext
iie=SIZE(xxhat)-jphext
ijb=1+jphext
ije=SIZE(xyhat)-jphext

ALLOCATE(xslopang(nnx,nny,4))
ALLOCATE(xslopazi(nnx,nny,4))
ALLOCATE(xsurf_triangle(nnx,nny,4))

!          4.1 Loop on 4 triangles
!
!------------------------------------------------------------------------------------------
DO jt=1,4
!
  DO jj=ijb,ije
    DO ji=iib,iie
!
!-------------------------------------------------------------------------------
!
!*       4.2  Compute local slope
!              --------------------------------------------
!
!* slopes in x and y
!
      SELECT CASE (jt)
        CASE (1)
          zdzsdx=(    2.* xzsl(ji,jj)                   &
                   - (xzs_xy(ji,jj)+xzs_xy(ji,jj+1)) )    &
                 / (xxhat(ji+1)-xxhat(ji))  
          zdzsdy=(  xzs_xy(ji,jj+1) - xzs_xy(ji,jj) )     &
                 / (xyhat(jj+1)-xyhat(jj))  
        CASE (2)
           zdzsdx=(  xzs_xy(ji+1,jj+1) -xzs_xy(ji,jj+1))  &
                 / (xxhat(ji+1)-xxhat(ji))  
           zdzsdy=(  (xzs_xy(ji+1,jj+1)+xzs_xy(ji,jj+1))  &
                     - 2.* xzsl(ji,jj) )                  &
                 / (xyhat(jj+1)-xyhat(jj))  
        CASE (3)
          zdzsdx=(  (xzs_xy(ji+1,jj)+xzs_xy(ji+1,jj+1))   &
                   - 2.* xzsl(ji,jj)                    )  &
                 / (xxhat(ji+1)-xxhat(ji))  
          zdzsdy=(  xzs_xy(ji+1,jj+1) - xzs_xy(ji+1,jj) ) &
                 / (xyhat(jj+1)-xyhat(jj))  
        CASE (4)
           zdzsdx=(  xzs_xy(ji+1,jj) - xzs_xy(ji,jj) )    &
                 / (xxhat(ji+1)-xxhat(ji))  
           zdzsdy=(  2.* xzsl(ji,jj)                       &
                   - (xzs_xy(ji+1,jj)+xzs_xy(ji,jj)) )    &
                 / (xyhat(jj+1)-xyhat(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))
!* slope angles
!
      xslopang(ji,jj,jt) = atan(sqrt(zdzsdx**2+zdzsdy**2))
      xslopazi(ji,jj,jt) = 1.5*xpi - atan2( zdzsdy, zdzsdx + sign(1.e-30,zdzsdx) )
!

! surface of each triangle
      xsurf_triangle(ji,jj,jt)=sqrt(1. + zdzsdx**2 + zdzsdy**2)

!
!* normalizes received radiation by the surface of the triangle to obtain
!  radiation representative of an horizontal surface.
!
!       XSSO_SURF(JI,JJ) = XSSO_SURF(JI,JJ) + 0.25*XSURF_TRIANGLE(JI,JJ,JT)
!
    END DO
  END DO
END DO

! 2d arrays for processor domain

!------------------------------------------------------------------------------------------


END SUBROUTINE init_slope_param