Skip to content

Commit

Permalink
Fix Lambert Conformal Conic mapping in the write component on the sou…
Browse files Browse the repository at this point in the history 
…thern hemisphere (#497)

Subroutine `lambert` in the write component has been fixed to do the mapping on the southern hemisphere correctly.
  • Loading branch information
@DusanJovic-NOAA
DusanJovic-NOAA authored Mar 18, 2022
1 parent 8c582f7 commit 1cebcf1
Showing 1 changed file with 39 additions and 29 deletions.
68 changes: 39 additions & 29 deletions io/module_wrt_grid_comp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3540,45 +3540,55 @@ end subroutine rtll
subroutine lambert(stlat1,stlat2,c_lat,c_lon,glon,glat,x,y,inv)

!-------------------------------------------------------------------------------
real(ESMF_KIND_R8), intent(in) :: stlat1,stlat2,c_lat,c_lon
real(ESMF_KIND_R8), intent(inout) :: glon, glat
real(ESMF_KIND_R8), intent(in) :: stlat1,stlat2,c_lat,c_lon
real(ESMF_KIND_R8), intent(inout) :: glon, glat
real(ESMF_KIND_R8), intent(inout) :: x, y
integer, intent(in) :: inv
integer, intent(in) :: inv
!-------------------------------------------------------------------------------
! real(ESMF_KIND_R8), parameter :: pi=3.14159265358979323846
real(ESMF_KIND_R8), parameter :: dtor=pi/180.0
real(ESMF_KIND_R8), parameter :: rtod=180.0/pi
! real(ESMF_KIND_R8), parameter :: pi = 3.14159265358979323846
real(ESMF_KIND_R8), parameter :: dtor = pi/180.0
real(ESMF_KIND_R8), parameter :: rtod = 180.0/pi
real(ESMF_KIND_R8), parameter :: a = 6371200.0
!-------------------------------------------------------------------------------
! inv == 1 (glon,glat) ---> (x,y) lat/lon to grid
! inv == -1 (x,y) ---> (glon,glat) grid to lat/lon
! inv == 1 (glon,glat) ---> (x,y)
! inv == -1 (x,y) ---> (glon,glat)

real(ESMF_KIND_R8) :: en,f,rho,rho0, dlon, theta
real(ESMF_KIND_R8) :: xp, yp, en, de, rho, rho0, rho2, dlon, theta, dr2
real(ESMF_KIND_R8) :: h = 1.0

IF (stlat1 == stlat2) THEN
en=sin(stlat1*dtor)
ELSE
en=log(cos(stlat1*dtor)/cos(stlat2*dtor))/ &
log(tan((45+0.5*stlat2)*dtor)/tan((45+0.5*stlat1)*dtor))
ENDIF
! For reference see:
! John P. Snyder (1987), Map projections: A working manual (pp. 104-110)
! https://doi.org/10.3133/pp1395

f=(cos(stlat1*dtor)*tan((45+0.5*stlat1)*dtor)**en)/en
rho0=a*f/(tan((45+0.5*c_lat)*dtor)**en)
if (stlat1 == stlat2) then
en = sin(stlat1*dtor)
else
en = log(cos(stlat1*dtor)/cos(stlat2*dtor)) / &
log(tan((45+0.5*stlat2)*dtor)/tan((45+0.5*stlat1)*dtor)) ! (15-3)
endif
h = sign(1.0_ESMF_KIND_R8,en)

de = a*(cos(stlat1*dtor)*tan((45+0.5*stlat1)*dtor)**en)/en ! (15-2)
rho0 = de/(tan((45+0.5*c_lat)*dtor)**en) ! (15-1a)

if (inv == 1) then ! FORWARD TRANSFORMATION
rho=a*f/(tan((45+0.5*glat)*dtor)**en)
dlon=modulo(glon-c_lon+180+3600,360.)-180.D0
theta=en*dlon*dtor
x=rho*sin(theta)
y=rho0-rho*cos(theta)
rho = de/(tan((45+0.5*glat)*dtor)**en) ! (15-1)
dlon = modulo(glon-c_lon+180.0+3600.0,360.0)-180.0
theta = en*dlon*dtor ! (14-4)
x = rho*sin(theta) ! (14-1)
y = rho0-rho*cos(theta) ! (14-2)
else if (inv == -1) then ! INVERSE TRANSFORMATION
y=rho0-y
rho = sqrt(x*x+y*y)
theta=atan2(x,y)
glon=c_lon+(theta/en)*rtod
glon=modulo(glon+180+3600,360.)-180.D0
! glat=(2.0*atan((a*f/rho)**(1.0/en))-0.5*pi)*rtod
glat=(0.5*pi-2.0*atan((rho/(a*f))**(1.0/en)))*rtod
xp = h*x;
yp = h*(rho0-y)
theta = atan2(xp,yp) ! (14-11)
glon = c_lon+(theta/en)*rtod ! (14-9)
glon = modulo(glon+180.0+3600.0,360.0)-180.0
rho2 = xp*xp+yp*yp ! (14-10)
if (rho2 == 0.0) then
glat = h*90.0
else
glat = 2.0*atan((de*de/rho2)**(1.0/(2.0*en)))*rtod-90.0 ! (15-5)
endif
else
write (unit=*,fmt=*) " lambert: unknown inv argument"
return
Expand Down

0 comments on commit 1cebcf1

Please sign in to comment.