ip_rot_equid_cylind_egrid_mod.F90

Go to the documentation of this file.

 
 
module ip_rot_equid_cylind_egrid_mod
  use iso_fortran_env, only: real64
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use constants_mod, only: dpr, pi
  use earth_radius_mod
  implicit none
 
  private
  public :: ip_rot_equid_cylind_egrid
 
  integer, parameter :: kd = real64 
 
  type, extends(ip_grid) :: ip_rot_equid_cylind_egrid
     real(kd) :: rlon0
     real(kd) :: rlon1
     real(kd) :: rlat1
     real(kd) :: clat0
     real(kd) :: slat0
     real(kd) :: dlats
     real(kd) :: dlons
     real(kd) :: hi
     integer :: irot 
   contains
     procedure :: init_grib1
     procedure :: init_grib2
     procedure :: gdswzd => gdswzd_rot_equid_cylind_egrid
  end type ip_rot_equid_cylind_egrid
 
  INTEGER :: irot
 
  REAL(kind=kd) :: clat 
  REAL(kind=kd) :: clat0 
  REAL(kind=kd) :: clatr 
  REAL(kind=kd) :: clon 
  REAL(kind=kd) :: dlats 
  REAL(kind=kd) :: dlons 
  REAL(kind=kd) :: rerth 
  REAL(kind=kd) :: rlon0 
  REAL(kind=kd) :: slat 
  REAL(kind=kd) :: slat0 
  REAL(kind=kd) :: slatr 
 
contains
 
  subroutine init_grib1(self, g1_desc)
    class(ip_rot_equid_cylind_egrid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    integer :: iscan
    real(kd) :: rlat0
 
    real(kd) :: rlat1, rlon1, rlon0, slat1, clat1, slat0, clat0, clon1
    real(kd) :: slatr, clatr, clonr, rlatr, rlonr, dlats, dlons, hs, hi
    integer :: im, jm
 
    integer :: is1, kscan,  irot
 
    associate(kgds => g1_desc%gds)
      self%rerth = 6.3712e6_kd
      self%eccen_squared = 0.0
 
      im=kgds(2)
      jm=kgds(3)
 
      self%nscan_field_pos = 3
      self%nscan = mod(kgds(11)/32,2)
 
      rlat1=kgds(4)*1.e-3_kd
      rlon1=kgds(5)*1.e-3_kd
      rlat0=kgds(7)*1.e-3_kd
      rlon0=kgds(8)*1.e-3_kd
 
      irot=mod(kgds(6)/8,2)
      kscan=mod(kgds(11)/256,2)
      iscan=mod(kgds(11)/128,2)
      hi=(-1.)**iscan
      slat1=sin(rlat1/dpr)
      clat1=cos(rlat1/dpr)
      slat0=sin(rlat0/dpr)
      clat0=cos(rlat0/dpr)
      hs=sign(1._kd,mod(rlon1-rlon0+180+3600,360._kd)-180)
      clon1=cos((rlon1-rlon0)/dpr)
      slatr=clat0*slat1-slat0*clat1*clon1
      clatr=sqrt(1-slatr**2)
      clonr=(clat0*clat1*clon1+slat0*slat1)/clatr
      rlatr=dpr*asin(slatr)
      rlonr=hs*dpr*acos(clonr)
      dlats=rlatr/(-(jm-1)/2)
      dlons=rlonr/(-((im * 2 - 1) -1)/2)
 
      IF(kscan.EQ.0) THEN
         is1=(jm+1)/2
      ELSE
         is1=jm/2
      ENDIF
 
      self%im = im
      self%jm = jm
      self%rlon0 = rlon0
      self%rlon1 = rlon1
      self%rlat1 = rlat1
      self%clat0 = clat0
      self%slat0 = slat0
      self%dlats = dlats
      self%dlons = dlons
      self%hi = hi
      self%irot = irot
      self%kscan = kscan
 
 
    end associate
 
  end subroutine init_grib1
 
 
  subroutine init_grib2(self, g2_desc)
    class(ip_rot_equid_cylind_egrid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    integer :: iscale, iscan
    real(kd) :: rlat0
    integer :: i_offset_odd!, i_offset_even
 
    associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
      CALL earth_radius(igdtmpl,igdtlen,self%rerth,self%eccen_squared)
 
      ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      ! ROUTINE ONLY WORKS FOR "E"-STAGGER GRIDS.
      !   "V" GRID WHEN BIT 5 IS '1' AND BIT 6 IS '0'.
      !   "H" GRID WHEN BIT 5 IS '0' AND BIT 6 IS '1'.
      ! I_OFFSET_ODD=MOD(IGDTMPL(19)/8,2)
      ! I_OFFSET_EVEN=MOD(IGDTMPL(19)/4,2)
      ! IF(I_OFFSET_ODD==I_OFFSET_EVEN) THEN
      !    CALL ROT_EQUID_CYLIND_EGRID_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
      !    RETURN
      ! ENDIF
      ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 
      self%IM=igdtmpl(8)
      self%JM=igdtmpl(9)
 
      self%NSCAN=mod(igdtmpl(16)/32,2)
      self%nscan_field_pos = 3
 
      iscale=igdtmpl(10)*igdtmpl(11)
      IF(iscale==0) iscale=10**6
 
      self%RLON0=float(igdtmpl(21))/float(iscale)
      self%DLATS=float(igdtmpl(18))/float(iscale)
      ! THE GRIB2 CONVENTION FOR "I" RESOLUTION IS TWICE WHAT THIS ROUTINE ASSUMES.
      self%DLONS=float(igdtmpl(17))/float(iscale) * 0.5_kd
 
      self%IROT=mod(igdtmpl(14)/8,2)
 
      i_offset_odd=mod(igdtmpl(19)/8,2)
      self%KSCAN=i_offset_odd
      iscan=mod(igdtmpl(19)/128,2)
 
      self%HI=(-1.)**iscan
 
      rlat0=float(igdtmpl(20))/float(iscale)
      rlat0=rlat0+90.0_kd
 
      self%SLAT0=sin(rlat0/dpr)
      self%CLAT0=cos(rlat0/dpr)
 
      self%RLAT1=float(igdtmpl(12))/float(iscale)
      self%RLON1=float(igdtmpl(13))/float(iscale)
    end associate
  end subroutine init_grib2
 
 
  SUBROUTINE gdswzd_rot_equid_cylind_egrid(self,IOPT,NPTS,&
       FILL,XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    IMPLICIT NONE
    !
    class(ip_rot_equid_cylind_egrid), intent(in) :: self
 
    INTEGER,         INTENT(IN   ) :: IOPT, NPTS
    INTEGER,         INTENT(  OUT) :: NRET
    !
    REAL,            INTENT(IN   ) :: FILL
    REAL,            INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
    REAL,            INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
    REAL, OPTIONAL,  INTENT(  OUT) :: CROT(NPTS),SROT(NPTS)
    REAL, OPTIONAL,  INTENT(  OUT) :: XLON(NPTS),XLAT(NPTS)
    REAL, OPTIONAL,  INTENT(  OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
    !
    INTEGER                        :: IM, JM, IS1, N
    INTEGER                        :: KSCAN
    !    INTEGER                        :: I_OFFSET_ODD, I_OFFSET_EVEN
    !
    LOGICAL                        :: LROT, LMAP, LAREA
    !
    REAL(KIND=kd)                  :: rlat1, rlon1
    REAL(KIND=kd)                  :: clonr
    REAL(KIND=kd)                  :: rlatr, rlonr, sbd, wbd, hs, hi
    REAL                           :: XMAX, XMIN, YMAX, YMIN, XPTF, YPTF
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    IF(PRESENT(crot)) crot=fill
    IF(PRESENT(srot)) srot=fill
    IF(PRESENT(xlon)) xlon=fill
    IF(PRESENT(xlat)) xlat=fill
    IF(PRESENT(ylon)) ylon=fill
    IF(PRESENT(ylat)) ylat=fill
    IF(PRESENT(area)) area=fill
 
    rlon0=self%rlon0
    irot=self%irot
    im=self%im * 2 - 1
    jm=self%jm
    dlats=self%dlats
    dlons=self%dlons
    kscan=self%kscan
    hi=self%hi
    slat0=self%slat0
    clat0=self%clat0
    rlat1=self%rlat1
    rlon1=self%rlon1
 
    rerth = self%rerth
 
    ! IS THE EARTH RADIUS DEFINED?
    IF(rerth<0.)THEN
       CALL rot_equid_cylind_egrid_error(iopt,fill,rlat,rlon,xpts,ypts,npts)
       RETURN
    ENDIF
 
    sbd=rlat1
    wbd=rlon1
 
    IF (wbd > 180.0) wbd = wbd - 360.0
    IF(kscan.EQ.0) THEN
       is1=(jm+1)/2
    ELSE
       is1=jm/2
    ENDIF
 
    xmin=0
    xmax=im+2
    ymin=0
    ymax=jm+1
    nret=0
 
    IF(PRESENT(crot).AND.PRESENT(srot))THEN
       lrot=.true.
    ELSE
       lrot=.false.
    ENDIF
    IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
       lmap=.true.
    ELSE
       lmap=.false.
    ENDIF
    IF(PRESENT(area))THEN
       larea=.true.
    ELSE
       larea=.false.
    ENDIF
 
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
    IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
       DO n=1,npts
          xptf=ypts(n)+(xpts(n)-is1)
          yptf=ypts(n)-(xpts(n)-is1)+kscan
          IF(xptf.GE.xmin.AND.xptf.LE.xmax.AND. &
               yptf.GE.ymin.AND.yptf.LE.ymax) THEN
             hs=hi*sign(1.,xptf-(im+1)/2)
             select type(desc => self%descriptor)
             type is(grib1_descriptor)
                rlonr=(xptf-(im+1)/2)*dlons
                rlatr=(yptf-(jm+1)/2)*dlats
             type is(grib2_descriptor)
                rlonr=(xptf-1.0_kd)*dlons + wbd
                rlatr=(yptf-1.0_kd)*dlats + sbd 
             end select
             clonr=cos(rlonr/dpr)
             slatr=sin(rlatr/dpr)
             clatr=cos(rlatr/dpr)
             slat=clat0*slatr+slat0*clatr*clonr
             IF(slat.LE.-1) THEN
                clat=0.
                clon=cos(rlon0/dpr)
                rlon(n)=0
                rlat(n)=-90
             ELSEIF(slat.GE.1) THEN
                clat=0.
                clon=cos(rlon0/dpr)
                rlon(n)=0
                rlat(n)=90
             ELSE
                clat=sqrt(1-slat**2)
                clon=(clat0*clatr*clonr-slat0*slatr)/clat
                clon=min(max(clon,-1._kd),1._kd)
                rlon(n)=real(mod(rlon0+hs*dpr*acos(clon)+3600,360._kd))
                rlat(n)=real(dpr*asin(slat))
             ENDIF
             nret=nret+1
             IF(lrot) CALL rot_equid_cylind_egrid_vect_rot(rlon(n),crot(n),srot(n))
             IF(lmap) CALL rot_equid_cylind_egrid_map_jacob(fill, rlon(n), &
                  xlon(n),xlat(n),ylon(n),ylat(n))
             IF(larea) CALL rot_equid_cylind_egrid_grid_area(fill, area(n))
          ELSE
             rlon(n)=fill
             rlat(n)=fill
          ENDIF
       ENDDO
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
    ELSEIF(iopt.EQ.-1) THEN
       DO n=1,npts
          IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
             hs=sign(1._kd,mod(rlon(n)-rlon0+180+3600,360._kd)-180)
             clon=cos((rlon(n)-rlon0)/dpr)
             slat=sin(rlat(n)/dpr)
             clat=cos(rlat(n)/dpr)
             slatr=clat0*slat-slat0*clat*clon
             IF(slatr.LE.-1) THEN
                clatr=0.
                rlonr=0
                rlatr=-90
             ELSEIF(slatr.GE.1) THEN
                clatr=0.
                rlonr=0
                rlatr=90
             ELSE
                clatr=sqrt(1-slatr**2)
                clonr=(clat0*clat*clon+slat0*slat)/clatr
                clonr=min(max(clonr,-1._kd),1._kd)
                rlonr=hs*dpr*acos(clonr)
                rlatr=dpr*asin(slatr)
             ENDIF
             select type(desc => self%descriptor)
             type is(grib1_descriptor)
                xptf=real(((rlonr-wbd)/dlons)+1.0_kd)
                yptf=real(((rlatr-sbd)/dlats)+1.0_kd)
             type is(grib2_descriptor)
                xptf=real((im+1)/2+rlonr/dlons)
                yptf=real((jm+1)/2+rlatr/dlats)
             end select
 
             IF(xptf.GE.xmin.AND.xptf.LE.xmax.AND. &
                  yptf.GE.ymin.AND.yptf.LE.ymax) THEN
                xpts(n)=is1+(xptf-(yptf-kscan))/2
                ypts(n)=(xptf+(yptf-kscan))/2
                nret=nret+1
                IF(lrot) CALL rot_equid_cylind_egrid_vect_rot(rlon(n),crot(n),srot(n))
                IF(lmap) CALL rot_equid_cylind_egrid_map_jacob(fill, rlon(n), &
                     xlon(n),xlat(n),ylon(n),ylat(n))
                IF(larea) CALL rot_equid_cylind_egrid_grid_area(fill, area(n))
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ELSE
             xpts(n)=fill
             ypts(n)=fill
          ENDIF
       ENDDO
    ENDIF
  END SUBROUTINE gdswzd_rot_equid_cylind_egrid
 
  SUBROUTINE rot_equid_cylind_egrid_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
    IMPLICIT NONE
    !
    INTEGER, INTENT(IN   ) :: IOPT, NPTS
    !
    REAL,    INTENT(IN   ) :: FILL
    REAL,    INTENT(  OUT) :: RLAT(NPTS),RLON(NPTS)
    REAL,    INTENT(  OUT) :: XPTS(NPTS),YPTS(NPTS)
 
    IF(iopt>=0) THEN
       rlon=fill
       rlat=fill
    ENDIF
    IF(iopt<=0) THEN
       xpts=fill
       ypts=fill
    ENDIF
  END SUBROUTINE rot_equid_cylind_egrid_error
 
  SUBROUTINE rot_equid_cylind_egrid_vect_rot(RLON, CROT, SROT)
    IMPLICIT NONE
 
    REAL         ,    INTENT(IN   ) :: RLON
    REAL         ,    INTENT(  OUT) :: CROT, SROT
 
    REAL(KIND=kd)                   :: slon
 
    IF(irot.EQ.1) THEN
       IF(clatr.LE.0) THEN
          crot=real(-sign(1._kd,slatr*slat0))
          srot=0.
       ELSE
          slon=sin((rlon-rlon0)/dpr)
          crot=real((clat0*clat+slat0*slat*clon)/clatr)
          srot=real(slat0*slon/clatr)
       ENDIF
    ELSE
       crot=1.
       srot=0.
    ENDIF
 
  END SUBROUTINE rot_equid_cylind_egrid_vect_rot
 
  SUBROUTINE rot_equid_cylind_egrid_map_jacob(FILL, RLON, &
       XLON, XLAT, YLON, YLAT)
    IMPLICIT NONE
 
    REAL         ,    INTENT(IN   ) :: FILL, RLON
    REAL         ,    INTENT(  OUT) :: XLON, XLAT, YLON, YLAT
 
    REAL(KIND=kd)                   :: slon, term1, term2
    REAL(KIND=kd)                   :: xlatf, xlonf, ylatf, ylonf
 
    IF(clatr.LE.0._kd) THEN
       xlon=fill
       xlat=fill
       ylon=fill
       ylat=fill
    ELSE
       slon=sin((rlon-rlon0)/dpr)
       term1=(clat0*clat+slat0*slat*clon)/clatr
       term2=slat0*slon/clatr
       xlonf=term1*clat/(dlons*clatr)
       xlatf=-term2/(dlons*clatr)
       ylonf=term2*clat/dlats
       ylatf=term1/dlats
       xlon=real(xlonf-ylonf)
       xlat=real(xlatf-ylatf)
       ylon=real(xlonf+ylonf)
       ylat=real(xlatf+ylatf)
    ENDIF
 
  END SUBROUTINE rot_equid_cylind_egrid_map_jacob
 
  SUBROUTINE rot_equid_cylind_egrid_grid_area(FILL, AREA)
    IMPLICIT NONE
 
    REAL,             INTENT(IN   ) :: FILL
    REAL,             INTENT(  OUT) :: AREA
 
    IF(clatr.LE.0._kd) THEN
       area=fill
    ELSE
       area=real(rerth**2*clatr*dlats*dlons)*2/dpr**2
    ENDIF
 
  END SUBROUTINE rot_equid_cylind_egrid_grid_area
 
end module ip_rot_equid_cylind_egrid_mod