Fortran Code dealing with longitude

! Copyright by http://scriptdemo.blogspot.com MODULE  mod_longcheck   USE mod_precision   IMPLICIT NONE CONTAINS   FUNCTION Check180(mylon)    ! check whether it is close to the International Date Line (IDL)    logical                    :: Check180    real(kind=rprec)           :: mylon    if ( abs(abs(mylon)-180.0) < 15.0) then    ! here the threshold value is arbitrary       Check180=.TRUE.    else       Check180=.FALSE.    endif   END FUNCTION Check180   FUNCTION CheckZero(mylon)    ! check if whether it is close to the Primary Line    logical                    :: CheckZero    real(kind=rprec)           :: mylon    if ( abs(mylon) < 5.0) then    ! here the threshold value is arbitrary       CheckZero=.TRUE.    else       if (abs(abs(mylon)-360.0)< 5.) then          CheckZero=.TRUE.       else          CheckZero=.FALSE.       endif    endif   END FUNCTION CheckZero   FUNCTION Back180(mylon)   ! convert a longitude into a range of [-180 180]    real(kind=rprec)   :: mylon    real(kind=rprec)   :: Back180    Back180=mylon    if (Back180<0.0) then       do while (Back180 < -180.0)          Back180=Back180+360.0       enddo    else      do while (Back180>180.0)          Back180=Back180-360.0      enddo    endif    return   END FUNCTION Back180   FUNCTION Back360(mylon)   ! convert a longitude into a range of [0 360]    real(kind=rprec)   :: mylon    real(kind=rprec)   :: Back360       Back360=mylon    if (Back360<0.0) then       do while (Back360 < 0.0)           Back360=Back360+360.0       enddo    else       do while (Back360>360.0)         Back360=Back360-360.0       enddo    endif    return   END FUNCTION Back360   SUBROUTINE ll2xy(mylon,mylat,myx,myy)    ! convert given long/latitude into x/y in stereographic projection    real(kind=rprec)   :: mylon,mylat    real(kind=rprec)   :: myx,myy    ! local variables    real(kind=rprec)               :: mypi,deg2r,myphi,mylam    real(kind=rprec)               :: rho,Re,phi0,lam0    real(kind=rprec)               :: a,b,k,f,e,e2,e4,e6,t       Re=6371000.0          ! the earth radius    mypi=acos(-1.0)    deg2r=mypi/180.0    myphi=mylat*deg2r    mylam=mylon*deg2r    phi0=0.    lam0=0.       if (mylat>84.0) then      rho=2.0*Re*(sqrt((1.0-sin(myphi))/(1.0+sin(myphi))))      myx=rho*sin(mylam-lam0)      myy=-rho*cos(mylam-lam0)      myx=myx+2000000.0      myy=myy+2000000.0    elseif (mylat<-84.0) then      ! south-hemisphere stereographic projection      myphi=-myphi      mylam=-mylam      lam0=-lam0           rho=2.0*Re*(sqrt((1.0-sin(myphi))/(1.0+sin(myphi))))      myx=rho*sin(mylam-lam0)      myy=-rho*cos(mylam-lam0)      myx=-myx; myy=-myy      myx=myx+2000000.0      myy=myy+2000000.0    else       stop "NOT DEFINED YET, maybe just not necessary"   endif   END SUBROUTINE ll2xy   SUBROUTINE xy2ll(myx,myy,mylon,mylat)    ! convert given x/y in the stereographic projection into long/latitude    real(kind=rprec)    :: myx,myy    real(kind=rprec)   :: mylon,mylat    ! local variables    real(kind=rprec)               :: mypi,deg2r,myphi,mylam,lat0    real(kind=rprec)               :: xx,yy    real(kind=rprec)               :: rho,Re,phi0,lam0    real(kind=rprec)               :: a,b,k,f,e,e2,e4,e6,t    mypi=acos(1.0)    deg2r=mypi/180.0    Re=6371000.0       lam0=0.0    phi0=0.0    lat0=90.0    if (lat0>84.0) then      xx=myx-2000000.0      yy=myy-2000000.0      rho = sqrt (xx*xx + yy*yy)        myphi=mypi*0.5-2*atan(0.5*rho/Re)      mylam=lam0+atan2(xx,-yy)    elseif (lat0<-84.0) then      stop 'Should not see this msg'      xx=-(myx-2000000.0)      yy=-(myy-2000000.0)      rho=sqrt(xx*xx+yy*yy)      myphi=0.5*mypi-2*atan(0.5*rho/Re)      myphi=-myphi      mylam=-(-lam0+atan2(xx,-yy))    else      stop 'Should not see this msg'    endif      END SUBROUTINE xy2ll END MODULE mod_longcheck