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(HPROGRAM,KGRID_PAR,KL,PGRID_PAR)
!################################################################
!
!!****  *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_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 yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
 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
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
INTEGER :: iluout ! output listing logical unit
INTEGER :: ilunam ! namelist file  logical unit
REAL,    DIMENSION(:),   ALLOCATABLE :: zlat_xy     ! pseudo-latitudes
REAL,    DIMENSION(:),   ALLOCATABLE :: zlon_xy     ! pseudo-longitudes
REAL,    DIMENSION(:),   ALLOCATABLE :: zlat        ! latitudes
REAL,    DIMENSION(:),   ALLOCATABLE :: zlon        ! longitudes
REAL,    DIMENSION(:),   ALLOCATABLE :: zmesh_size  ! Mesh size
!                                                                              _____ Sup
REAL,    DIMENSION(:),   ALLOCATABLE :: zlatsup     ! Grid corner Latitude    |     |
REAL,    DIMENSION(:),   ALLOCATABLE :: zlonsup     ! Grid corner Longitude   |     |
REAL,    DIMENSION(:),   ALLOCATABLE :: zlatinf     ! Grid corner Latitude    |_____|
REAL,    DIMENSION(:),   ALLOCATABLE :: zloninf     ! 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)
!
 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
!
!---------------------------------------------------------------------------
!
!*       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)  
!
!---------------------------------------------------------------------------
!
!*       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)
DEALLOCATE(zxinf)
DEALLOCATE(zyinf)
DEALLOCATE(zxsup)
DEALLOCATE(zysup)
!
!---------------------------------------------------------------------------
!
!* 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