ol_define_dim.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 ol_define_dim (UG, KSIZE_FULL, HPROGRAM, KLUOUT, KNI, &
                          KDIM1, HUNIT1, HUNIT2, PX, PY, KDIMS,  &
                          KDDIM, HNAME_DIM, KNPATCH, KNSNLAYER, PLAT, PLON)
!     #######################################################
!!****  *OL_DEFINE_DIM* - 
!!
!!    PURPOSE
!!    -------
!!
      !!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      S. Faroux   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    06/2010 
!!      07/2011     add specific computation for IGN grid (B. Decharme)
!!      09/2015     M. Lafaysse : snow layer dimension
!-------------------------------------------------------------------------------                         
!
USE modd_surf_atm_grid_n, ONLY : surf_atm_grid_t
!
USE modd_io_surf_ol, ONLY: nmask_ign
!
USE modn_io_offline, ONLY : lwrite_coord
!
USE modi_get_grid_dim
USE modi_get_grid_coord
USE modi_get_ign_maskall
!
USE mode_gridtype_ign
USE mode_gridtype_conf_proj
USE mode_gridtype_lonlat_rot
USE mode_gridtype_lonlat_reg
USE mode_gridtype_gauss
USE mode_gridtype_cartesian
USE mode_gridtype_lonlatval
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE netcdf
!
IMPLICIT NONE
!
!
TYPE(surf_atm_grid_t), INTENT(INOUT) :: UG
!
INTEGER, INTENT(IN) :: KSIZE_FULL
!
 CHARACTER(LEN=6),  INTENT(IN)    :: HPROGRAM
INTEGER, INTENT(IN)              :: KLUOUT
INTEGER, INTENT(IN)              :: KNI
INTEGER, INTENT(OUT)             :: KDIM1
 CHARACTER(LEN=13) , DIMENSION(:), INTENT(OUT) :: HUNIT1, HUNIT2
!
REAL,DIMENSION(:), POINTER                 :: PX, PY
INTEGER, DIMENSION(:), POINTER             :: KDIMS, KDDIM
 CHARACTER(LEN=100), DIMENSION(:), POINTER :: HNAME_DIM
INTEGER, OPTIONAL, INTENT(IN)              :: KNPATCH
INTEGER, OPTIONAL, INTENT(IN)              :: KNSNLAYER
REAL, DIMENSION(:), OPTIONAL, POINTER      :: PLAT, PLON
!
REAL, DIMENSION(:), ALLOCATABLE  :: ZXFULL, ZYFULL, ZDXFULL, ZDYFULL, ZLAT_XY, ZLON_XY
REAL, DIMENSION(KNI)             :: ZXX, ZYY
REAL                             :: ZLAT0,ZLON0,ZRPK,ZBETA,ZLATOR,ZLONOR,ZCODIL,&
                                    ZLAPO,ZLOPO
 CHARACTER(LEN=3)                :: YTYPE
INTEGER, DIMENSION(:), ALLOCATABLE :: INLOPA
INTEGER                          :: JJ, ILAMBERT, IFULL, IDIM1, IL, IIMAX, IJMAX
INTEGER                          :: INDIMS, INDIMSMAX, IDIM2, INLATI
INTEGER                          :: INSNLAYER, INPATCH, ID
LOGICAL                          :: GRECT     ! T if rectangular grid
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
IF (lhook) CALL dr_hook('OL_DEFINE_DIM',0,zhook_handle)
!
NULLIFY(px)
NULLIFY(py)
!
insnlayer = 0
IF (PRESENT(knsnlayer)) insnlayer = knsnlayer
!
inpatch = 0
IF (PRESENT(knpatch)) inpatch = knpatch
!
kdim1 = 0
idim2 = 0
!
ytype='XY '
IF (ug%G%CGRID.EQ.'LONLAT REG'.OR.ug%G%CGRID.EQ.'LONLATVAL '.OR.ug%G%CGRID.EQ.'LONLAT ROT') ytype='LON'
!
IF (.NOT.lwrite_coord) THEN
  !
  IF ( ug%G%CGRID.EQ.'CONF PROJ ' .OR. ug%G%CGRID.EQ.'CARTESIAN '&
  .OR. ug%G%CGRID.EQ.'LONLAT REG' .OR. ug%G%CGRID.EQ.'IGN' ) THEN
    IF (ASSOCIATED(ug%XGRID_FULL_PAR)) THEN
      CALL get_grid_dim(ug%G%CGRID,ug%NGRID_FULL_PAR,ug%XGRID_FULL_PAR,grect,kdim1,idim2)
    ELSEIF (ASSOCIATED(ug%G%XGRID_PAR)) THEN
      CALL get_grid_dim(ug%G%CGRID,ug%G%NGRID_PAR,ug%G%XGRID_PAR,grect,kdim1,idim2)
    ENDIF
  ENDIF
  !
ENDIF
!
IF ( ug%G%CGRID.EQ.'IGN' ) CALL get_gridtype_ign(ug%XGRID_FULL_PAR,kl=il,kdimx=iimax,kdimy=ijmax)
!
indims = 0
IF ( inpatch/=0 .AND. insnlayer/=0 ) THEN
  indims = 2
ELSEIF ( inpatch/=0 .OR. insnlayer/=0 ) THEN
  indims = 1
ENDIF
!
!
! case LWRITE_COORD = F and regular grid
IF ( kdim1.NE.0 ) THEN
  indims = indims + 3
ELSEIF ( ug%G%CGRID=="GAUSS     " ) THEN
  indims = indims + 4
ELSEIF (ug%G%CGRID/="LONLATVAL ".AND.(ug%G%CGRID/="IGN       ".OR.iimax*ijmax==il) ) THEN
  indims = indims + 4
ELSE   
  indims = indims + 2
ENDIF
!
ALLOCATE(kdims(indims))
ALLOCATE(kddim(indims))
ALLOCATE(hname_dim(indims))
!
! case LWRITE_COORD = F and regular grid
IF ( kdim1.NE.0 ) THEN

  kdims(1) = kdim1
  kdims(2) = idim2
  IF (ytype.EQ.'LON') THEN
    hname_dim(1) = 'lon'
    hname_dim(2) = 'lat'
    hunit1(1)    = 'degrees_east'
    hunit2(1)    = 'degrees_north'
  ELSE
    hname_dim(1) = 'xx'
    hname_dim(2) = 'yy'
    hunit1(1)    = 'meters'
    hunit2(1)    = 'meters'
  ENDIF
  ALLOCATE(px(kdim1))
  ALLOCATE(py(idim2))
  idim1 = kdim1
  ifull = kdim1*idim2

! case LWRITE_COORD = T
ELSE

  kdims(1) = kni
  hname_dim(1) = 'Number_of_points' 

  IF (ug%G%CGRID=="GAUSS     ") THEN
   
    hname_dim(2) = "latitude"
    hname_dim(3) = "longitude"
    hunit1(1)    = 'degrees_east'
    hunit2(1)    = 'degrees_north'     
    CALL get_gridtype_gauss(ug%XGRID_FULL_PAR,knlati=inlati)
    ALLOCATE(inlopa(inlati))
    CALL get_gridtype_gauss(ug%XGRID_FULL_PAR,knlopa=inlopa)
    kdims(2) = inlati
    kdims(3) = maxval(inlopa)
    DEALLOCATE(inlopa)

  ELSEIF (ug%G%CGRID/="LONLATVAL ".AND.(ug%G%CGRID/="IGN       ".OR.iimax*ijmax==il)) THEN
    
    IF (ytype.EQ.'LON') THEN
      hname_dim(2) = 'lon'
      hname_dim(3) = 'lat'
      hunit1(1)    = 'degrees_east'
      hunit2(1)    = 'degrees_north'      
    ELSE
      hname_dim(2) = 'xx'
      hname_dim(3) = 'yy'
      hunit1(1)    = 'meters'
      hunit2(1)    = 'meters'

    ENDIF
          
    IF (ug%G%CGRID=="CONF PROJ ") THEN
      CALL get_gridtype_conf_proj(ug%XGRID_FULL_PAR,kimax=iimax,kjmax=ijmax)
    ELSEIF(ug%G%CGRID=="CARTESIAN ") THEN
      CALL get_gridtype_cartesian(ug%XGRID_FULL_PAR,kimax=iimax,kjmax=ijmax)            
    ELSEIF(ug%G%CGRID=="LONLAT REG") THEN
      CALL get_gridtype_lonlat_reg(ug%XGRID_FULL_PAR,klon=iimax,klat=ijmax)
    ELSEIF(ug%G%CGRID=="LONLATROT ") THEN  
      CALL get_gridtype_lonlat_rot(ug%XGRID_FULL_PAR,klon=iimax,klat=ijmax)            
    ENDIF

    kdims(2) = iimax
    kdims(3) = ijmax

  ENDIF

  IF (lwrite_coord) THEN
    ALLOCATE(px(kni))
    ALLOCATE(py(kni))
  ENDIF
  idim1 = kni
  idim2 = kni
  ifull = kni

ENDIF
!
!
IF (lwrite_coord) THEN
  !
  IF (ASSOCIATED(ug%XGRID_FULL_PAR)) THEN
    CALL get_grid_coord(ug%G%CGRID, ug%G%NGRID_PAR, ug%G%XGRID_PAR, ksize_full, &
                        kluout,px=px,py=py,kl=kni,hgrid=ug%G%CGRID,pgrid_par=ug%XGRID_FULL_PAR)
  ELSEIF (ASSOCIATED(ug%G%XGRID_PAR)) THEN
    CALL get_grid_coord(ug%G%CGRID, ug%G%NGRID_PAR, ug%G%XGRID_PAR, ksize_full, &
                        kluout,px=px,py=py,kl=kni,hgrid=ug%G%CGRID,pgrid_par=ug%G%XGRID_PAR)
  ENDIF
  !
ELSEIF ( ug%G%CGRID.EQ.'CONF PROJ '.OR. ug%G%CGRID.EQ.'CARTESIAN '.OR. &
         ug%G%CGRID.EQ.'LONLAT REG' ) THEN
  !
  IF (ASSOCIATED(ug%XGRID_FULL_PAR)) THEN
    CALL get_grid_coord(ug%G%CGRID, ug%G%NGRID_PAR, ug%G%XGRID_PAR, ksize_full, &
                        kluout,px=zxx,py=zyy,kl=kni,hgrid=ug%G%CGRID,pgrid_par=ug%XGRID_FULL_PAR)
  ELSEIF (ASSOCIATED(ug%G%XGRID_PAR)) THEN
    CALL get_grid_coord(ug%G%CGRID, ug%G%NGRID_PAR, ug%G%XGRID_PAR, ksize_full, &
                        kluout,px=zxx,py=zyy,kl=kni,hgrid=ug%G%CGRID,pgrid_par=ug%G%XGRID_PAR)
  ENDIF
  !
  IF (ASSOCIATED(px)) THEN
    DO jj=1,SIZE(px)
      px(jj)=zxx(jj)
    ENDDO
  ENDIF
  IF (ASSOCIATED(py)) THEN
    DO jj=1,SIZE(py)
      py(jj)=zyy((jj-1)*(kni/SIZE(py))+1)
    ENDDO
  ENDIF
!
ELSEIF(ug%G%CGRID.EQ.'IGN       ')THEN
  !
  CALL get_ign_maskall(ug, kni, px, py)
  !
ENDIF
!
IF (PRESENT(plat) .AND. PRESENT(plon)) THEN
  !
  IF (ug%G%CGRID=="IGN       ".AND.iimax*ijmax/=il.AND.lwrite_coord) THEN
    !
    ALLOCATE(zxfull(ifull),zyfull(ifull))
    ALLOCATE(plat(ifull),plon(ifull))
    !
    CALL get_gridtype_ign(ug%XGRID_FULL_PAR,klambert=ilambert,px=zxfull,py=zyfull)
    CALL latlon_ign(ilambert,zxfull,zyfull,plat,plon)

    DEALLOCATE(zxfull, zyfull)
    !
  ELSEIF (ug%G%CGRID=="LONLATVAL ") THEN
    !
    ALLOCATE(plat(ifull),plon(ifull))
    CALL get_gridtype_lonlatval(ug%XGRID_FULL_PAR,px=plon,py=plat)
    !
  ELSEIF (ug%G%CGRID=="GAUSS     ") THEN
    !
    ALLOCATE(plat(ifull),plon(ifull))
    CALL get_gridtype_gauss(ug%XGRID_FULL_PAR,plat=plat,plon=plon) 
    !
  ELSE
    !
    ALLOCATE(plat(0),plon(0))
    !
  ENDIF
  !
ENDIF
!
!
IF ( inpatch/=0 ) THEN
  kdims(indims-1) = knpatch
  hname_dim(indims-1) = 'Number_of_Patches'
ENDIF
!
IF ( insnlayer/=0 ) THEN
  id = 1
  IF (inpatch/=0) id = 2
  kdims(indims-id) = insnlayer
  hname_dim(indims-id) = 'snow_layer'
ENDIF
!
!
IF (hprogram/='NOTIME ') THEN
  kdims(indims) = nf90_unlimited
  hname_dim(indims) = 'time'
ELSE
  kdims(indims) = 40
  hname_dim(indims) = 'char_len'
ENDIF
!
IF (lhook) CALL dr_hook('OL_DEFINE_DIM',1,zhook_handle)
!
END SUBROUTINE ol_define_dim