forked from 3rdparty/wrf-python
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
511 lines
16 KiB
511 lines
16 KiB
C NCLFORTSTART |
|
SUBROUTINE DLLTOIJ(MAP_PROJ,TRUELAT1,TRUELAT2,STDLON,LAT1,LON1, |
|
+ POLE_LAT,POLE_LON,KNOWNI,KNOWNJ,DX,DY,LATINC, |
|
+ LONINC,LAT,LON,LOC) |
|
DOUBLE PRECISION DELTALON1 |
|
DOUBLE PRECISION TL1R |
|
|
|
|
|
ccc Converts input lat/lon values to the cartesian (i,j) value |
|
ccc for the given projection. |
|
|
|
INTEGER MAP_PROJ |
|
DOUBLE PRECISION TRUELAT1,TRUELAT2,STDLON |
|
DOUBLE PRECISION LAT1,LON1,POLE_LAT,POLE_LON,KNOWNI,KNOWNJ |
|
DOUBLE PRECISION DX,DY,LATINC,LONINC,LAT,LON,LOC(2) |
|
C NCLEND |
|
|
|
DOUBLE PRECISION CLAIN,DLON,RSW,DELTALON,DELTALAT |
|
DOUBLE PRECISION REFLON,SCALE_TOP,ALA1,ALO1,ALA,ALO,RM,POLEI,POLEJ |
|
c Earth radius divided by dx |
|
DOUBLE PRECISION REBYDX |
|
DOUBLE PRECISION DELTALON1TL1R,CTL1R,ARG,CONE,HEMI |
|
DOUBLE PRECISION I,J |
|
DOUBLE PRECISION LAT1N,LON1N,OLAT,OLON |
|
|
|
DOUBLE PRECISION PI,RAD_PER_DEG,DEG_PER_RAD,RE_M |
|
|
|
ccc lat1 ! SW latitude (1,1) in degrees (-90->90N) |
|
ccc lon1 ! SW longitude (1,1) in degrees (-180->180E) |
|
ccc dx ! Grid spacing in meters at truelats |
|
ccc dlat ! Lat increment for lat/lon grids |
|
ccc dlon ! Lon increment for lat/lon grids |
|
ccc stdlon ! Longitude parallel to y-axis (-180->180E) |
|
ccc truelat1 ! First true latitude (all projections) |
|
ccc truelat2 ! Second true lat (LC only) |
|
ccc hemi ! 1 for NH, -1 for SH |
|
ccc cone ! Cone factor for LC projections |
|
ccc polei ! Computed i-location of pole point |
|
ccc polej ! Computed j-location of pole point |
|
ccc rsw ! Computed radius to SW corner |
|
ccc knowni ! X-location of known lat/lon |
|
ccc knownj ! Y-location of known lat/lon |
|
ccc RE_M ! Radius of spherical earth, meters |
|
ccc REbydx ! Earth radius divided by dx |
|
|
|
PI = 3.141592653589793D0 |
|
RAD_PER_DEG = PI/180.D0 |
|
DEG_PER_RAD = 180.D0/PI |
|
c Radius of spherical earth, meters |
|
RE_M = 6370000.D0 |
|
REBYDX = RE_M/DX |
|
|
|
HEMI = 1.0D0 |
|
IF (TRUELAT1.LT.0.0D0) THEN |
|
HEMI = -1.0D0 |
|
END IF |
|
|
|
|
|
ccc !MERCATOR |
|
IF (MAP_PROJ.EQ.3) THEN |
|
|
|
ccc ! Preliminary variables |
|
CLAIN = COS(RAD_PER_DEG*TRUELAT1) |
|
DLON = DX/ (RE_M*CLAIN) |
|
|
|
ccc ! Compute distance from equator to origin, and store in |
|
ccc ! the rsw tag. |
|
RSW = 0.D0 |
|
IF (LAT1.NE.0.D0) THEN |
|
RSW = (DLOG(TAN(0.5D0* ((LAT1+90.D0)*RAD_PER_DEG))))/DLON |
|
END IF |
|
|
|
DELTALON = LON - LON1 |
|
IF (DELTALON.LT.-180.D0) DELTALON = DELTALON + 360.D0 |
|
IF (DELTALON.GT.180.D0) DELTALON = DELTALON - 360.D0 |
|
I = KNOWNI + (DELTALON/ (DLON*DEG_PER_RAD)) |
|
J = KNOWNJ + (DLOG(TAN(0.5D0* ((LAT+90.D0)*RAD_PER_DEG))))/ |
|
+ DLON - RSW |
|
|
|
ccc !PS |
|
ELSE IF (MAP_PROJ.EQ.2) THEN |
|
|
|
REFLON = STDLON + 90.D0 |
|
|
|
ccc ! Compute numerator term of map scale factor |
|
SCALE_TOP = 1.D0 + HEMI*SIN(TRUELAT1*RAD_PER_DEG) |
|
|
|
ccc ! Compute radius to lower-left (SW) corner |
|
ALA1 = LAT1*RAD_PER_DEG |
|
RSW = REBYDX*COS(ALA1)*SCALE_TOP/ (1.D0+HEMI*SIN(ALA1)) |
|
|
|
ccc ! Find the pole point |
|
ALO1 = (LON1-REFLON)*RAD_PER_DEG |
|
POLEI = KNOWNI - RSW*COS(ALO1) |
|
POLEJ = KNOWNJ - HEMI*RSW*SIN(ALO1) |
|
|
|
ccc ! Find radius to desired point |
|
ALA = LAT*RAD_PER_DEG |
|
RM = REBYDX*COS(ALA)*SCALE_TOP/ (1.D0+HEMI*SIN(ALA)) |
|
ALO = (LON-REFLON)*RAD_PER_DEG |
|
I = POLEI + RM*COS(ALO) |
|
J = POLEJ + HEMI*RM*SIN(ALO) |
|
|
|
ccc !LAMBERT |
|
ELSE IF (MAP_PROJ.EQ.1) THEN |
|
|
|
IF (ABS(TRUELAT2).GT.90.D0) THEN |
|
TRUELAT2 = TRUELAT1 |
|
END IF |
|
|
|
IF (ABS(TRUELAT1-TRUELAT2).GT.0.1D0) THEN |
|
CONE = (DLOG(COS(TRUELAT1*RAD_PER_DEG))- |
|
+ DLOG(COS(TRUELAT2*RAD_PER_DEG)))/ |
|
+ (DLOG(TAN((90.D0-ABS(TRUELAT1))*RAD_PER_DEG* |
|
+ 0.5D0))-DLOG(TAN((90.D0-ABS(TRUELAT2))*RAD_PER_DEG* |
|
+ 0.5D0))) |
|
ELSE |
|
CONE = SIN(ABS(TRUELAT1)*RAD_PER_DEG) |
|
END IF |
|
|
|
ccc ! Compute longitude differences and ensure we stay |
|
ccc ! out of the forbidden "cut zone" |
|
DELTALON1 = LON1 - STDLON |
|
IF (DELTALON1.GT.+180.D0) DELTALON1 = DELTALON1 - 360.D0 |
|
IF (DELTALON1.LT.-180.D0) DELTALON1 = DELTALON1 + 360.D0 |
|
|
|
ccc ! Convert truelat1 to radian and compute COS for later use |
|
TL1R = TRUELAT1*RAD_PER_DEG |
|
CTL1R = COS(TL1R) |
|
|
|
ccc ! Compute the radius to our known lower-left (SW) corner |
|
RSW = REBYDX*CTL1R/CONE* (TAN((90.D0*HEMI- |
|
+ LAT1)*RAD_PER_DEG/2.D0)/TAN((90.D0*HEMI- |
|
+ TRUELAT1)*RAD_PER_DEG/2.D0))**CONE |
|
|
|
ccc ! Find pole point |
|
ARG = CONE* (DELTALON1*RAD_PER_DEG) |
|
POLEI = HEMI*KNOWNI - HEMI*RSW*SIN(ARG) |
|
POLEJ = HEMI*KNOWNJ + RSW*COS(ARG) |
|
|
|
ccc ! Compute deltalon between known longitude and standard |
|
ccc ! lon and ensure it is not in the cut zone |
|
DELTALON = LON - STDLON |
|
IF (DELTALON.GT.+180.D0) DELTALON = DELTALON - 360.D0 |
|
IF (DELTALON.LT.-180.D0) DELTALON = DELTALON + 360.D0 |
|
|
|
ccc ! Radius to desired point |
|
RM = REBYDX*CTL1R/CONE* (TAN((90.D0*HEMI- |
|
+ LAT)*RAD_PER_DEG/2.D0)/TAN((90.D0*HEMI- |
|
+ TRUELAT1)*RAD_PER_DEG/2.D0))**CONE |
|
|
|
ARG = CONE* (DELTALON*RAD_PER_DEG) |
|
I = POLEI + HEMI*RM*SIN(ARG) |
|
J = POLEJ - RM*COS(ARG) |
|
|
|
ccc ! Finally, if we are in the southern hemisphere, flip the |
|
ccc ! i/j values to a coordinate system where (1,1) is the SW |
|
ccc ! corner (what we assume) which is different than the |
|
ccc ! original NCEP algorithms which used the NE corner as |
|
ccc ! the origin in the southern hemisphere (left-hand vs. |
|
ccc ! right-hand coordinate?) |
|
I = HEMI*I |
|
J = HEMI*J |
|
|
|
|
|
ccc !lat-lon |
|
ELSE IF (MAP_PROJ.EQ.6) THEN |
|
|
|
IF (POLE_LAT.NE.90.D0) THEN |
|
CALL ROTATECOORDS(LAT,LON,OLAT,OLON,POLE_LAT,POLE_LON, |
|
+ STDLON,-1) |
|
LAT = OLAT |
|
LON = OLON + STDLON |
|
END IF |
|
|
|
c ! make sure center lat/lon is good |
|
IF (POLE_LAT.NE.90.D0) THEN |
|
CALL ROTATECOORDS(LAT1,LON1,OLAT,OLON,POLE_LAT,POLE_LON, |
|
+ STDLON,-1) |
|
LAT1N = OLAT |
|
LON1N = OLON + STDLON |
|
DELTALAT = LAT - LAT1N |
|
DELTALON = LON - LON1N |
|
ELSE |
|
DELTALAT = LAT - LAT1 |
|
DELTALON = LON - LON1 |
|
END IF |
|
|
|
c ! Compute i/j |
|
I = DELTALON/LONINC |
|
J = DELTALAT/LATINC |
|
|
|
I = I + KNOWNI |
|
J = J + KNOWNJ |
|
|
|
ELSE |
|
|
|
PRINT *,'ERROR: Do not know map projection ',MAP_PROJ |
|
|
|
END IF |
|
|
|
LOC(1) = J |
|
LOC(2) = I |
|
|
|
RETURN |
|
END |
|
|
|
|
|
C NCLFORTSTART |
|
SUBROUTINE DIJTOLL(MAP_PROJ,TRUELAT1,TRUELAT2,STDLON,LAT1,LON1, |
|
+ POLE_LAT,POLE_LON,KNOWNI,KNOWNJ,DX,DY,LATINC, |
|
+ LONINC,AI,AJ,LOC) |
|
DOUBLE PRECISION GI2 |
|
DOUBLE PRECISION ARCCOS |
|
DOUBLE PRECISION DELTALON1 |
|
DOUBLE PRECISION TL1R |
|
|
|
ccc ! Converts input lat/lon values to the cartesian (i,j) value |
|
ccc ! for the given projection. |
|
|
|
INTEGER MAP_PROJ |
|
DOUBLE PRECISION TRUELAT1,TRUELAT2,STDLON |
|
DOUBLE PRECISION LAT1,LON1,POLE_LAT,POLE_LON,KNOWNI,KNOWNJ |
|
DOUBLE PRECISION DX,DY,LATINC,LONINC,AI,AJ,LOC(2) |
|
C NCLEND |
|
|
|
DOUBLE PRECISION CLAIN,DLON,RSW,DELTALON,DELTALAT |
|
DOUBLE PRECISION REFLON,SCALE_TOP,ALA1,ALO1,ALA,ALO,RM,POLEI,POLEJ |
|
c Earth radius divided by dx |
|
DOUBLE PRECISION REBYDX |
|
DOUBLE PRECISION DELTALON1TL1R,CTL1R,ARG,CONE,HEMI |
|
|
|
DOUBLE PRECISION PI,RAD_PER_DEG,DEG_PER_RAD,RE_M |
|
|
|
DOUBLE PRECISION INEW,JNEW,R,R2 |
|
DOUBLE PRECISION CHI,CHI1,CHI2 |
|
DOUBLE PRECISION XX,YY,LAT,LON |
|
|
|
DOUBLE PRECISION RLAT,RLON,OLAT,OLON,LAT1N,LON1N |
|
DOUBLE PRECISION PHI_NP,LAM_NP,LAM_0,DLAM |
|
DOUBLE PRECISION SINPHI,COSPHI,COSLAM,SINLAM |
|
|
|
|
|
ccc lat1 ! SW latitude (1,1) in degrees (-90->90N) |
|
ccc lon1 ! SW longitude (1,1) in degrees (-180->180E) |
|
ccc dx ! Grid spacing in meters at truelats |
|
ccc dlat ! Lat increment for lat/lon grids |
|
ccc dlon ! Lon increment for lat/lon grids |
|
ccc stdlon ! Longitude parallel to y-axis (-180->180E) |
|
ccc truelat1 ! First true latitude (all projections) |
|
ccc truelat2 ! Second true lat (LC only) |
|
ccc hemi ! 1 for NH, -1 for SH |
|
ccc cone ! Cone factor for LC projections |
|
ccc polei ! Computed i-location of pole point |
|
ccc polej ! Computed j-location of pole point |
|
ccc rsw ! Computed radius to SW corner |
|
ccc knowni ! X-location of known lat/lon |
|
ccc knownj ! Y-location of known lat/lon |
|
ccc RE_M ! Radius of spherical earth, meters |
|
ccc REbydx ! Earth radius divided by dx |
|
|
|
PI = 3.141592653589793D0 |
|
RAD_PER_DEG = PI/180.D0 |
|
DEG_PER_RAD = 180.D0/PI |
|
c Radius of spherical earth, meters |
|
RE_M = 6370000.D0 |
|
REBYDX = RE_M/DX |
|
|
|
HEMI = 1.0D0 |
|
IF (TRUELAT1.LT.0.0D0) THEN |
|
HEMI = -1.0D0 |
|
END IF |
|
|
|
|
|
ccc !MERCATOR |
|
IF (MAP_PROJ.EQ.3) THEN |
|
|
|
ccc ! Preliminary variables |
|
CLAIN = COS(RAD_PER_DEG*TRUELAT1) |
|
DLON = DX/ (RE_M*CLAIN) |
|
|
|
ccc ! Compute distance from equator to origin, and store in |
|
ccc ! the rsw tag. |
|
RSW = 0.D0 |
|
IF (LAT1.NE.0.D0) THEN |
|
RSW = (DLOG(TAN(0.5D0* ((LAT1+90.D0)*RAD_PER_DEG))))/DLON |
|
END IF |
|
|
|
LAT = 2.0D0*ATAN(EXP(DLON* (RSW+AJ-KNOWNJ)))*DEG_PER_RAD - |
|
+ 90.D0 |
|
LON = (AI-KNOWNI)*DLON*DEG_PER_RAD + LON1 |
|
IF (LON.GT.180.D0) LON = LON - 360.D0 |
|
IF (LON.LT.-180.D0) LON = LON + 360.D0 |
|
|
|
|
|
ccc !PS |
|
ELSE IF (MAP_PROJ.EQ.2) THEN |
|
|
|
ccc ! Compute the reference longitude by rotating 90 degrees to |
|
ccc ! the east to find the longitude line parallel to the |
|
ccc ! positive x-axis. |
|
REFLON = STDLON + 90.D0 |
|
|
|
ccc ! Compute numerator term of map scale factor |
|
SCALE_TOP = 1.D0 + HEMI*SIN(TRUELAT1*RAD_PER_DEG) |
|
|
|
ccc ! Compute radius to known point |
|
ALA1 = LAT1*RAD_PER_DEG |
|
RSW = REBYDX*COS(ALA1)*SCALE_TOP/ (1.D0+HEMI*SIN(ALA1)) |
|
|
|
ccc ! Find the pole point |
|
ALO1 = (LON1-REFLON)*RAD_PER_DEG |
|
POLEI = KNOWNI - RSW*COS(ALO1) |
|
POLEJ = KNOWNJ - HEMI*RSW*SIN(ALO1) |
|
|
|
ccc ! Compute radius to point of interest |
|
XX = AI - POLEI |
|
YY = (AJ-POLEJ)*HEMI |
|
R2 = XX**2 + YY**2 |
|
|
|
ccc ! Now the magic code |
|
IF (R2.EQ.0.D0) THEN |
|
LAT = HEMI*90.D0 |
|
LON = REFLON |
|
ELSE |
|
GI2 = (REBYDX*SCALE_TOP)**2.D0 |
|
LAT = DEG_PER_RAD*HEMI*ASIN((GI2-R2)/ (GI2+R2)) |
|
ARCCOS = ACOS(XX/SQRT(R2)) |
|
IF (YY.GT.0) THEN |
|
LON = REFLON + DEG_PER_RAD*ARCCOS |
|
ELSE |
|
LON = REFLON - DEG_PER_RAD*ARCCOS |
|
END IF |
|
END IF |
|
|
|
ccc ! Convert to a -180 -> 180 East convention |
|
IF (LON.GT.180.D0) LON = LON - 360.D0 |
|
IF (LON.LT.-180.D0) LON = LON + 360.D0 |
|
|
|
ccc !LAMBERT |
|
ELSE IF (MAP_PROJ.EQ.1) THEN |
|
|
|
IF (ABS(TRUELAT2).GT.90.D0) THEN |
|
TRUELAT2 = TRUELAT1 |
|
END IF |
|
|
|
IF (ABS(TRUELAT1-TRUELAT2).GT.0.1D0) THEN |
|
CONE = (DLOG(COS(TRUELAT1*RAD_PER_DEG))- |
|
+ DLOG(COS(TRUELAT2*RAD_PER_DEG)))/ |
|
+ (DLOG(TAN((90.D0-ABS(TRUELAT1))*RAD_PER_DEG* |
|
+ 0.5D0))-DLOG(TAN((90.D0-ABS(TRUELAT2))*RAD_PER_DEG* |
|
+ 0.5D0))) |
|
ELSE |
|
CONE = SIN(ABS(TRUELAT1)*RAD_PER_DEG) |
|
END IF |
|
|
|
ccc ! Compute longitude differences and ensure we stay out of the |
|
ccc ! forbidden "cut zone" |
|
DELTALON1 = LON1 - STDLON |
|
IF (DELTALON1.GT.+180.D0) DELTALON1 = DELTALON1 - 360.D0 |
|
IF (DELTALON1.LT.-180.D0) DELTALON1 = DELTALON1 + 360.D0 |
|
|
|
ccc ! Convert truelat1 to radian and compute COS for later use |
|
TL1R = TRUELAT1*RAD_PER_DEG |
|
CTL1R = COS(TL1R) |
|
|
|
ccc ! Compute the radius to our known point |
|
RSW = REBYDX*CTL1R/CONE* (TAN((90.D0*HEMI- |
|
+ LAT1)*RAD_PER_DEG/2.D0)/TAN((90.D0*HEMI- |
|
+ TRUELAT1)*RAD_PER_DEG/2.D0))**CONE |
|
|
|
ccc ! Find pole point |
|
ALO1 = CONE* (DELTALON1*RAD_PER_DEG) |
|
POLEI = HEMI*KNOWNI - HEMI*RSW*SIN(ALO1) |
|
POLEJ = HEMI*KNOWNJ + RSW*COS(ALO1) |
|
|
|
CHI1 = (90.D0-HEMI*TRUELAT1)*RAD_PER_DEG |
|
CHI2 = (90.D0-HEMI*TRUELAT2)*RAD_PER_DEG |
|
|
|
ccc ! See if we are in the southern hemispere and flip the |
|
ccc ! indices if we are. |
|
INEW = HEMI*AI |
|
JNEW = HEMI*AJ |
|
|
|
ccc ! Compute radius**2 to i/j location |
|
REFLON = STDLON + 90.D0 |
|
XX = INEW - POLEI |
|
YY = POLEJ - JNEW |
|
R2 = (XX*XX+YY*YY) |
|
R = SQRT(R2)/REBYDX |
|
|
|
ccc ! Convert to lat/lon |
|
IF (R2.EQ.0.D0) THEN |
|
LAT = HEMI*90.D0 |
|
LON = STDLON |
|
ELSE |
|
LON = STDLON + DEG_PER_RAD*ATAN2(HEMI*XX,YY)/CONE |
|
LON = DMOD(LON+360.D0,360.D0) |
|
IF (CHI1.EQ.CHI2) THEN |
|
CHI = 2.0D0*ATAN((R/TAN(CHI1))** (1.D0/CONE)* |
|
+ TAN(CHI1*0.5D0)) |
|
ELSE |
|
CHI = 2.0D0*ATAN((R*CONE/SIN(CHI1))** (1.D0/CONE)* |
|
+ TAN(CHI1*0.5D0)) |
|
END IF |
|
LAT = (90.0D0-CHI*DEG_PER_RAD)*HEMI |
|
END IF |
|
|
|
IF (LON.GT.+180.D0) LON = LON - 360.D0 |
|
IF (LON.LT.-180.D0) LON = LON + 360.D0 |
|
|
|
|
|
ccc !lat-lon |
|
ELSE IF (MAP_PROJ.EQ.6) THEN |
|
|
|
INEW = AI - KNOWNI |
|
JNEW = AJ - KNOWNJ |
|
|
|
IF (INEW.LT.0.D0) INEW = INEW + 360.D0/LONINC |
|
IF (INEW.GE.360.D0/DX) INEW = INEW - 360.D0/LONINC |
|
c |
|
ccc ! Compute deltalat and deltalon |
|
DELTALAT = JNEW*LATINC |
|
DELTALON = INEW*LONINC |
|
|
|
IF (POLE_LAT.NE.90.D0) THEN |
|
CALL ROTATECOORDS(LAT1,LON1,OLAT,OLON,POLE_LAT,POLE_LON, |
|
+ STDLON,-1) |
|
LAT1N = OLAT |
|
LON1N = OLON + STDLON |
|
LAT = DELTALAT + LAT1N |
|
LON = DELTALON + LON1N |
|
ELSE |
|
LAT = DELTALAT + LAT1 |
|
LON = DELTALON + LON1 |
|
END IF |
|
|
|
|
|
IF (POLE_LAT.NE.90.D0) THEN |
|
LON = LON - STDLON |
|
CALL ROTATECOORDS(LAT,LON,OLAT,OLON,POLE_LAT,POLE_LON, |
|
+ STDLON,1) |
|
LAT = OLAT |
|
LON = OLON |
|
END IF |
|
|
|
IF (LON.LT.-180.D0) LON = LON + 360.D0 |
|
IF (LON.GT.180.D0) LON = LON - 360.D0 |
|
|
|
ELSE |
|
|
|
PRINT *,'ERROR: Do not know map projection ',MAP_PROJ |
|
|
|
END IF |
|
|
|
LOC(1) = LAT |
|
LOC(2) = LON |
|
RETURN |
|
|
|
END |
|
|
|
|
|
C NCLFORTSTART |
|
SUBROUTINE ROTATECOORDS(ILAT,ILON,OLAT,OLON,LAT_NP,LON_NP,LON_0, |
|
+ DIRECTION) |
|
DOUBLE PRECISION ILAT,ILON |
|
DOUBLE PRECISION OLAT,OLON |
|
DOUBLE PRECISION LAT_NP,LON_NP,LON_0 |
|
INTEGER DIRECTION |
|
C NCLEND |
|
|
|
c ! >=0, default : computational -> geographical |
|
c ! < 0 : geographical -> computational |
|
|
|
DOUBLE PRECISION RLAT,RLON |
|
DOUBLE PRECISION PHI_NP,LAM_NP,LAM_0,DLAM |
|
DOUBLE PRECISION SINPHI,COSPHI,COSLAM,SINLAM |
|
DOUBLE PRECISION PI,RAD_PER_DEG,DEG_PER_RAD |
|
|
|
PI = 3.141592653589793D0 |
|
RAD_PER_DEG = PI/180.D0 |
|
DEG_PER_RAD = 180.D0/PI |
|
|
|
c ! Convert all angles to radians |
|
PHI_NP = LAT_NP*RAD_PER_DEG |
|
LAM_NP = LON_NP*RAD_PER_DEG |
|
LAM_0 = LON_0*RAD_PER_DEG |
|
RLAT = ILAT*RAD_PER_DEG |
|
RLON = ILON*RAD_PER_DEG |
|
|
|
IF (DIRECTION.LT.0) THEN |
|
c ! The equations are exactly the same except for one |
|
c ! small difference with respect to longitude ... |
|
DLAM = PI - LAM_0 |
|
ELSE |
|
DLAM = LAM_NP |
|
END IF |
|
SINPHI = COS(PHI_NP)*COS(RLAT)*COS(RLON-DLAM) + |
|
+ SIN(PHI_NP)*SIN(RLAT) |
|
COSPHI = SQRT(1.D0-SINPHI*SINPHI) |
|
COSLAM = SIN(PHI_NP)*COS(RLAT)*COS(RLON-DLAM) - |
|
+ COS(PHI_NP)*SIN(RLAT) |
|
SINLAM = COS(RLAT)*SIN(RLON-DLAM) |
|
IF (COSPHI.NE.0.D0) THEN |
|
COSLAM = COSLAM/COSPHI |
|
SINLAM = SINLAM/COSPHI |
|
END IF |
|
OLAT = DEG_PER_RAD*ASIN(SINPHI) |
|
OLON = DEG_PER_RAD* (ATAN2(SINLAM,COSLAM)-DLAM-LAM_0+LAM_NP) |
|
|
|
END
|
|
|