ip_lambert_conf_grid_mod.F90

Go to the documentation of this file.

 
 
module ip_lambert_conf_grid_mod
  use ip_grid_descriptor_mod
  use ip_grid_mod
  use earth_radius_mod
  use ip_constants_mod
  implicit none
 
  private
  public :: ip_lambert_conf_grid
 
  type, extends(ip_grid) :: ip_lambert_conf_grid
     real :: rlat1
     real :: rlon1
     real :: rlati1
     real :: rlati2
     real :: orient
     real :: dxs
     real :: dys
     real :: h
     integer :: irot
   contains
     procedure :: init_grib1
     procedure :: init_grib2
     procedure :: gdswzd => gdswzd_lambert_conf
  end type ip_lambert_conf_grid
 
  INTEGER :: irot
  REAL :: an
  REAL :: dxs
  REAL :: dys
  REAL :: h
  REAL :: rerth
  REAL :: tinyreal=tiny(1.0) 
 
contains
 
  subroutine init_grib1(self, g1_desc)
    class(ip_lambert_conf_grid), intent(inout) :: self
    type(grib1_descriptor), intent(in) :: g1_desc
 
    real :: dx, dy, hi, hj
    integer :: iproj, iscan, jscan
 
    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%IROT=mod(kgds(6)/8,2)
      self%ORIENT=kgds(7)*1.e-3
 
      dx=kgds(8)
      dy=kgds(9)
 
      iproj=mod(kgds(10)/128,2)
      iscan=mod(kgds(11)/128,2)
      jscan=mod(kgds(11)/64,2)
 
      self%RLATI1=kgds(12)*1.e-3
      self%RLATI2=kgds(13)*1.e-3
      self%H=(-1.)**iproj
 
      hi=(-1.)**iscan
      hj=(-1.)**(1-jscan)
      self%DXS=dx*hi
      self%DYS=dy*hj
 
      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
    end associate
 
  end subroutine init_grib1
 
  subroutine init_grib2(self, g2_desc)
    class(ip_lambert_conf_grid), intent(inout) :: self
    type(grib2_descriptor), intent(in) :: g2_desc
 
    real :: dx, dy, hi, hj
    integer :: iproj, iscan, jscan
 
 
    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)
 
      self%RLAT1=float(igdtmpl(10))*1.0e-6
      self%RLON1=float(igdtmpl(11))*1.0e-6
 
      self%IROT=mod(igdtmpl(12)/8,2)
      self%ORIENT=float(igdtmpl(14))*1.0e-6
 
      dx=float(igdtmpl(15))*1.0e-3
      dy=float(igdtmpl(16))*1.0e-3
 
      iproj=mod(igdtmpl(17)/128,2)
      iscan=mod(igdtmpl(18)/128,2)
      jscan=mod(igdtmpl(18)/64,2)
 
      self%RLATI1=float(igdtmpl(19))*1.0e-6
      self%RLATI2=float(igdtmpl(20))*1.0e-6
 
      self%H=(-1.)**iproj
      hi=(-1.)**iscan
      hj=(-1.)**(1-jscan)
      self%DXS=dx*hi
      self%DYS=dy*hj
 
      self%nscan = mod(igdtmpl(18) / 32, 2)
      self%nscan_field_pos = self%nscan
      self%iwrap = 0
      self%jwrap1 = 0
      self%jwrap2 = 0
      self%kscan = 0
    end associate
  end subroutine init_grib2
 
  SUBROUTINE gdswzd_lambert_conf(self,IOPT,NPTS,FILL, &
       XPTS,YPTS,RLON,RLAT,NRET, &
       CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
    IMPLICIT NONE
    !
    class(ip_lambert_conf_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                          :: ANTR, DI, DJ
    REAL                          :: DLON1
    REAL                          :: DE, DE2, DR2
    REAL                          :: ORIENT, RLAT1, RLON1
    REAL                          :: RLATI1, RLATI2
    REAL                          :: XMAX, XMIN, YMAX, YMIN, XP, YP
    REAL                          :: DLON, DR
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    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
 
    irot=self%irot
    orient=self%orient
 
    rlati1=self%rlati1
    rlati2=self%rlati2
 
    h=self%h
    dxs=self%dxs
    dys=self%dys
 
    rerth = self%rerth
 
    IF(abs(rlati1-rlati2).LT.tinyreal) THEN
       an=sin(rlati1/dpr)
    ELSE
       an=log(cos(rlati1/dpr)/cos(rlati2/dpr))/ &
            log(tan((90-rlati1)/2/dpr)/tan((90-rlati2)/2/dpr))
    ENDIF
    de=rerth*cos(rlati1/dpr)*tan((rlati1+90)/2/dpr)**an/an
    IF(abs(h*rlat1-90).LT.tinyreal) THEN
       xp=1
       yp=1
    ELSE
       dr=de/tan((rlat1+90)/2/dpr)**an
       dlon1=mod(rlon1-orient+180+3600,360.)-180
       xp=1-sin(an*dlon1/dpr)*dr/dxs
       yp=1+cos(an*dlon1/dpr)*dr/dys
    ENDIF
    antr=1/(2*an)
    de2=de**2
    xmin=0
    xmax=im+1
    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,DI,DJ,DR2,DR,DLON) 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
             di=h*(xpts(n)-xp)*dxs
             dj=h*(ypts(n)-yp)*dys
             dr2=di**2+dj**2
             dr=sqrt(dr2)
             IF(dr2.LT.de2*1.e-6) THEN
                rlon(n)=0.
                rlat(n)=h*90.
             ELSE
                rlon(n)=mod(orient+1./an*dpr*atan2(di,-dj)+3600,360.)
                rlat(n)=(2*dpr*atan((de2/dr2)**antr)-90)
             ENDIF
             nret=nret+1
             dlon=mod(rlon(n)-orient+180+3600,360.)-180
             IF(lrot)  CALL lambert_conf_vect_rot(dlon,crot(n),srot(n))
             IF(lmap)  CALL lambert_conf_map_jacob(rlat(n),fill, dlon, dr, &
                  xlon(n),xlat(n),ylon(n),ylat(n))
             IF(larea) CALL lambert_conf_grid_area(rlat(n),fill,dr,area(n))
          ELSE
             rlon(n)=fill
             rlat(n)=fill
          ENDIF
       ENDDO
       ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       !  TRANSLATE EARTH COORDINATES TO GRID COORDINATES
    ELSEIF(iopt.EQ.-1) THEN
       !$OMP PARALLEL DO PRIVATE(N,DR,DLON) REDUCTION(+:NRET) SCHEDULE(STATIC)
       DO n=1,npts
          IF(abs(rlon(n)).LT.(360.+tinyreal).AND.abs(rlat(n)).LT.(90.+tinyreal).AND. &
               abs(h*rlat(n)+90).GT.tinyreal) THEN
             dr=h*de*tan((90-rlat(n))/2/dpr)**an
             dlon=mod(rlon(n)-orient+180+3600,360.)-180
             xpts(n)=xp+h*sin(an*dlon/dpr)*dr/dxs
             ypts(n)=yp-h*cos(an*dlon/dpr)*dr/dys
             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 lambert_conf_vect_rot(dlon,crot(n),srot(n))
                IF(lmap)  CALL lambert_conf_map_jacob(rlat(n),fill,dlon,dr, &
                     xlon(n),xlat(n),ylon(n),ylat(n))
                IF(larea) CALL lambert_conf_grid_area(rlat(n),fill,dr,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_lambert_conf
 
  SUBROUTINE lambert_conf_vect_rot(DLON,CROT,SROT)
    IMPLICIT NONE
    REAL,           INTENT(   IN) :: DLON
    REAL,           INTENT(  OUT) :: CROT, SROT
 
    IF(irot.EQ.1) THEN
       crot=cos(an*dlon/dpr)
       srot=sin(an*dlon/dpr)
    ELSE
       crot=1.
       srot=0.
    ENDIF
 
  END SUBROUTINE lambert_conf_vect_rot
 
  SUBROUTINE lambert_conf_map_jacob(RLAT,FILL,DLON,DR,XLON,XLAT,YLON,YLAT)
    IMPLICIT NONE
 
    REAL,           INTENT(IN   ) :: RLAT, FILL, DLON, DR
    REAL,           INTENT(  OUT) :: XLON, XLAT, YLON, YLAT
 
    REAL                          :: CLAT
 
    clat=cos(rlat/dpr)
    IF(clat.LE.0.OR.dr.LE.0) THEN
       xlon=fill
       xlat=fill
       ylon=fill
       ylat=fill
    ELSE
       xlon=h*cos(an*dlon/dpr)*an/dpr*dr/dxs
       xlat=-h*sin(an*dlon/dpr)*an/dpr*dr/dxs/clat
       ylon=h*sin(an*dlon/dpr)*an/dpr*dr/dys
       ylat=h*cos(an*dlon/dpr)*an/dpr*dr/dys/clat
    ENDIF
 
  END SUBROUTINE lambert_conf_map_jacob
 
  SUBROUTINE lambert_conf_grid_area(RLAT,FILL,DR,AREA)
    IMPLICIT NONE
 
    REAL,           INTENT(IN   ) :: RLAT
    REAL,           INTENT(IN   ) :: FILL
    REAL,           INTENT(IN   ) :: DR
    REAL,           INTENT(  OUT) :: AREA
 
    REAL                          :: CLAT
 
    clat=cos(rlat/dpr)
    IF(clat.LE.0.OR.dr.LE.0) THEN
       area=fill
    ELSE
       area=rerth**2*clat**2*abs(dxs)*abs(dys)/(an*dr)**2
    ENDIF
 
  END SUBROUTINE lambert_conf_grid_area
 
end module ip_lambert_conf_grid_mod