latlonmask_conf_proj.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 latlonmask_conf_proj(KGRID_PAR,PGRID_PAR,OLATLONMASK)
!     ##################################
!
!!**** *LATLONMASK* builds the latiude and longitude mask including the grid
!!
!!    PURPOSE
!!    -------
!!
!!    METHOD
!!    ------
!!
!!    Two tests are performed:
!!
!!   1) test if the points of the mask are in the domain
!!
!!   2) fills the mask points corresponding to points scanning
!!   
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!
!!    REFERENCE
!!    ---------
!!
!!    AUTHOR
!!    ------
!!
!!      V. Masson       Meteo-France
!!
!!    MODIFICATION
!!    ------------
!!      
!!      Original        19/07/95
!----------------------------------------------------------------------------
!
!*    0.     DECLARATION
!            -----------
!
USE mode_gridtype_conf_proj
!
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*    0.1    Declaration of arguments
!            ------------------------
!
INTEGER,                       INTENT(IN)  :: kgrid_par   ! size of PGRID_PAR
REAL,    DIMENSION(KGRID_PAR), INTENT(IN)  :: pgrid_par   ! parameters defining this grid
LOGICAL, DIMENSION(720,360),   INTENT(OUT) :: olatlonmask ! mask where data are to be read
!
!*    0.2    Declaration of local variables
!            ------------------------------
!
INTEGER                           :: iimax    ! x coordinates array dimension
INTEGER                           :: ijmax    ! y coordinates array dimension
INTEGER                           :: ji, jj   ! loop counters
INTEGER                           :: jlon,jlat! loop counters
REAL                              :: zxmin    ! minimum of X for domain
REAL                              :: zxmax    ! maximum of X for domain
REAL                              :: zymin    ! minimum of Y for domain
REAL                              :: zymax    ! maximum of Y for domain
REAL, DIMENSION(:),   ALLOCATABLE :: zdx      ! X grid mesh size
REAL, DIMENSION(:),   ALLOCATABLE :: zdy      ! Y grid mesh size
REAL, DIMENSION(:),   ALLOCATABLE :: zx       ! X conformal coordinate of center of grid mesh
REAL, DIMENSION(:),   ALLOCATABLE :: zy       ! Y conformal coordinate of center of grid mesh
REAL, DIMENSION(:),   ALLOCATABLE :: zxcorner ! X conformal coordinate of corner of grid mesh
REAL, DIMENSION(:),   ALLOCATABLE :: zycorner ! Y conformal coordinate of corner of grid mesh
REAL, DIMENSION(:),   ALLOCATABLE :: zlon     ! corner points longitudes
REAL, DIMENSION(:),   ALLOCATABLE :: zlat     ! corner points latitudes
REAL, DIMENSION(:,:), ALLOCATABLE :: zlon2d   ! corner points longitudes
REAL, DIMENSION(:,:), ALLOCATABLE :: zlat2d   ! corner points latitudes
REAL, DIMENSION(720,360)          :: zx_mask  ! mask points X value
REAL, DIMENSION(720,360)          :: zy_mask  ! mask points Y value
INTEGER, DIMENSION(720,360)       :: icount1  ! counter
INTEGER, DIMENSION(720,360)       :: icount2  ! counter
REAL, DIMENSION(720)              :: zlon_mask! mask points longitudes
REAL, DIMENSION(720)              :: zlat_mask! mask points latitudes
REAL                              :: zlat0    ! reference latitude
REAL                              :: zlon0    ! reference longitude
REAL                              :: zrpk     ! projection parameter 
!                                             !   K=1 : stereographic north pole
!                                             ! 0<K<1 : Lambert, north hemisphere
!                                             !   K=0 : Mercator
!                                             !-1<K<0 : Lambert, south hemisphere
!                                             !   K=-1: stereographic south pole
REAL                              :: zbeta    ! angle between grid and reference longitude
REAL                              :: zlator   ! latitude  of point of coordinates X=0, Y=0
REAL                              :: zlonor   ! longitude of point of coordinates X=0, Y=0
!
INTEGER                           :: iverb=1  ! verbosity level
REAL(KIND=JPRB) :: zhook_handle
!----------------------------------------------------------------------------
!
IF (lhook) CALL dr_hook('LATLONMASK_CONF_PROJ',0,zhook_handle)
 CALL get_gridtype_conf_proj(pgrid_par,kimax=iimax,kjmax=ijmax)
!
!
ALLOCATE(zx(iimax*ijmax))
ALLOCATE(zy(iimax*ijmax))
ALLOCATE(zdx(iimax*ijmax))
ALLOCATE(zdy(iimax*ijmax))
ALLOCATE(zxcorner((iimax+1)*(ijmax+1)))
ALLOCATE(zycorner((iimax+1)*(ijmax+1)))
ALLOCATE(zlon((iimax+1)*(ijmax+1)))
ALLOCATE(zlat((iimax+1)*(ijmax+1)))
!
!-------------------------------------------------------------------------------
!
olatlonmask(:,:) = .false.
!
!-------------------------------------------------------------------------------
!
!*      1.   Limits of the domain conformal plane coordinates
!            ------------------------------------------------
!
 CALL get_gridtype_conf_proj(pgrid_par,px=zx,py=zy,pdx=zdx,pdy=zdy)
!
zxmin =      minval(zx) - 3.* maxval(zdx)/2.
zxmax =      maxval(zx) + 3.* maxval(zdx)/2.
!
zymin =      minval(zy) - 3.* maxval(zdy)/2.
zymax =      maxval(zy) + 3.* maxval(zdy)/2.
!
DO jj=1,ijmax
  DO ji=1,iimax
    zxcorner(ji+(jj-1)*(iimax+1)) = zx(ji+(jj-1)*iimax) - 0.5*zdx(ji+(jj-1)*iimax)
    zycorner(ji+(jj-1)*(iimax+1)) = zy(ji+(jj-1)*iimax) - 0.5*zdy(ji+(jj-1)*iimax)
  END DO
  zxcorner(iimax+1+(jj-1)*(iimax+1)) = zx(iimax+(jj-1)*iimax) + 0.5*zdx(iimax+(jj-1)*iimax)
  zycorner(iimax+1+(jj-1)*(iimax+1)) = zy(iimax+(jj-1)*iimax) - 0.5*zdx(iimax+(jj-1)*iimax)
END DO
DO ji=1,iimax
  zxcorner(ji+ijmax*(iimax+1)) = zx(ji+(ijmax-1)*iimax) - 0.5*zdx(ji+(ijmax-1)*iimax)
  zycorner(ji+ijmax*(iimax+1)) = zy(ji+(ijmax-1)*iimax) + 0.5*zdy(ji+(ijmax-1)*iimax)
END DO
zxcorner((ijmax+1)*(iimax+1)) = zx(ijmax*iimax) + 0.5*zdx(ijmax*iimax)
zycorner((ijmax+1)*(iimax+1)) = zy(ijmax*iimax) + 0.5*zdy(ijmax*iimax)
!
DEALLOCATE(zdx)
DEALLOCATE(zdy)
DEALLOCATE(zx)
DEALLOCATE(zy)
!
!-------------------------------------------------------------------------------
!
!*      2.   Definition of the coordinates at center of the mask meshes
!            ----------------------------------------------------------
!
zlon_mask(:)= (/ (  jlon     /2. - 0.25 , jlon=1,720 ) /)
!
!*      3.   Longitude correction / LON0
!            ---------------------------
!
 CALL get_gridtype_conf_proj(pgrid_par,                   &
                              plat0=zlat0,plon0=zlon0,prpk=zrpk,        &
                              pbeta=zbeta,plator=zlator,plonor=zlonor   )  
!
zlon_mask(:)=zlon_mask(:)+nint((zlon0-zlon_mask(:))/360.)*360.
!
!*      4.   X and Y of the points of the mask
!            ---------------------------------
!
DO jlat=1,SIZE(olatlonmask,2)
  zlat_mask(:) = (jlat-180)/2. - 0.25
  CALL xy_conf_proj(zlat0,zlon0,zrpk,zbeta,zlator,zlonor, &
                      zx_mask(:,jlat),zy_mask(:,jlat),    &
                      zlat_mask(:),zlon_mask(:)   )
END DO
!
!*      5.   Are the points in the domain?
!            ----------------------------
!
WHERE (        zx_mask(:,:) >= zxmin .AND. zx_mask(:,:) <= zxmax   &
           .AND. zy_mask(:,:) >= zymin .AND. zy_mask(:,:) <= zymax )  
  olatlonmask(:,:) = .true.
END WHERE
!
!-------------------------------------------------------------------------------
!
!*      6.   Latitude and longitude of the points of the domain
!            --------------------------------------------------
!
 CALL latlon_conf_proj(zlat0,zlon0,zrpk,zbeta,zlator,zlonor, &
                      zxcorner,zycorner,zlat,zlon           )  
!
!*      7.   Longitudes between 0. and 360.
!            ------------------------------
!
zlon(:) = zlon(:) + nint((180.-zlon(:))/360.)*360.
!
!*      8.   Loop on grid points
!            --------------------
!
icount1(:,:) = 0
icount2(:,:) = 0
!
ALLOCATE(zlat2d(iimax+1,ijmax+1))
ALLOCATE(zlon2d(iimax+1,ijmax+1))
!
DO jj=1,ijmax+1
  DO ji=1,iimax+1
    zlat2d(ji,jj) = zlat(ji+(jj-1)*(iimax+1))
    zlon2d(ji,jj) = zlon(ji+(jj-1)*(iimax+1))
  END DO
END DO
!
DO jj=1,ijmax+1
  DO ji=1,iimax+1
    !
    !*      8.1  localisation of the point
    !            -------------------------
    !
    jlat = min( 1 + int( ( zlat2d(ji,jj) + 90. ) * 2. ) ,360)
    jlon = min( 1 + int( ( zlon2d(ji,jj)       ) * 2. ) ,720)
    !
    icount1(jlon,jlat) = icount1(jlon,jlat) + 1
    !
    !*      8.2  Does point contain data?
    !            ------------------------
    !
    !
    !*      8.3  point contains data
    !            -------------------
    !
    icount2(jlon,jlat) = icount2(jlon,jlat) + 1
    !
    !*      8.4  Boundary effects
    !            ----------------
    !
    jlat = min( 1 + int( ( zlat2d(min(ji+1,iimax),jj) + 90. ) * 2. ) ,360)
    jlon = min( 1 + int( ( zlon2d(min(ji+1,iimax),jj)       ) * 2. ) ,720)
    icount1(jlon,jlat) = icount1(jlon,jlat) + 1
    icount2(jlon,jlat) = icount2(jlon,jlat) + 1
    !
    jlat = min( 1 + int( ( zlat2d(ji,min(jj+1,ijmax)) + 90. ) * 2. ) ,360)
    jlon = min( 1 + int( ( zlon2d(ji,min(jj+1,ijmax))       ) * 2. ) ,720)
    icount1(jlon,jlat) = icount1(jlon,jlat) + 1
    icount2(jlon,jlat) = icount2(jlon,jlat) + 1
    !
    jlat = min( 1 + int( ( zlat2d(max(ji-1,1),jj) + 90. ) * 2. ) ,360)
    jlon = min( 1 + int( ( zlon2d(max(ji-1,1),jj)       ) * 2. ) ,720)
    icount1(jlon,jlat) = icount1(jlon,jlat) + 1
    icount2(jlon,jlat) = icount2(jlon,jlat) + 1
    !
    jlat = min( 1 + int( ( zlat2d(ji,max(jj-1,1)) + 90. ) * 2. ) ,360)
    jlon = min( 1 + int( ( zlon2d(ji,max(jj-1,1))       ) * 2. ) ,720)
    icount1(jlon,jlat) = icount1(jlon,jlat) + 1
    icount2(jlon,jlat) = icount2(jlon,jlat) + 1
    !
  END DO
END DO
!
!*      9.   Surface type check (if points are present in mask mesh)
!            ------------------
!
WHERE (icount1(:,:) > 0 .AND. icount2(:,:) == 0)
  olatlonmask(:,:) = .false.
END WHERE
!
WHERE (icount1(:,:) > 0 .AND. icount2(:,:) > 0)
  olatlonmask(:,:) = .true.
END WHERE
!
zlat_mask(1:360)= (/ ( (jlat-180)/2. - 0.25 , jlat=1,360 ) /) 
!
DO jlon=1,720
  DO jlat=1,360
    IF ( (icount1(jlon,jlat) > 0 .OR. olatlonmask(jlon,jlat)) .AND.  iverb > 1) &
      WRITE(*,'(2(I3,1X),2(F6.2,1X),2(F8.0,1X),L1)') jlon,jlat,zlon_mask(jlon),zlat_mask(jlat), &
                zx_mask(jlon,jlat),zy_mask(jlon,jlat),olatlonmask(jlon,jlat)  

  END DO
END DO
!
!-------------------------------------------------------------------------------
!
DEALLOCATE(zlon  )
DEALLOCATE(zlat  )
DEALLOCATE(zlon2d)
DEALLOCATE(zlat2d)
DEALLOCATE(zxcorner)
DEALLOCATE(zycorner)
IF (lhook) CALL dr_hook('LATLONMASK_CONF_PROJ',1,zhook_handle)
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE latlonmask_conf_proj