ip_equid_cylind_grid_mod.F90

Go to the documentation of this file.

 
 
module ip_equid_cylind_grid_mod
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use earth_radius_mod
  implicit none
 
  private
  public :: ip_equid_cylind_grid
 
  type, extends(ip_grid) :: ip_equid_cylind_grid
     real :: hi
     real :: rlat1
     real :: rlon1
     real :: rlat2
     real :: rlon2
     real :: dlat
     real :: dlon
   contains
     procedure :: init_grib1
     procedure :: init_grib2
     procedure :: gdswzd => gdswzd_equid_cylind
  end type ip_equid_cylind_grid
 
  REAL :: dlat
  REAL :: dlon
  REAL :: rerth
 
contains
 
  subroutine init_grib1(self, g1_desc)
    class(ip_equid_cylind_grid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    integer :: iscan
 
    associate(kgds => g1_desc%gds)
      self%IM=kgds(2)
      self%JM=kgds(3)
      self%RLAT1=kgds(4)*1.e-3
      self%RLON1=kgds(5)*1.e-3
      self%RLAT2=kgds(7)*1.e-3
      self%RLON2=kgds(8)*1.e-3
      iscan=mod(kgds(11)/128,2)
      self%HI=(-1.)**iscan
      self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
      self%DLAT=(self%RLAT2-self%RLAT1)/(self%JM-1)
 
      ! defaults
      self%iwrap = 0
      self%jwrap1 = 0
      self%jwrap2 = 0
      self%nscan = mod(kgds(11) / 32, 2)
      self%nscan_field_pos = self%nscan
      self%kscan = 0
 
      self%iwrap = nint(360/abs(self%dlon))
 
      if(self%im < self%iwrap) self%iwrap=0
      self%jwrap1 = 0
      self%jwrap2 = 0
      if(self%iwrap > 0 .and. mod(self%iwrap,2) == 0) then
         if(abs(self%rlat1) > 90-0.25*self%dlat) then
            self%jwrap1 = 2
         elseif(abs(self%rlat1) > 90-0.75*self%dlat) then
            self%jwrap1 = 1
         endif
         if(abs(self%rlat2) > 90-0.25*self%dlat) then
            self%jwrap2 = 2 * self%jm
         elseif(abs(self%rlat2) > 90-0.75*self%dlat) then
            self%jwrap2 = 2 * self%jm+1
         endif
      endif
 
      self%rerth = 6.3712e6
      self%eccen_squared = 0.0
    end associate
 
  end subroutine init_grib1
  
  subroutine init_grib2(self, g2_desc)
    class(ip_equid_cylind_grid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    integer :: iscale, iscan
 
    associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
      self%IM=igdtmpl(8)
      self%JM=igdtmpl(9)
      iscale=igdtmpl(10)*igdtmpl(11)
      IF(iscale==0) iscale=10**6
      self%RLAT1=float(igdtmpl(12))/float(iscale)
      self%RLON1=float(igdtmpl(13))/float(iscale)
      self%RLAT2=float(igdtmpl(15))/float(iscale)
      self%RLON2=float(igdtmpl(16))/float(iscale)
      iscan=mod(igdtmpl(19)/128,2)
      self%HI=(-1.)**iscan
      self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
      self%DLAT=(self%RLAT2-self%RLAT1)/(self%JM-1)
 
      self%nscan = mod(igdtmpl(19)/32,2)
      self%nscan_field_pos = self%nscan
      self%kscan = 0
      self%iwrap = nint(360/abs(self%DLON))
 
      if(self%im.lt.self%iwrap) self%iwrap=0
      self%jwrap1=0
      self%jwrap2=0
 
      if(self%im < self%iwrap) self%iwrap=0
      self%jwrap1 = 0
      self%jwrap2 = 0
      if(self%iwrap > 0 .and. mod(self%iwrap,2) == 0) then
         if(abs(self%rlat1) > 90-0.25*self%dlat) then
            self%jwrap1 = 2
         elseif(abs(self%rlat1) > 90-0.75*self%dlat) then
            self%jwrap1 = 1
         endif
         if(abs(self%rlat2) > 90-0.25*self%dlat) then
            self%jwrap2 = 2 * self%jm
         elseif(abs(self%rlat2) > 90-0.75*self%dlat) then
            self%jwrap2 = 2 * self%jm+1
         endif
      endif
 
      call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
 
    end associate
  end subroutine init_grib2
 
  SUBROUTINE gdswzd_equid_cylind(self,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET,  &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    IMPLICIT NONE
    !
    class(ip_equid_cylind_grid), 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, N
    !
    LOGICAL                            :: LROT, LMAP, LAREA
    !
    REAL                               :: HI, RLAT1, RLON1, RLAT2, RLON2
    REAL                               :: XMAX, XMIN, YMAX, YMIN
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    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
 
    im=self%im
    jm=self%jm
 
    rlat1=self%rlat1
    rlon1=self%rlon1
    rlat2=self%rlat2
    rlon2=self%rlon2
 
    hi=self%hi
 
    rerth = self%rerth
    dlat = self%dlat
    dlon = self%dlon
 
    xmin=0
    xmax=im+1
    IF(im.EQ.nint(360/abs(dlon))) 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
       !$OMP PARALLEL DO PRIVATE(N) REDUCTION(+:NRET) SCHEDULE(STATIC)
       DO n=1,npts
          IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
               ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
             rlon(n)=mod(rlon1+dlon*(xpts(n)-1)+3600,360.)
             rlat(n)=min(max(rlat1+dlat*(ypts(n)-1),-90.),90.)
             nret=nret+1
             IF(lrot)  CALL equid_cylind_vect_rot(crot(n),srot(n))
             IF(lmap)  CALL equid_cylind_map_jacob(xlon(n),xlat(n),ylon(n),ylat(n))
             IF(larea) CALL equid_cylind_grid_area(rlat(n),area(n))
          ELSE
             rlon(n)=fill
             rlat(n)=fill
          ENDIF
       ENDDO
       !$OMP END PARALLEL DO
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
    ELSEIF(iopt.EQ.-1) THEN
       !$OMP PARALLEL DO PRIVATE(N) REDUCTION(+:NRET) SCHEDULE(STATIC)
       DO n=1,npts
          IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90) THEN
             xpts(n)=1+hi*mod(hi*(rlon(n)-rlon1)+3600,360.)/dlon
             ypts(n)=1+(rlat(n)-rlat1)/dlat
             IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND.  &
                  ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
                nret=nret+1
                IF(lrot)  CALL equid_cylind_vect_rot(crot(n),srot(n))
                IF(lmap)  CALL equid_cylind_map_jacob(xlon(n),xlat(n),ylon(n),ylat(n))
                IF(larea) CALL equid_cylind_grid_area(rlat(n),area(n))
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ELSE
             xpts(n)=fill
             ypts(n)=fill
          ENDIF
       ENDDO
       !$OMP END PARALLEL DO
    ENDIF
  END SUBROUTINE gdswzd_equid_cylind
 
  SUBROUTINE equid_cylind_vect_rot(CROT,SROT)
    IMPLICIT NONE
 
    REAL,                INTENT(  OUT) :: CROT, SROT
 
    crot=1.0
    srot=0.0
 
  END SUBROUTINE equid_cylind_vect_rot
 
  SUBROUTINE equid_cylind_map_jacob(XLON,XLAT,YLON,YLAT)
    REAL,                INTENT(  OUT) :: XLON,XLAT,YLON,YLAT
 
    xlon=1.0/dlon
    xlat=0.
    ylon=0.
    ylat=1.0/dlat
 
  END SUBROUTINE equid_cylind_map_jacob
 
  SUBROUTINE equid_cylind_grid_area(RLAT,AREA)
    IMPLICIT NONE
 
    REAL,                INTENT(IN   ) :: RLAT
    REAL,                INTENT(  OUT) :: AREA
 
    REAL,                PARAMETER     :: PI=3.14159265358979
    REAL,                PARAMETER     :: DPR=180./pi
 
    REAL                               :: DSLAT, RLATU, RLATD
 
    rlatu=min(max(rlat+dlat/2,-90.),90.)
    rlatd=min(max(rlat-dlat/2,-90.),90.)
    dslat=sin(rlatu/dpr)-sin(rlatd/dpr)
    area=rerth**2*abs(dslat*dlon)/dpr
 
  END SUBROUTINE equid_cylind_grid_area
 
end module ip_equid_cylind_grid_mod