ip_gaussian_grid_mod.F90

Go to the documentation of this file.

 
 
module ip_gaussian_grid_mod
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use earth_radius_mod
  use constants_mod
  implicit none
 
  private
  public :: ip_gaussian_grid
 
  type, extends(ip_grid) :: ip_gaussian_grid
     integer :: jh
     real :: dlon
     real :: rlat1
     real :: rlon1
     real :: rlon2
     real :: hi
     integer :: jg
     integer :: jscan
   contains
     procedure :: init_grib1
     procedure :: init_grib2
     procedure :: gdswzd => gdswzd_gaussian
  end type ip_gaussian_grid
 
  INTEGER :: j1
  INTEGER :: jh
  REAL, ALLOCATABLE :: blat(:)
  REAL :: dlon
  REAL :: rerth
  REAL, ALLOCATABLE :: ylat_row(:)
 
contains
 
  subroutine init_grib1(self, g1_desc)
    class(ip_gaussian_grid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    integer :: iscan, jg
 
    associate(kgds => g1_desc%gds)
      self%rerth = 6.3712e6
      self%eccen_squared = 0.0
 
      self%IM=kgds(2)
      self%JM=kgds(3)
      self%RLAT1=kgds(4)*1.e-3
      self%RLON1=kgds(5)*1.e-3
      self%RLON2=kgds(8)*1.e-3
      self%JG=kgds(10)*2
      iscan=mod(kgds(11)/128,2)
      self%JSCAN=mod(kgds(11)/64,2)
      self%HI=(-1.)**iscan
      self%JH=(-1)**self%JSCAN
      self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
 
      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
 
      if(self%iwrap > 0 .and. mod(self%iwrap, 2) == 0) then
         jg=kgds(10)*2
         if(self%jm == self%jg) then
            self%jwrap1 = 1
            self%jwrap2 = 2 * self%jm + 1
         endif
      endif
 
    end associate
  end subroutine init_grib1
 
  subroutine init_grib2(self, g2_desc)
    class(ip_gaussian_grid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    integer :: iscale, iscan, jg
 
    associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
      call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
 
      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%RLON2=float(igdtmpl(16))/float(iscale)
      self%JG=igdtmpl(18)*2
      iscan=mod(igdtmpl(19)/128,2)
      self%JSCAN=mod(igdtmpl(19)/64,2)
      self%HI=(-1.)**iscan
      self%JH=(-1)**self%JSCAN
      self%DLON=self%HI*(mod(self%HI*(self%RLON2-self%RLON1)-1+3600,360.)+1)/(self%IM-1)
 
 
      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
         jg = igdtmpl(18) * 2
         if(self%jm == jg) then
            self%jwrap1=1
            self%jwrap2 = 2 * self%jm + 1
         endif
      endif
      self%nscan = mod(igdtmpl(19) / 32, 2)
      self%nscan_field_pos = self%nscan
      self%kscan = 0
    end associate
 
  end subroutine init_grib2
 
  SUBROUTINE gdswzd_gaussian(self,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    IMPLICIT NONE
    !
    class(ip_gaussian_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                        :: JSCAN, IM, JM
    INTEGER                        :: J, JA, JG
    INTEGER                        :: N
    !
    LOGICAL                        :: LROT, LMAP, LAREA
    !
    REAL,            ALLOCATABLE   :: ALAT(:), ALAT_JSCAN(:)
    REAL,            ALLOCATABLE   :: ALAT_TEMP(:),BLAT_TEMP(:)
    REAL                           :: HI, RLATA, RLATB, RLAT1, RLON1, RLON2
    REAL                           :: XMAX, XMIN, YMAX, YMIN, YPTSA, YPTSB
    REAL                           :: WB
 
    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
 
    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
 
    im=self%im
    jm=self%jm
 
    rlat1=self%rlat1
    rlon1=self%rlon1
    rlon2=self%rlon2
 
    jg=self%jg
    jscan=self%jscan
    hi=self%hi
 
    jh=self%jh
    dlon=self%dlon
    rerth = self%rerth
 
    ALLOCATE(alat_temp(jg))
    ALLOCATE(blat_temp(jg))
    CALL splat(4,jg,alat_temp,blat_temp)
    ALLOCATE(alat(0:jg+1))
    ALLOCATE(blat(0:jg+1))
    !$OMP PARALLEL DO PRIVATE(JA) SCHEDULE(STATIC)
    DO ja=1,jg
       alat(ja)=real(dpr*asin(alat_temp(ja)))
       blat(ja)=blat_temp(ja)
    ENDDO
    !$OMP END PARALLEL DO
    DEALLOCATE(alat_temp,blat_temp)
    alat(0)=180.-alat(1)
    alat(jg+1)=-alat(0)
    blat(0)=-blat(1)
    blat(jg+1)=blat(0)
    j1=1
    DO WHILE(j1.LT.jg.AND.rlat1.LT.(alat(j1)+alat(j1+1))/2)
       j1=j1+1
    ENDDO
    IF(lmap)THEN
       ALLOCATE(alat_jscan(jg))
       DO ja=1,jg
          alat_jscan(j1+jh*(ja-1))=alat(ja)
       ENDDO
       ALLOCATE(ylat_row(0:jg+1))
       DO ja=2,(jg-1)
          ylat_row(ja)=2.0/(alat_jscan(ja+1)-alat_jscan(ja-1))
       ENDDO
       ylat_row(1)=1.0/(alat_jscan(2)-alat_jscan(1))
       ylat_row(0)=ylat_row(1)
       ylat_row(jg)=1.0/(alat_jscan(jg)-alat_jscan(jg-1))
       ylat_row(jg+1)=ylat_row(jg)
       DEALLOCATE(alat_jscan)
    ENDIF
    xmin=0
    xmax=im+1
    IF(im.EQ.nint(360/abs(dlon))) xmax=im+2
    ymin=0.5
    ymax=jm+0.5
    nret=0
 
    !  TRANSLATE GRID COORDINATES TO EARTH COORDINATES
    IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
       !$OMP PARALLEL DO PRIVATE(N,J,WB,RLATA,RLATB) 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.)
             j=int(ypts(n))
             wb=ypts(n)-j
             rlata=alat(j1+jh*(j-1))
             rlatb=alat(j1+jh*j)
             rlat(n)=rlata+wb*(rlatb-rlata)
             nret=nret+1
             IF(lrot) CALL gaussian_vect_rot(crot(n),srot(n))
             IF(lmap) CALL gaussian_map_jacob(ypts(n),&
                  xlon(n),xlat(n),ylon(n),ylat(n))
             IF(larea) CALL gaussian_grid_area(ypts(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,JA,YPTSA, YPTSB, WB) REDUCTION(+:NRET) SCHEDULE(STATIC)
       DO n=1,npts
          xpts(n)=fill
          ypts(n)=fill
          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
             ja=min(int((jg+1)/180.*(90-rlat(n))),jg)
             IF(rlat(n).GT.alat(ja)) ja=max(ja-2,0)
             IF(rlat(n).LT.alat(ja+1)) ja=min(ja+2,jg)
             IF(rlat(n).GT.alat(ja)) ja=ja-1
             IF(rlat(n).LT.alat(ja+1)) ja=ja+1
             yptsa=1+jh*(ja-j1)
             yptsb=1+jh*(ja+1-j1)
             wb=(alat(ja)-rlat(n))/(alat(ja)-alat(ja+1))
             ypts(n)=yptsa+wb*(yptsb-yptsa)
             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 gaussian_vect_rot(crot(n),srot(n))
                IF(lmap) CALL gaussian_map_jacob(ypts(n), &
                     xlon(n),xlat(n),ylon(n),ylat(n))
                IF(larea) CALL gaussian_grid_area(ypts(n),area(n))
             ELSE
                xpts(n)=fill
                ypts(n)=fill
             ENDIF
          ENDIF
       ENDDO
       !$OMP END PARALLEL DO
    ENDIF
    DEALLOCATE(alat, blat)
    IF (ALLOCATED(ylat_row)) DEALLOCATE(ylat_row)
 
  END SUBROUTINE gdswzd_gaussian
 
  SUBROUTINE gaussian_vect_rot(CROT,SROT)
    IMPLICIT NONE
 
    REAL,                INTENT(  OUT) :: CROT, SROT
 
    crot=1.0
    srot=0.0
 
  END SUBROUTINE gaussian_vect_rot
 
  SUBROUTINE gaussian_map_jacob(YPTS, XLON, XLAT, YLON, YLAT)
    IMPLICIT NONE
 
    REAL,                INTENT(IN   ) :: YPTS
    REAL,                INTENT(  OUT) :: XLON, XLAT, YLON, YLAT
 
    xlon=1/dlon
    xlat=0.
    ylon=0.
    ylat=ylat_row(nint(ypts))
 
  END SUBROUTINE gaussian_map_jacob
 
  SUBROUTINE gaussian_grid_area(YPTS,AREA)
    IMPLICIT NONE
 
    REAL,            INTENT(IN   ) :: YPTS
    REAL,            INTENT(  OUT) :: AREA
 
    INTEGER                        :: J
 
    REAL                           :: WB, WLAT, WLATA, WLATB
 
    j = int(ypts)
    wb=ypts-j
    wlata=blat(j1+jh*(j-1))
    wlatb=blat(j1+jh*j)
    wlat=wlata+wb*(wlatb-wlata)
    area=real(rerth**2*wlat*dlon/dpr)
 
  END SUBROUTINE gaussian_grid_area
end module ip_gaussian_grid_mod