! Program Name: ! Author(s)/Contact(s): ! Abstract: ! History Log: ! ! Usage: ! Parameters: ! Input Files: ! ! Output Files: ! ! ! Condition codes: ! ! If appropriate, descriptive troubleshooting instructions or ! likely causes for failures could be mentioned here with the ! appropriate error code ! ! User controllable options: module module_llxy_generic real, parameter, private :: pi = 3.141592653589793 ! you know! pi = 180 degrees real, parameter, private :: piovr2 = pi/2. ! real, parameter, private :: re = 6370.949 ! the radius of the earth in km real, parameter, private :: re = 6371.2 ! the radius of the earth in km real, parameter, private :: twopi = pi*2. real, parameter, private :: degrad = pi/180. real, parameter, private :: ce = 2.*pi*re ! the circumference of the earth in km contains !======================================================================================= !======================================================================================= subroutine lltoxy_generic (xlat,xlon,x,y, & project,dskm,reflat,reflon,refx,refy,& cenlon,truelat1,truelat2) !*****************************************************************************! ! Notes - Modeled after XYTOLL in the plots.o library ! !*****************************************************************************! implicit none ! Input: ---------------------------------------------------------------------! real :: xlat ! latitude of the point of interest. real :: xlon ! longitude of the point of interest. character(LEN=2) :: project ! projection indicator ("ME", "CE", "LC", "ST") real :: dskm ! grid distance in km real :: reflat real :: reflon real :: refx real :: refy real :: cenlon ! Grid ratio with respect to MOAD real :: truelat1 ! True latitude 1, closest to equator real :: truelat2 ! True latitude 2, closest to pole ! Output: --------------------------------------------------------------------! real :: x ! x location of the given (lat,lon) point real :: y ! y location of the given (lat,lon) point ! Real variables: -----------------------------------------------------------! real :: flat1 real :: confac ! cone factor real :: rcln ! center longitude in radians (local) real :: dj ! distance from pole to point (local) real :: di ! distance from the central ! meridian to the point (local) real :: bm ! calculation variable (local) real :: rflt, rfln, rlat, rlon, diovrdj, disdjs, ri, rj real :: ct1, st1, tt1, drp real :: isn, pi2 real :: londiff !**************************** subroutine begin *****************************! rlat = xlat * degrad rlon = xlon * degrad rcln = cenlon * degrad flat1 = truelat1 * degrad rflt = reflat * degrad rfln = reflon * degrad if ((xlon - cenlon) < -180) then londiff = ((xlon-cenlon)+360.)*degrad else if ((xlon - cenlon) > 180) then londiff = ((xlon-cenlon)-360.)*degrad else londiff = (xlon-cenlon)*degrad endif if (project(1:2) .eq. 'ME') then ct1 = re*cos(flat1) dj = ct1 * log(tan (0.5*(rlat + piovr2))) di = ct1 * (rlon - rfln) y = refy +(dj + ct1 * log(cos(rflt)/(1 + sin(rflt))))/dskm x = refx + di/dskm else if (project(1:2) .eq. 'CE') then !KWM dj = re*(rlat-rflt) !KWM di = re*(rlon-rfln) !KWM y = refy + dj/dskm !KWM x = refx + di/dskm ! NOTE: Cylindrical Equidistant grid increment expressed in terms ! of thousandths of degrees! ! Calculate the distance from the horizontal axis to (J,I) dj = xlat-reflat di = xlon-reflon y = refy + (dj/dskm) x = refx + (di/dskm) else if (project(1:2) .eq. 'LC') then isn = sign(1.0, truelat2) pi2 = piovr2*isn call lccone(truelat1,truelat2,int(sign(1.0, truelat2)),confac) tt1 = tan((pi2 - flat1)*0.5) ! Tangent Term 1. st1 = sin (pi2 - flat1) * re/(confac*dskm) ! Sine Term 1. bm = tan((pi2 - rlat)*0.5) if ((rlon-rcln) > pi) then diovrdj = -tan(((rlon-rcln)-2.*pi)*confac) elseif ((rlon-rcln) < -pi) then diovrdj = -tan(((rlon-rcln)+2.*pi)*confac) else diovrdj = -tan((rlon-rcln)*confac) endif disdjs = ( (bm/tt1)**confac * st1)**2 ! Dj: y distance (km) from pole to given x/y point. dj = -sqrt(disdjs/(1+diovrdj*diovrdj)) ! Di: x distance (km) from central longitude to given x/y point. di = dj * diovrdj bm = tan((pi2-rflt)*0.5) if ((rfln-rcln) > pi) then diovrdj = -tan(((rfln-rcln)-2.*pi)*confac) else if ((rfln-rcln) < -pi) then diovrdj = -tan(((rfln-rcln)+2.*pi)*confac) else diovrdj = -tan((rfln - rcln)*confac) endif disdjs = ( (bm/tt1)**confac * st1)**2 ! Rj: y distance (km) from pole to reference point. rj = -sqrt(disdjs/(1+diovrdj*diovrdj)) ! Ri: x distance (km) from central longitude to reference x/y point. ri = rj * diovrdj y = refy + isn*(dj - rj) x = refx + (di - ri) else if (project(1:2) .eq. 'ST') then isn = sign(1.0, truelat1) pi2 = piovr2*isn diovrdj = -tan(rfln-rcln) bm = (re/dskm)*tan(0.5*(pi2-rflt)) disdjs = (bm*(1.0 + cos(pi2-flat1)))**2 ! RJ: Distance from pole to reference latitude along the center lon rj = -sign(sqrt(disdjs/(1+diovrdj**2)),cos(rfln-rcln)) ! RI: Distance from center reference to requested point rlat, rlon. ri = rj * diovrdj diovrdj = -tan(rlon-rcln) bm = (re/dskm)*tan(0.5*(pi2-rlat)) disdjs = (bm*(1.0 + cos(pi2-flat1)))**2 dj = -sign(sqrt(disdjs/(1+diovrdj**2)),cos(rlon-rcln)) di = dj * diovrdj y = refy + isn*(dj-rj) x = refx + (di-ri) endif !***************************** Subroutine End ******************************! end subroutine lltoxy_generic !======================================================================================= !======================================================================================= subroutine xytoll_generic (x,y,xlat,xlon,& project,dskm,reflat,reflon,refx,refy,& cenlon, truelat1,truelat2) !*****************************************************************************! ! xytoll ! ! Purpose - To transform mesoscale grid point coordinates (x,y) into ! ! latitude and longitude coordinates. ! ! ! ! On entry - X and Y are an ordered pair representing a grid point in the ! ! mesoscale grid. ! ! ! ! On exit - XLAT, XLON contain the latitude and longitude respectively ! ! that resulted from the transformation. ! ! ! !*****************************************************************************! implicit none ! Input real :: x ! x (i) coordinate of point of interest real :: y ! y (j) coordinate of point of interest character(LEN=2) :: project ! projection indicator ("ME", "CE", "LC", "ST") real :: dskm ! grid distance in km real :: reflat ! latitude of the reference point real :: reflon ! longitude of the reference point real :: refx real :: refy real :: cenlon ! longitude of the center of the projection real :: truelat1 ! True latitude 1. real :: truelat2 ! True latitude 2. ! Output real :: xlat ! latitude of point (x,y) real :: xlon ! longitude of point (x,y) ! Real variables real :: confac ! cone factor real :: rfln ! reference longitude in radians (local) real :: rflt ! reference latitude in radians (local) real :: rcln ! center longitude in radians (local) real :: dj, drp, djj ! distance from the central ! meridian to the point (local) real :: di,dii ! distance from pole to point (local) real :: bm ! calculation variable (local) real :: flat1 real :: ct1, tt1, tt2, pi2 integer :: isn !**************************** subroutine begin *****************************! rflt = reflat * degrad rfln = reflon * degrad rcln = cenlon * degrad flat1 = truelat1 * degrad ! If the projection is mercator ('ME') then ... if (project(1:2) .eq. 'ME') then di = (x-refx) * dskm ! Calculate the distance the point in question is from the pole dj = -re * cos(flat1)*log(cos(rflt)/(1 + sin(rflt))) + & (y - refy) * dskm ! Calculate the latitude desired in radians xlat = 2.0 * atan(exp(dj/(re*cos(flat1)))) - piovr2 ! Calculate the longitude desired in radians xlon = rfln + di/(re*cos(flat1)) ! Convert the calculated lat,lon pair into degrees xlat = xlat * 180.0/pi xlon = xlon * 180.0/pi ! If the projection is cylindrical equidistant ('CE') then ... else if (project(1:2) .eq. 'CE') then ! NOTE: Cylindrical Equidistant grid increment expressed in terms ! of thousandths of degrees! di = (x-refx) * dskm ! Calculate the distance from the horizontal axis to (J,I) dj = (y-refy) * dskm ! Determine the shift north-south xlat = reflat + dj ! Determine the shift east-west xlon = reflon + di ! If the projection is lambert conic conformal ('LC') then ... else if (project(1:2) .eq. 'LC') then isn = sign(1.0, truelat2) pi2 = piovr2*isn call lccone(truelat1,truelat2,isn,confac) tt1 = tan((pi2 - flat1)*0.5) ! Tangent Term 1. tt2 = tan((pi2 - rflt )*0.5) ! Tangent Term 2. ct1 = -cos(flat1) * re/confac ! cosine Term 1. ! Calculate the (projected) distance from the pole to ! the reference lat/lon. drp = ct1 * (tt2/tt1)**confac ! Now from the pole to the reference y along the center lon. djj = drp * cos((rcln-rfln)*confac) ! Now from the pole to the requested y along the center lon. dj = djj + isn * ((y-refy)*dskm) ! Now the (projected) distance from center longitude to reference x dii = drp*sin((rcln-rfln)*confac) ! And now from center longitude to requested X di = dii + ((x-refx)*dskm) ! Calculate the Big Messy equation bm = tt1 * (sqrt(di**2+dj**2) / abs(ct1))**(1.0/confac) ! Calculate the desired latitude in radians xlat = pi2 - 2.0*atan(bm) ! Calculate the desired longitude in radians xlon = rcln + (1.0/confac) * atan2(di,-dj) ! Convert the calculated lat,lon pair into degrees xlat = xlat * 180.0/pi xlon = xlon * 180.0/pi ! If the projection is polar stereographic ('ST') then ... else if (project(1:2) .eq. 'ST') then isn = sign(1.0, truelat1) pi2 = piovr2*isn ! Calculate the (projected) y-distance I,J lies from the pole drp = -re*cos(rflt) * (1.0 + cos(pi2-flat1)) / (1.0 + cos(pi2-rflt)) djj = drp * cos(rcln-rfln) dj = djj + isn*( (y-refy) * dskm ) dii = drp * sin(rcln-rfln) di = dii + ((x-refx)*dskm) ! Calculate the Big Messy quantity as would be done for LC ! projections. This quantity is different in value, same ! in purpose of BM above bm = (1.0/re) * sqrt(di*di + dj*dj) / (1.0 + cos(pi2-flat1)) ! Calculate the desired latitude in radians xlat = pi2 - isn*2.0 * atan(bm) ! Calculate the desired longitude in radians xlon = rcln + atan2(di,-dj) ! Convert the calculated lat,lon pair into degrees xlat = xlat * 180.0/pi xlon = xlon * 180.0/pi else print*, 'Unrecognized project: ', project stop end if ! Make sure no values are greater than 180 degrees and none ! are less than -180 degrees if (xlon .gt. 180.0) xlon = xlon - 360.0 if (xlon .lt. -180.0) xlon = xlon + 360.0 !***************************** Subroutine End ******************************! end subroutine xytoll_generic !======================================================================================= !======================================================================================= subroutine lccone ( fsplat, ssplat, sign, confac ) implicit none integer, intent(in) :: sign real, intent(in) :: fsplat real, intent(in) :: ssplat real, intent(out) :: confac if (abs(fsplat-ssplat).lt.1.E-2) then confac = sin(fsplat*degrad) else confac = log10(cos(fsplat*degrad))-log10(cos(ssplat*degrad)) confac = confac/(log10(tan((45.-float(sign)*fsplat/2.)*degrad))-& log10(tan((45.-float(sign)*ssplat/2.)*degrad))) endif end subroutine lccone end module module_llxy_generic