read_nam_grid_gauss.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 read_nam_grid_gauss(PGRID_FULL_PAR,KDIM_FULL,HPROGRAM,KGRID_PAR,KL,PGRID_PAR,HDIR)
!################################################################
!
!!****  *READ_NAM_GRID_GAUSS* - routine to read in namelist the horizontal grid
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004
!!      B. Decharme    2008  Comput and save the Mesh size
!!                     2013  Bug lat and lon for non rotat-strech grid
!!                           Grid corner (used with oasis)
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_surfex_mpi, ONLY : nrank, nsize_task
!
USE modd_csts,     ONLY : xpi
USE modd_surf_par, ONLY : xundef
! 
USE mode_pos_surf
!
USE modi_open_namelist
USE modi_close_namelist
USE modi_get_luout
!
USE mode_gridtype_gauss
!
USE eggangles , ONLY : p_asin
!
USE modi_read_and_send_mpi
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
REAL, DIMENSION(:), POINTER :: PGRID_FULL_PAR
INTEGER, INTENT(IN) :: KDIM_FULL
!
 CHARACTER(LEN=6),           INTENT(IN)    :: HPROGRAM   ! calling program
INTEGER,                    INTENT(INOUT) :: KGRID_PAR  ! size of PGRID_PAR
INTEGER,                    INTENT(OUT)   :: KL         ! number of points
REAL, DIMENSION(KGRID_PAR), INTENT(OUT)   :: PGRID_PAR  ! parameters defining this grid
 CHARACTER(LEN=1), INTENT(IN) :: HDIR
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER :: ILUOUT ! output listing logical unit
INTEGER :: ILUNAM ! namelist file  logical unit
REAL,    DIMENSION(:), ALLOCATABLE :: ZLAT_XY, ZLAT_XY0 ! pseudo-latitudes
REAL,    DIMENSION(:), ALLOCATABLE :: ZLON_XY, ZLON_XY0 ! pseudo-longitudes
REAL,    DIMENSION(:), ALLOCATABLE :: ZLAT, ZLAT0    ! latitudes
REAL,    DIMENSION(:), ALLOCATABLE :: ZLON, ZLON0    ! longitudes
REAL,    DIMENSION(:), ALLOCATABLE :: ZMESH_SIZE, ZMESH_SIZE0 ! Mesh size
!                                                                              _____ Sup
REAL,    DIMENSION(:),   ALLOCATABLE :: ZLATSUP, ZLATSUP0    ! Grid corner Latitude    |     |
REAL,    DIMENSION(:),   ALLOCATABLE :: ZLONSUP, ZLONSUP0     ! Grid corner Longitude   |     |
REAL,    DIMENSION(:),   ALLOCATABLE :: ZLATINF, ZLATINF0     ! Grid corner Latitude    |_____|
REAL,    DIMENSION(:),   ALLOCATABLE :: ZLONINF, ZLONINF0     ! Grid corner Longitude  Inf
!
REAL,    DIMENSION(:),   ALLOCATABLE :: ZXINF    ! pseudo-longitude western limit of grid mesh
REAL,    DIMENSION(:),   ALLOCATABLE :: ZXSUP    ! pseudo-longitude eastern limit of grid mesh
REAL,    DIMENSION(:),   ALLOCATABLE :: ZYINF    ! pseudo-latitude southern limit of grid mesh
REAL,    DIMENSION(:),   ALLOCATABLE :: ZYSUP    ! pseudo-latitude northern limit of grid mesh
!
!*       0.3   Declarations of namelist
!              ------------------------
!
INTEGER :: NDGLG    ! number of pseudo-latitudes
REAL    :: RMUCEN   ! sine of the latitude of the rotated pole
REAL    :: RLOCEN   ! longitude of the rotated pole (radian)
REAL    :: RSTRET   ! stretching factor (must be greater than or equal to 1)
INTEGER, DIMENSION(1000) :: NRGRI ! number of pseudo-longitudes on each
                                  ! pseudo-latitude circle on pseau
                                  ! northern hemisphere (starting from
                                  ! the rotated pole)
!
REAL    :: ZLAPO    ! latitude  of the rotated pole (deg)
REAL    :: ZLOPO    ! longitude of the rotated pole (deg)
REAL    :: ZCODIL   ! stretching factor (must be greater than or equal to 1)
INTEGER                            :: ITYP   ! type of transform (0 --> no rotation, 1 otherwise)
INTEGER                            :: INLATI ! number of latitudes
INTEGER, DIMENSION(:), ALLOCATABLE :: INLOPA ! number of pseudo-longitudes on each
                                             ! pseudo-latitude circle

INTEGER :: JSTGLO

!
REAL, DIMENSION(:), POINTER :: ZGRID_PAR
!
LOGICAL :: GFOUND
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
NAMELIST/namdim/ndglg
NAMELIST/namgem/rmucen, rlocen, rstret
NAMELIST/namrgri/nrgri
!
!------------------------------------------------------------------------------
!
!*       1.    Default values
! 
IF (lhook) CALL dr_hook('READ_NAM_GRID_GAUSS',0,zhook_handle)
ndglg = 0
rmucen = 1.
rlocen = xpi
rstret = 1.
!
nrgri(:) = 0
!------------------------------------------------------------------------------
!
!*       2.    opening of namelist
! 
 CALL get_luout(hprogram,iluout)
!
IF (hdir/='H') THEN
  !
  CALL open_namelist(hprogram,ilunam)
  !
  !---------------------------------------------------------------------------
  !
  !*       3.    Reading of projection parameters
  !              --------------------------------
  !
  CALL posnam(ilunam,'NAMGEM',gfound,iluout)
  IF (gfound) READ(unit=ilunam,nml=namgem)
  !
  IF (rstret<1.) THEN
    WRITE(iluout,*) '****************************************************'
    WRITE(iluout,*) 'stretching factor RSTRET for the Gaussian grid'
    WRITE(iluout,*) 'definition must be greater than or equal to 1'
    WRITE(iluout,*) 'You have set RSTRET=', rstret
    WRITE(iluout,*) 'Please modify the value of RSTRET in namelist NAMGEM'
    WRITE(iluout,*) '****************************************************'
    CALL abor1_sfx('READ_NAM_GRID_GAUSS: STRETCHING FACTOR MUST BE >= 1.')
  END IF
  !
  zlapo = 180. / xpi * p_asin(rmucen)
  zlopo = 180. / xpi * rlocen
  !
  zcodil = rstret
  !
  !---------------------------------------------------------------------------
  !
  !*       4.    Reading parameters of the grid
  !              ------------------------------
  !
  CALL posnam(ilunam,'NAMDIM',gfound,iluout)
  IF (gfound) READ(unit=ilunam,nml=namdim)
  CALL posnam(ilunam,'NAMRGRI',gfound,iluout)
  IF (gfound) READ(unit=ilunam,nml=namrgri)
  !
  inlati = ndglg
  ALLOCATE(inlopa(inlati))
  inlopa(1:inlati/2) = nrgri(1:inlati/2)
  inlopa(inlati/2+1:inlati) = nrgri(inlati/2:1:-1)
  !
  !---------------------------------------------------------------------------
  CALL close_namelist(hprogram,ilunam)
  !---------------------------------------------------------------------------
  !
  !*       5.    Computes pseudo-latitudes and pseudo-longitudes of all points
  !              -------------------------------------------------------------
  !
  !* number of points
  kl = sum(inlopa)
  !
  !
  !* type of transform
  IF (zlapo>89.99 .AND. abs(zlopo)<0.00001) THEN
    ityp=0
  ELSE
    ityp=1
  ENDIF
  !
  ALLOCATE(zlat_xy(kl))
  ALLOCATE(zlon_xy(kl))
  zlat_xy(:) = xundef
  zlon_xy(:) = xundef
  !
  CALL comp_gridtype_gauss(inlati,inlopa,kl,ityp,zlat_xy,zlon_xy)
  !
  !---------------------------------------------------------------------------
  !
  !*       6.    Computes latitudes and longitudes
  !              ---------------------------------
  !
  !* all points are used
  ALLOCATE(zlat(kl))
  ALLOCATE(zlon(kl))
  !
  zlat(:) = xundef
  zlon(:) = xundef
  !
  IF(zcodil==1.0.AND.ityp==0)THEN
    zlon(:)=zlon_xy(:)
    zlat(:)=zlat_xy(:)
  ELSE
    CALL latlon_gauss(zlon_xy,zlat_xy,kl,zlopo,zlapo,zcodil,zlon,zlat)
  ENDIF
  !
  !---------------------------------------------------------------------------
  !
  !*       7.    Computes grid corner latitudes and longitudes
  !              ---------------------------------------------
  !
  ALLOCATE(zxinf(kl))
  ALLOCATE(zyinf(kl))
  ALLOCATE(zxsup(kl))
  ALLOCATE(zysup(kl))
  !
  ALLOCATE(zloninf(kl))
  ALLOCATE(zlatinf(kl))
  ALLOCATE(zlonsup(kl))
  ALLOCATE(zlatsup(kl))
  !
  zxinf(:) = xundef
  zyinf(:) = xundef
  zxsup(:) = xundef
  zysup(:) = xundef
  zloninf(:) = xundef
  zlatinf(:) = xundef
  zlonsup(:) = xundef
  zlatsup(:) = xundef
  !
  CALL gauss_grid_limits(inlati,inlopa,zxinf,zxsup,zyinf,zysup)
  !
  IF(zcodil==1.0.AND.ityp==0)THEN
    zloninf(:) = zxinf(:)
    zlatinf(:) = zyinf(:)
    zlonsup(:) = zxsup(:)
    zlatsup(:) = zysup(:)
  ELSE        
    CALL latlon_gauss(zxinf,zyinf,kl,zlopo,zlapo,zcodil,zloninf,zlatinf)  
    CALL latlon_gauss(zxsup,zysup,kl,zlopo,zlapo,zcodil,zlonsup,zlatsup)  
  ENDIF
  !
  DEALLOCATE(zxinf)
  DEALLOCATE(zyinf)
  DEALLOCATE(zxsup)
  DEALLOCATE(zysup)
  !
  !---------------------------------------------------------------------------
  !
  !*       8.    Computes mesh size
  !              ---------------------------------
  !
  ALLOCATE(zmesh_size(kl))
  zmesh_size(:) = xundef
  !
  CALL mesh_size_gauss(kl,inlati,inlopa,zlapo,zlopo,zcodil,&
                               zlat_xy,zlon,zlat,zmesh_size)  
!
ELSE
  !
  ALLOCATE(zlon0(kdim_full),zlat0(kdim_full),zlon_xy0(kdim_full),zlat_xy0(kdim_full),&
           zmesh_size0(kdim_full),zloninf0(kdim_full),zlatinf0(kdim_full),&
           zlonsup0(kdim_full),zlatsup0(kdim_full))
  !
 CALL get_gridtype_gauss(pgrid_full_par,inlati)
  !
  ALLOCATE(inlopa(inlati))  
  !
  CALL get_gridtype_gauss(pgrid_full_par,knlati=inlati,plapo=zlapo,&
                          plopo=zlopo,pcodil=zcodil,knlopa=inlopa, &
                          plat=zlat0,plon=zlon0,plat_xy=zlat_xy0,&
                          plon_xy=zlon_xy0,pmesh_size=zmesh_size0,&
                          ploninf=zloninf0,platinf=zlatinf0,&
                          plonsup=zlonsup0,platsup=zlatsup0)
  !
  kl = nsize_task(nrank)
  ALLOCATE(zlat(kl),zlon(kl),zlat_xy(kl),zlon_xy(kl),zmesh_size(kl))
  ALLOCATE(zloninf(kl))
  ALLOCATE(zlatinf(kl))
  ALLOCATE(zlonsup(kl))
  ALLOCATE(zlatsup(kl))
  !
  CALL read_and_send_mpi(zlon0,zlon)
  CALL read_and_send_mpi(zlat0,zlat)
  CALL read_and_send_mpi(zlon_xy0,zlon_xy)
  CALL read_and_send_mpi(zlat_xy0,zlat_xy)
  CALL read_and_send_mpi(zmesh_size0,zmesh_size)
  CALL read_and_send_mpi(zloninf0,zloninf)
  CALL read_and_send_mpi(zlatinf0,zlatinf)
  CALL read_and_send_mpi(zlonsup0,zlonsup)
  CALL read_and_send_mpi(zlatsup0,zlatsup)
  !
  DEALLOCATE(zlon0,zlat0,zlon_xy0,zlat_xy0,zmesh_size0,zloninf0,zlatinf0,&
          zlonsup0,zlatsup0)    
  !
ENDIF
!---------------------------------------------------------------------------
!
!*       9.    All this information stored into pointer PGRID_PAR
!              --------------------------------------------------
!
 CALL put_gridtype_gauss(zgrid_par,inlati,zlapo,zlopo,zcodil,inlopa, &
                          kl,zlat,zlon,zlat_xy,zlon_xy,zmesh_size,   &
                          zloninf,zlatinf,zlonsup,zlatsup            )  
!
!---------------------------------------------------------------------------
!
!*       9.    All this information stored into pointer PGRID_PAR
!              --------------------------------------------------
!
DEALLOCATE(zlat)
DEALLOCATE(zlon)
DEALLOCATE(zlat_xy)
DEALLOCATE(zlon_xy)
DEALLOCATE(inlopa)
DEALLOCATE(zmesh_size)
DEALLOCATE(zlatinf)
DEALLOCATE(zloninf)
DEALLOCATE(zlatsup)
DEALLOCATE(zlonsup)
!
!---------------------------------------------------------------------------
!
!* 1st call : initializes dimension
!
IF (kgrid_par==0) THEN
  kgrid_par = SIZE(zgrid_par)
!
ELSE
!
!* 2nd call : initializes grid array
!
  pgrid_par(:) = 0.
  pgrid_par(:) = zgrid_par
END IF
!
DEALLOCATE(zgrid_par)
IF (lhook) CALL dr_hook('READ_NAM_GRID_GAUSS',1,zhook_handle)

!
!---------------------------------------------------------------------------
!
END SUBROUTINE read_nam_grid_gauss