regular_grid_spawn.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 regular_grid_spawn(KLUOUT,                               &
                                      kl1, kimax1,kjmax1,px1,py1,pdx1,pdy1, &
                                      kxor, kyor, kdxratio, kdyratio,       &
                                      kxsize, kysize,                       &
                                      kl2, kimax2,kjmax2,px2,py2,pdx2,pdy2  )  
!     ################################################################
!
!!****  *REGULAR_GRID_SPAWN* - routine to read in namelist the horizontal grid
!!
!!    PURPOSE
!!    -------
!!
!!**  METHOD
!!    ------
!!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!    REFERENCE
!!    ---------
!!
!!
!!    AUTHOR
!!    ------
!!      V. Masson   *Meteo France*
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    01/2004 
!-------------------------------------------------------------------------------
!
!*       0.    DECLARATIONS
!              ------------
!
USE modd_surf_par, ONLY : nundef
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
USE modi_abor1_sfx
!
IMPLICIT NONE
!
!*       0.1   Declarations of arguments
!              -------------------------
!
INTEGER,                      INTENT(IN)    :: kluout     ! output listing logical unit
INTEGER,                      INTENT(IN)    :: kl1        ! total number of points KIMAX1 * KJMAX1
INTEGER,                      INTENT(IN)    :: kimax1     ! number of points in x direction
INTEGER,                      INTENT(IN)    :: kjmax1     ! number of points in y direction
REAL, DIMENSION(KL1),         INTENT(IN)    :: px1        ! X coordinate of all points
REAL, DIMENSION(KL1),         INTENT(IN)    :: py1        ! Y coordinate of all points
REAL, DIMENSION(KL1),         INTENT(IN)    :: pdx1       ! X mesh size of all points
REAL, DIMENSION(KL1),         INTENT(IN)    :: pdy1       ! Y mesh size of all points
INTEGER,                      INTENT(IN)    :: kxor       ! position of modified bottom left point
INTEGER,                      INTENT(IN)    :: kyor       ! according to initial grid
INTEGER,                      INTENT(IN)    :: kxsize     ! number of grid meshes in initial grid to be
INTEGER,                      INTENT(IN)    :: kysize     ! covered by the modified grid
INTEGER,                      INTENT(IN)    :: kdxratio   ! resolution ratio between modified grid
INTEGER,                      INTENT(IN)    :: kdyratio   ! and initial grid
INTEGER,                      INTENT(IN)    :: kl2        ! total number of points KIMAX2 * KJMAX2
INTEGER,                      INTENT(IN)    :: kimax2     ! number of points in x direction
INTEGER,                      INTENT(IN)    :: kjmax2     ! number of points in y direction
REAL, DIMENSION(KL2),         INTENT(OUT)   :: px2        ! X coordinate of all points
REAL, DIMENSION(KL2),         INTENT(OUT)   :: py2        ! Y coordinate of all points
REAL, DIMENSION(KL2),         INTENT(OUT)   :: pdx2       ! X mesh size of all points
REAL, DIMENSION(KL2),         INTENT(OUT)   :: pdy2       ! Y mesh size of all points
!
!*       0.2   Declarations of local variables
!              -------------------------------
!
!* initial grid
!
REAL, DIMENSION(:),   ALLOCATABLE :: zxm1     ! X coordinate of center of mesh (IIMAX1   points)
REAL, DIMENSION(:),   ALLOCATABLE :: zym1     ! Y coordinate of center of mesh (IJMAX1   points)
REAL, DIMENSION(:),   ALLOCATABLE :: zxhat1   ! X coordinate of left side      (IIMAX1+1 points)
REAL, DIMENSION(:),   ALLOCATABLE :: zyhat1   ! Y coordinate of bottom side    (IJMAX1+1 points)
!
!* new grid
!
REAL, DIMENSION(:),   ALLOCATABLE :: zxhat2   ! X coordinate of left side      (IIMAX2 points)
REAL, DIMENSION(:),   ALLOCATABLE :: zyhat2   ! Y coordinate of bottom side    (IJMAX2 points)
!
!* other variables
!
INTEGER     :: jl            ! loop counter
INTEGER     :: ji,jj         ! loop controls relatively to modified grid
INTEGER     :: jibox,jjbox   ! grid mesh relatively to initial grid
REAL        :: zcoef         ! ponderation coefficient for linear interpolation
REAL(KIND=JPRB) :: zhook_handle
!
!------------------------------------------------------------------------------
!
!*       1.    Coherence tests
!              ---------------
!
!* tests
!
IF (lhook) CALL dr_hook('REGULAR_GRID_SPAWN',0,zhook_handle)
IF ( kxor+kxsize-1 > kimax1 ) THEN
  WRITE(kluout,*) 'spawned domain is not contained in the input domain'
  WRITE(kluout,*) 'IXOR = ', kxor, ' IXSIZE = ', kxsize,&
                    ' with NIMAX(file) = ', kimax1  
  CALL abor1_sfx('REGULAR_GRID_SPAWN: (1) SPAWNED DOMAIN NOT CONTAINED IN INPUT DOMAIN')
END IF
IF ( kyor+kysize-1 > kjmax1 ) THEN
  WRITE(kluout,*) 'spawned domain is not contained in the input domain'
  WRITE(kluout,*) 'IYOR = ', kyor, ' IYSIZE = ', kysize,&
                    ' with NJMAX(file) = ', kjmax1  
  CALL abor1_sfx('REGULAR_GRID_SPAWN: (2) SPAWNED DOMAIN NOT CONTAINED IN INPUT DOMAIN')
END IF
!
!------------------------------------------------------------------------------
!
!*       2.    Center of mesh coordinate arrays for each direction separately
!              --------------------------------------------------------------
!
ALLOCATE(zxm1(kimax1))
ALLOCATE(zym1(kjmax1))
ALLOCATE(zxhat1(kimax1+1))
ALLOCATE(zyhat1(kjmax1+1))
ALLOCATE(zxhat2(kimax2+1))
ALLOCATE(zyhat2(kjmax2+1))
!
zxm1(:) = px1(1:kimax1)
DO jl=1,kl1
  IF (mod(jl,kimax1)==0) zym1(jl/kimax1) = py1(jl)
END DO
!
!------------------------------------------------------------------------------
!
!*       3.    side of mesh coordinate arrays for each direction separately
!              ------------------------------------------------------------
!
!
IF (kimax1==1) THEN
  zxhat1(1) = zxm1(1) - 0.5 * pdx1(1)
  zxhat1(2) = zxm1(1) + 0.5 * pdx1(1)
ELSE
  zxhat1(1) = 1.5 * zxm1(1) - 0.5 * zxm1(2)
  DO ji=2,kimax1
    zxhat1(ji) = 0.5 * zxm1(ji-1) + 0.5 * zxm1(ji)
  END DO
  zxhat1(kimax1+1) = 1.5 * zxm1(kimax1) - 0.5 * zxm1(kimax1-1)
END IF
!
IF (kjmax1==1) THEN
  zyhat1(1) = zym1(1) - 0.5 * pdy1(1)
  zyhat1(2) = zym1(1) + 0.5 * pdy1(1)
ELSE
  zyhat1(1) = 1.5 * zym1(1) - 0.5 * zym1(2)
  DO jj=2,kjmax1
    zyhat1(jj) = 0.5 * zym1(jj-1) + 0.5 * zym1(jj)
  END DO
  zyhat1(kjmax1+1) = 1.5 * zym1(kjmax1) - 0.5 * zym1(kjmax1-1)
END IF
!
!------------------------------------------------------------------------------
!
!*       5.    Interpolation of coordinate arrays for each direction separately
!              ----------------------------------------------------------------
!
!* X coordinate array
!
DO ji=1,kimax2
  jibox=(ji-1)/kdxratio + kxor
  zcoef= float(mod(ji-1,kdxratio))/float(kdxratio)
  zxhat2(ji)=(1.-zcoef)*zxhat1(jibox)+zcoef*zxhat1(jibox+1)
END DO
IF (kimax2==1) THEN
  zxhat2(kimax2+1) = zxhat2(kimax2) + zxhat1(jibox+1) - zxhat1(jibox)
ELSE
  zxhat2(kimax2+1) = 2. * zxhat2(kimax2) - zxhat2(kimax2-1)
END IF
!
!
!* Y coordinate array
!
DO jj=1,kjmax2
  jjbox=(jj-1)/kdyratio + kyor
  zcoef= float(mod(jj-1,kdyratio))/float(kdyratio)
  zyhat2(jj)=(1.-zcoef)*zyhat1(jjbox)+zcoef*zyhat1(jjbox+1)
END DO
IF (kjmax2==1) THEN
  zyhat2(kjmax2+1) = zyhat2(kjmax2) + zyhat1(jjbox+1) - zyhat1(jjbox)
ELSE
  zyhat2(kjmax2+1) = 2. * zyhat2(kjmax2) - zyhat2(kjmax2-1)
END IF
!---------------------------------------------------------------------------
DEALLOCATE(zxm1)
DEALLOCATE(zym1)
DEALLOCATE(zxhat1)
DEALLOCATE(zyhat1)
!------------------------------------------------------------------------------
!
!*       5.    Coordinate arrays of all points
!              -------------------------------
!
DO jj=1,kjmax2
  DO ji=1,kimax2
    jl = (jj-1) * kimax2 + ji
      px2(jl) = 0.5 * zxhat2(ji) + 0.5 * zxhat2(ji+1)
      pdx2(jl) = zxhat2(ji+1) - zxhat2(ji)
      py2(jl) = 0.5 * zyhat2(jj) + 0.5 * zyhat2(jj+1)
      pdy2(jl) = zyhat2(jj+1) - zyhat2(jj)
  END DO
END DO
!
!---------------------------------------------------------------------------
DEALLOCATE(zxhat2)
DEALLOCATE(zyhat2)
IF (lhook) CALL dr_hook('REGULAR_GRID_SPAWN',1,zhook_handle)
!---------------------------------------------------------------------------
!
END SUBROUTINE regular_grid_spawn