Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_llxy.F @ 3567

Last change on this file since 3567 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 73.3 KB

Rev	Line
[2759]	1	MODULE module_llxy
	2
	3	! Module that defines constants, data structures, and
	4	! subroutines used to convert grid indices to lat/lon
	5	! and vice versa.
	6	!
	7	! SUPPORTED PROJECTIONS
	8	! ---------------------
	9	! Cylindrical Lat/Lon (code = PROJ_LATLON)
	10	! Mercator (code = PROJ_MERC)
	11	! Lambert Conformal (code = PROJ_LC)
	12	! Gaussian (code = PROJ_GAUSS)
	13	! Polar Stereographic (code = PROJ_PS)
	14	! Rotated Lat/Lon (code = PROJ_ROTLL)
	15	!
	16	! REMARKS
	17	! -------
	18	! The routines contained within were adapted from routines
	19	! obtained from NCEP's w3 library. The original NCEP routines were less
	20	! flexible (e.g., polar-stereo routines only supported truelat of 60N/60S)
	21	! than what we needed, so modifications based on equations in Hoke, Hayes, and
	22	! Renninger (AFGWC/TN/79-003) were added to improve the flexibility.
	23	! Additionally, coding was improved to F90 standards and the routines were
	24	! combined into this module.
	25	!
	26	! ASSUMPTIONS
	27	! -----------
	28	! Grid Definition:
	29	! For mercator, lambert conformal, and polar-stereographic projections,
	30	! the routines within assume the following:
	31	!
	32	! 1. Grid is dimensioned (i,j) where i is the East-West direction,
	33	! positive toward the east, and j is the north-south direction,
	34	! positive toward the north.
	35	! 2. Origin is at (1,1) and is located at the southwest corner,
	36	! regardless of hemispere.
	37	! 3. Grid spacing (dx) is always positive.
	38	! 4. Values of true latitudes must be positive for NH domains
	39	! and negative for SH domains.
	40	!
	41	! For the latlon and Gaussian projection, the grid origin may be at any
	42	! of the corners, and the deltalat and deltalon values can be signed to
	43	! account for this using the following convention:
	44	! Origin Location Deltalat Sign Deltalon Sign
	45	! --------------- ------------- -------------
	46	! SW Corner + +
	47	! NE Corner - -
	48	! NW Corner - +
	49	! SE Corner + -
	50	!
	51	! Data Definitions:
	52	! 1. Any arguments that are a latitude value are expressed in
	53	! degrees north with a valid range of -90 -> 90
	54	! 2. Any arguments that are a longitude value are expressed in
	55	! degrees east with a valid range of -180 -> 180.
	56	! 3. Distances are in meters and are always positive.
	57	! 4. The standard longitude (stdlon) is defined as the longitude
	58	! line which is parallel to the grid's y-axis (j-direction), along
	59	! which latitude increases (NOT the absolute value of latitude, but
	60	! the actual latitude, such that latitude increases continuously
	61	! from the south pole to the north pole) as j increases.
	62	! 5. One true latitude value is required for polar-stereographic and
	63	! mercator projections, and defines at which latitude the
	64	! grid spacing is true. For lambert conformal, two true latitude
	65	! values must be specified, but may be set equal to each other to
	66	! specify a tangent projection instead of a secant projection.
	67	!
	68	! USAGE
	69	! -----
	70	! To use the routines in this module, the calling routines must have the
	71	! following statement at the beginning of its declaration block:
	72	! USE map_utils
	73	!
	74	! The use of the module not only provides access to the necessary routines,
	75	! but also defines a structure of TYPE (proj_info) that can be used
	76	! to declare a variable of the same type to hold your map projection
	77	! information. It also defines some integer parameters that contain
	78	! the projection codes so one only has to use those variable names rather
	79	! than remembering the acutal code when using them. The basic steps are
	80	! as follows:
	81	!
	82	! 1. Ensure the "USE map_utils" is in your declarations.
	83	! 2. Declare the projection information structure as type(proj_info):
	84	! TYPE(proj_info) :: proj
	85	! 3. Populate your structure by calling the map_set routine:
	86	! CALL map_set(code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,proj)
	87	! where:
	88	! code (input) = one of PROJ_LATLON, PROJ_MERC, PROJ_LC, PROJ_PS,
	89	! PROJ_GAUSS, or PROJ_ROTLL
	90	! lat1 (input) = Latitude of grid origin point (i,j)=(1,1)
	91	! (see assumptions!)
	92	! lon1 (input) = Longitude of grid origin
	93	! knowni (input) = origin point, x-location
	94	! knownj (input) = origin point, y-location
	95	! dx (input) = grid spacing in meters (ignored for LATLON projections)
	96	! stdlon (input) = Standard longitude for PROJ_PS and PROJ_LC,
	97	! deltalon (see assumptions) for PROJ_LATLON,
	98	! ignored for PROJ_MERC
	99	! truelat1 (input) = 1st true latitude for PROJ_PS, PROJ_LC, and
	100	! PROJ_MERC, deltalat (see assumptions) for PROJ_LATLON
	101	! truelat2 (input) = 2nd true latitude for PROJ_LC,
	102	! ignored for all others.
	103	! proj (output) = The structure of type (proj_info) that will be fully
	104	! populated after this call
	105	!
	106	! 4. Now that the proj structure is populated, you may call either
	107	! of the following routines:
	108	!
	109	! latlon_to_ij(proj, lat, lon, i, j)
	110	! ij_to_latlon(proj, i, j, lat, lon)
	111	!
	112	! It is incumbent upon the calling routine to determine whether or
	113	! not the values returned are within your domain's bounds. All values
	114	! of i, j, lat, and lon are REAL values.
	115	!
	116	!
	117	! REFERENCES
	118	! ----------
	119	! Hoke, Hayes, and Renninger, "Map Preojections and Grid Systems for
	120	! Meteorological Applications." AFGWC/TN-79/003(Rev), Air Weather
	121	! Service, 1985.
	122	!
	123	! NCAR MM5v3 Modeling System, REGRIDDER program, module_first_guess_map.F
	124	! NCEP routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
	125	!
	126	! HISTORY
	127	! -------
	128	! 27 Mar 2001 - Original Version
	129	! Brent L. Shaw, NOAA/FSL (CSU/CIRA)
	130	!
	131	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	132
	133	USE module_wrf_error
	134
	135	INTEGER, PARAMETER :: HH=4, VV=5
	136
	137	REAL, PARAMETER :: PI = 3.141592653589793
	138	REAL, PARAMETER :: OMEGA_E = 7.292e-5
	139
	140	REAL, PARAMETER :: DEG_PER_RAD = 180./PI
	141	REAL, PARAMETER :: RAD_PER_DEG = PI/180.
	142
	143	REAL, PARAMETER :: A_WGS84 = 6378137.
	144	REAL, PARAMETER :: B_WGS84 = 6356752.314
	145	REAL, PARAMETER :: RE_WGS84 = A_WGS84
	146	REAL, PARAMETER :: E_WGS84 = 0.081819192
	147
	148	REAL, PARAMETER :: A_NAD83 = 6378137.
	149	REAL, PARAMETER :: RE_NAD83 = A_NAD83
	150	REAL, PARAMETER :: E_NAD83 = 0.0818187034
	151
	152	REAL, PARAMETER :: EARTH_RADIUS_M = 6370000.
	153	REAL, PARAMETER :: EARTH_CIRC_M = 2.PIEARTH_RADIUS_M
	154
	155	INTEGER, PUBLIC, PARAMETER :: PROJ_LATLON = 0
	156	INTEGER, PUBLIC, PARAMETER :: PROJ_LC = 1
	157	INTEGER, PUBLIC, PARAMETER :: PROJ_PS = 2
	158	INTEGER, PUBLIC, PARAMETER :: PROJ_PS_WGS84 = 102
	159	INTEGER, PUBLIC, PARAMETER :: PROJ_MERC = 3
	160	INTEGER, PUBLIC, PARAMETER :: PROJ_GAUSS = 4
	161	INTEGER, PUBLIC, PARAMETER :: PROJ_CYL = 5
	162	INTEGER, PUBLIC, PARAMETER :: PROJ_CASSINI = 6
	163	INTEGER, PUBLIC, PARAMETER :: PROJ_ALBERS_NAD83 = 105
	164	INTEGER, PUBLIC, PARAMETER :: PROJ_ROTLL = 203
	165
	166	! Define some private constants
	167	INTEGER, PRIVATE, PARAMETER :: HIGH = 8
	168
	169	TYPE proj_info
	170
	171	INTEGER :: code ! Integer code for projection TYPE
	172	INTEGER :: nlat ! For Gaussian -- number of latitude points
	173	! north of the equator
	174	INTEGER :: nlon !
	175	!
	176	INTEGER :: ixdim ! For Rotated Lat/Lon -- number of mass points
	177	! in an odd row
	178	INTEGER :: jydim ! For Rotated Lat/Lon -- number of rows
	179	INTEGER :: stagger ! For Rotated Lat/Lon -- mass or velocity grid
	180	REAL :: phi ! For Rotated Lat/Lon -- domain half-extent in
	181	! degrees latitude
	182	REAL :: lambda ! For Rotated Lat/Lon -- domain half-extend in
	183	! degrees longitude
	184	REAL :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
	185	REAL :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
	186	REAL :: lat0 ! For Cassini, latitude of projection pole
	187	REAL :: lon0 ! For Cassini, longitude of projection pole
	188	REAL :: dx ! Grid spacing in meters at truelats, used
	189	! only for ps, lc, and merc projections
	190	REAL :: dy ! Grid spacing in meters at truelats, used
	191	! only for ps, lc, and merc projections
	192	REAL :: latinc ! Latitude increment for cylindrical lat/lon
	193	REAL :: loninc ! Longitude increment for cylindrical lat/lon
	194	! also the lon increment for Gaussian grid
	195	REAL :: dlat ! Lat increment for lat/lon grids
	196	REAL :: dlon ! Lon increment for lat/lon grids
	197	REAL :: stdlon ! Longitude parallel to y-axis (-180->180E)
	198	REAL :: truelat1 ! First true latitude (all projections)
	199	REAL :: truelat2 ! Second true lat (LC only)
	200	REAL :: hemi ! 1 for NH, -1 for SH
	201	REAL :: cone ! Cone factor for LC projections
	202	REAL :: polei ! Computed i-location of pole point
	203	REAL :: polej ! Computed j-location of pole point
	204	REAL :: rsw ! Computed radius to SW corner
	205	REAL :: rebydx ! Earth radius divided by dx
	206	REAL :: knowni ! X-location of known lat/lon
	207	REAL :: knownj ! Y-location of known lat/lon
	208	REAL :: re_m ! Radius of spherical earth, meters
	209	REAL :: rho0 ! For Albers equal area
	210	REAL :: nc ! For Albers equal area
	211	REAL :: bigc ! For Albers equal area
	212	LOGICAL :: init ! Flag to indicate if this struct is
	213	! ready for use
	214	LOGICAL :: wrap ! For Gaussian -- flag to indicate wrapping
	215	! around globe?
	216	REAL, POINTER, DIMENSION(:) :: gauss_lat ! Latitude array for Gaussian grid
	217
	218	END TYPE proj_info
	219
	220	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	221	CONTAINS
	222	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	223
	224	SUBROUTINE map_init(proj)
	225	! Initializes the map projection structure to missing values
	226
	227	IMPLICIT NONE
	228	TYPE(proj_info), INTENT(INOUT) :: proj
	229
	230	proj%lat1 = -999.9
	231	proj%lon1 = -999.9
	232	proj%lat0 = -999.9
	233	proj%lon0 = -999.9
	234	proj%dx = -999.9
	235	proj%dy = -999.9
	236	proj%latinc = -999.9
	237	proj%loninc = -999.9
	238	proj%stdlon = -999.9
	239	proj%truelat1 = -999.9
	240	proj%truelat2 = -999.9
	241	proj%phi = -999.9
	242	proj%lambda = -999.9
	243	proj%ixdim = -999
	244	proj%jydim = -999
	245	proj%stagger = HH
	246	proj%nlat = 0
	247	proj%nlon = 0
	248	proj%hemi = 0.0
	249	proj%cone = -999.9
	250	proj%polei = -999.9
	251	proj%polej = -999.9
	252	proj%rsw = -999.9
	253	proj%knowni = -999.9
	254	proj%knownj = -999.9
	255	proj%re_m = EARTH_RADIUS_M
	256	proj%init = .FALSE.
	257	proj%wrap = .FALSE.
	258	proj%rho0 = 0.
	259	proj%nc = 0.
	260	proj%bigc = 0.
	261	nullify(proj%gauss_lat)
	262
	263	END SUBROUTINE map_init
	264
	265
	266	SUBROUTINE map_set(proj_code, proj, lat1, lon1, lat0, lon0, knowni, knownj, dx, latinc, &
	267	loninc, stdlon, truelat1, truelat2, nlat, nlon, ixdim, jydim, &
	268	stagger, phi, lambda, r_earth)
	269	! Given a partially filled proj_info structure, this routine computes
	270	! polei, polej, rsw, and cone (if LC projection) to complete the
	271	! structure. This allows us to eliminate redundant calculations when
	272	! calling the coordinate conversion routines multiple times for the
	273	! same map.
	274	! This will generally be the first routine called when a user wants
	275	! to be able to use the coordinate conversion routines, and it
	276	! will call the appropriate subroutines based on the
	277	! proj%code which indicates which projection type this is.
	278
	279	IMPLICIT NONE
	280
	281	! Declare arguments
	282	INTEGER, INTENT(IN) :: proj_code
	283	INTEGER, INTENT(IN), OPTIONAL :: nlat
	284	INTEGER, INTENT(IN), OPTIONAL :: nlon
	285	INTEGER, INTENT(IN), OPTIONAL :: ixdim
	286	INTEGER, INTENT(IN), OPTIONAL :: jydim
	287	INTEGER, INTENT(IN), OPTIONAL :: stagger
	288	REAL, INTENT(IN), OPTIONAL :: latinc
	289	REAL, INTENT(IN), OPTIONAL :: loninc
	290	REAL, INTENT(IN), OPTIONAL :: lat1
	291	REAL, INTENT(IN), OPTIONAL :: lon1
	292	REAL, INTENT(IN), OPTIONAL :: lat0
	293	REAL, INTENT(IN), OPTIONAL :: lon0
	294	REAL, INTENT(IN), OPTIONAL :: dx
	295	REAL, INTENT(IN), OPTIONAL :: stdlon
	296	REAL, INTENT(IN), OPTIONAL :: truelat1
	297	REAL, INTENT(IN), OPTIONAL :: truelat2
	298	REAL, INTENT(IN), OPTIONAL :: knowni
	299	REAL, INTENT(IN), OPTIONAL :: knownj
	300	REAL, INTENT(IN), OPTIONAL :: phi
	301	REAL, INTENT(IN), OPTIONAL :: lambda
	302	REAL, INTENT(IN), OPTIONAL :: r_earth
	303	TYPE(proj_info), INTENT(OUT) :: proj
	304
	305	INTEGER :: iter
	306	REAL :: dummy_lon1
	307	REAL :: dummy_lon0
	308	REAL :: dummy_stdlon
	309
	310	! First, verify that mandatory parameters are present for the specified proj_code
	311	IF ( proj_code == PROJ_LC ) THEN
	312	IF ( .NOT.PRESENT(truelat1) .OR. &
	313	.NOT.PRESENT(truelat2) .OR. &
	314	.NOT.PRESENT(lat1) .OR. &
	315	.NOT.PRESENT(lon1) .OR. &
	316	.NOT.PRESENT(knowni) .OR. &
	317	.NOT.PRESENT(knownj) .OR. &
	318	.NOT.PRESENT(stdlon) .OR. &
	319	.NOT.PRESENT(dx) ) THEN
	320	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	321	PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knowni, knownj, stdlon, dx'
	322	CALL wrf_error_fatal ( 'MAP_INIT' )
	323	END IF
	324	ELSE IF ( proj_code == PROJ_PS ) THEN
	325	IF ( .NOT.PRESENT(truelat1) .OR. &
	326	.NOT.PRESENT(lat1) .OR. &
	327	.NOT.PRESENT(lon1) .OR. &
	328	.NOT.PRESENT(knowni) .OR. &
	329	.NOT.PRESENT(knownj) .OR. &
	330	.NOT.PRESENT(stdlon) .OR. &
	331	.NOT.PRESENT(dx) ) THEN
	332	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	333	PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx'
	334	CALL wrf_error_fatal ( 'MAP_INIT' )
	335	END IF
	336	ELSE IF ( proj_code == PROJ_PS_WGS84 ) THEN
	337	IF ( .NOT.PRESENT(truelat1) .OR. &
	338	.NOT.PRESENT(lat1) .OR. &
	339	.NOT.PRESENT(lon1) .OR. &
	340	.NOT.PRESENT(knowni) .OR. &
	341	.NOT.PRESENT(knownj) .OR. &
	342	.NOT.PRESENT(stdlon) .OR. &
	343	.NOT.PRESENT(dx) ) THEN
	344	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	345	PRINT '(A)', ' truelat1, lat1, lon1, knonwi, knownj, stdlon, dx'
	346	CALL wrf_error_fatal ( 'MAP_INIT' )
	347	END IF
	348	ELSE IF ( proj_code == PROJ_ALBERS_NAD83 ) THEN
	349	IF ( .NOT.PRESENT(truelat1) .OR. &
	350	.NOT.PRESENT(truelat2) .OR. &
	351	.NOT.PRESENT(lat1) .OR. &
	352	.NOT.PRESENT(lon1) .OR. &
	353	.NOT.PRESENT(knowni) .OR. &
	354	.NOT.PRESENT(knownj) .OR. &
	355	.NOT.PRESENT(stdlon) .OR. &
	356	.NOT.PRESENT(dx) ) THEN
	357	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	358	PRINT '(A)', ' truelat1, truelat2, lat1, lon1, knonwi, knownj, stdlon, dx'
	359	CALL wrf_error_fatal ( 'MAP_INIT' )
	360	END IF
	361	ELSE IF ( proj_code == PROJ_MERC ) THEN
	362	IF ( .NOT.PRESENT(truelat1) .OR. &
	363	.NOT.PRESENT(lat1) .OR. &
	364	.NOT.PRESENT(lon1) .OR. &
	365	.NOT.PRESENT(knowni) .OR. &
	366	.NOT.PRESENT(knownj) .OR. &
	367	.NOT.PRESENT(dx) ) THEN
	368	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	369	PRINT '(A)', ' truelat1, lat1, lon1, knowni, knownj, dx'
	370	CALL wrf_error_fatal ( 'MAP_INIT' )
	371	END IF
	372	ELSE IF ( proj_code == PROJ_LATLON ) THEN
	373	IF ( .NOT.PRESENT(latinc) .OR. &
	374	.NOT.PRESENT(loninc) .OR. &
	375	.NOT.PRESENT(knowni) .OR. &
	376	.NOT.PRESENT(knownj) .OR. &
	377	.NOT.PRESENT(lat1) .OR. &
	378	.NOT.PRESENT(lon1) ) THEN
	379	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	380	PRINT '(A)', ' latinc, loninc, knowni, knownj, lat1, lon1'
	381	CALL wrf_error_fatal ( 'MAP_INIT' )
	382	END IF
	383	ELSE IF ( proj_code == PROJ_CYL ) THEN
	384	IF ( .NOT.PRESENT(latinc) .OR. &
	385	.NOT.PRESENT(loninc) .OR. &
	386	.NOT.PRESENT(stdlon) ) THEN
	387	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	388	PRINT '(A)', ' latinc, loninc, stdlon'
	389	CALL wrf_error_fatal ( 'MAP_INIT' )
	390	END IF
	391	ELSE IF ( proj_code == PROJ_CASSINI ) THEN
	392	IF ( .NOT.PRESENT(latinc) .OR. &
	393	.NOT.PRESENT(loninc) .OR. &
	394	.NOT.PRESENT(lat1) .OR. &
	395	.NOT.PRESENT(lon1) .OR. &
	396	.NOT.PRESENT(lat0) .OR. &
	397	.NOT.PRESENT(lon0) .OR. &
	398	.NOT.PRESENT(knowni) .OR. &
	399	.NOT.PRESENT(knownj) .OR. &
	400	.NOT.PRESENT(stdlon) ) THEN
	401	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	402	PRINT '(A)', ' latinc, loninc, lat1, lon1, knowni, knownj, lat0, lon0, stdlon'
	403	CALL wrf_error_fatal ( 'MAP_INIT' )
	404	END IF
	405	ELSE IF ( proj_code == PROJ_GAUSS ) THEN
	406	IF ( .NOT.PRESENT(nlat) .OR. &
	407	.NOT.PRESENT(lat1) .OR. &
	408	.NOT.PRESENT(lon1) .OR. &
	409	.NOT.PRESENT(loninc) ) THEN
	410	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	411	PRINT '(A)', ' nlat, lat1, lon1, loninc'
	412	CALL wrf_error_fatal ( 'MAP_INIT' )
	413	END IF
	414	ELSE IF ( proj_code == PROJ_ROTLL ) THEN
	415	IF ( .NOT.PRESENT(ixdim) .OR. &
	416	.NOT.PRESENT(jydim) .OR. &
	417	.NOT.PRESENT(phi) .OR. &
	418	.NOT.PRESENT(lambda) .OR. &
	419	.NOT.PRESENT(lat1) .OR. &
	420	.NOT.PRESENT(lon1) .OR. &
	421	.NOT.PRESENT(stagger) ) THEN
	422	PRINT '(A,I2)', 'The following are mandatory parameters for projection code : ', proj_code
	423	PRINT '(A)', ' ixdim, jydim, phi, lambda, lat1, lon1, stagger'
	424	CALL wrf_error_fatal ( 'MAP_INIT' )
	425	END IF
	426	ELSE
	427	PRINT '(A,I2)', 'Unknown projection code: ', proj_code
	428	CALL wrf_error_fatal ( 'MAP_INIT' )
	429	END IF
	430
	431	! Check for validity of mandatory variables in proj
	432	IF ( PRESENT(lat1) ) THEN
	433	IF ( ABS(lat1) .GT. 90. ) THEN
	434	PRINT '(A)', 'Latitude of origin corner required as follows:'
	435	PRINT '(A)', ' -90N <= lat1 < = 90.N'
	436	CALL wrf_error_fatal ( 'MAP_INIT' )
	437	ENDIF
	438	ENDIF
	439
	440	IF ( PRESENT(lon1) ) THEN
	441	dummy_lon1 = lon1
	442	IF ( ABS(dummy_lon1) .GT. 180.) THEN
	443	iter = 0
	444	DO WHILE (ABS(dummy_lon1) > 180. .AND. iter < 10)
	445	IF (dummy_lon1 < -180.) dummy_lon1 = dummy_lon1 + 360.
	446	IF (dummy_lon1 > 180.) dummy_lon1 = dummy_lon1 - 360.
	447	iter = iter + 1
	448	END DO
	449	IF (abs(dummy_lon1) > 180.) THEN
	450	PRINT '(A)', 'Longitude of origin required as follows:'
	451	PRINT '(A)', ' -180E <= lon1 <= 180W'
	452	CALL wrf_error_fatal ( 'MAP_INIT' )
	453	ENDIF
	454	ENDIF
	455	ENDIF
	456
	457	IF ( PRESENT(lon0) ) THEN
	458	dummy_lon0 = lon0
	459	IF ( ABS(dummy_lon0) .GT. 180.) THEN
	460	iter = 0
	461	DO WHILE (ABS(dummy_lon0) > 180. .AND. iter < 10)
	462	IF (dummy_lon0 < -180.) dummy_lon0 = dummy_lon0 + 360.
	463	IF (dummy_lon0 > 180.) dummy_lon0 = dummy_lon0 - 360.
	464	iter = iter + 1
	465	END DO
	466	IF (abs(dummy_lon0) > 180.) THEN
	467	PRINT '(A)', 'Longitude of pole required as follows:'
	468	PRINT '(A)', ' -180E <= lon0 <= 180W'
	469	CALL wrf_error_fatal ( 'MAP_INIT' )
	470	ENDIF
	471	ENDIF
	472	ENDIF
	473
	474	IF ( PRESENT(dx) ) THEN
	475	IF ((dx .LE. 0.).AND.(proj_code .NE. PROJ_LATLON)) THEN
	476	PRINT '(A)', 'Require grid spacing (dx) in meters be positive'
	477	CALL wrf_error_fatal ( 'MAP_INIT' )
	478	ENDIF
	479	ENDIF
	480
	481	IF ( PRESENT(stdlon) ) THEN
	482	dummy_stdlon = stdlon
	483	IF ((ABS(dummy_stdlon) > 180.).AND.(proj_code /= PROJ_MERC)) THEN
	484	iter = 0
	485	DO WHILE (ABS(dummy_stdlon) > 180. .AND. iter < 10)
	486	IF (dummy_stdlon < -180.) dummy_stdlon = dummy_stdlon + 360.
	487	IF (dummy_stdlon > 180.) dummy_stdlon = dummy_stdlon - 360.
	488	iter = iter + 1
	489	END DO
	490	IF (abs(dummy_stdlon) > 180.) THEN
	491	PRINT '(A)', 'Need orientation longitude (stdlon) as: '
	492	PRINT '(A)', ' -180E <= stdlon <= 180W'
	493	CALL wrf_error_fatal ( 'MAP_INIT' )
	494	ENDIF
	495	ENDIF
	496	ENDIF
	497
	498	IF ( PRESENT(truelat1) ) THEN
	499	IF (ABS(truelat1).GT.90.) THEN
	500	PRINT '(A)', 'Set true latitude 1 for all projections'
	501	CALL wrf_error_fatal ( 'MAP_INIT' )
	502	ENDIF
	503	ENDIF
	504
	505	CALL map_init(proj)
	506	proj%code = proj_code
	507	IF ( PRESENT(lat1) ) proj%lat1 = lat1
	508	IF ( PRESENT(lon1) ) proj%lon1 = dummy_lon1
	509	IF ( PRESENT(lat0) ) proj%lat0 = lat0
	510	IF ( PRESENT(lon0) ) proj%lon0 = dummy_lon0
	511	IF ( PRESENT(latinc) ) proj%latinc = latinc
	512	IF ( PRESENT(loninc) ) proj%loninc = loninc
	513	IF ( PRESENT(knowni) ) proj%knowni = knowni
	514	IF ( PRESENT(knownj) ) proj%knownj = knownj
	515	IF ( PRESENT(dx) ) proj%dx = dx
	516	IF ( PRESENT(stdlon) ) proj%stdlon = dummy_stdlon
	517	IF ( PRESENT(truelat1) ) proj%truelat1 = truelat1
	518	IF ( PRESENT(truelat2) ) proj%truelat2 = truelat2
	519	IF ( PRESENT(nlat) ) proj%nlat = nlat
	520	IF ( PRESENT(nlon) ) proj%nlon = nlon
	521	IF ( PRESENT(ixdim) ) proj%ixdim = ixdim
	522	IF ( PRESENT(jydim) ) proj%jydim = jydim
	523	IF ( PRESENT(stagger) ) proj%stagger = stagger
	524	IF ( PRESENT(phi) ) proj%phi = phi
	525	IF ( PRESENT(lambda) ) proj%lambda = lambda
	526	IF ( PRESENT(r_earth) ) proj%re_m = r_earth
	527
	528	IF ( PRESENT(dx) ) THEN
	529	IF ( (proj_code == PROJ_LC) .OR. (proj_code == PROJ_PS) .OR. &
	530	(proj_code == PROJ_PS_WGS84) .OR. (proj_code == PROJ_ALBERS_NAD83) .OR. &
	531	(proj_code == PROJ_MERC) ) THEN
	532	proj%dx = dx
	533	IF (truelat1 .LT. 0.) THEN
	534	proj%hemi = -1.0
	535	ELSE
	536	proj%hemi = 1.0
	537	ENDIF
	538	proj%rebydx = proj%re_m / dx
	539	ENDIF
	540	ENDIF
	541
	542	pick_proj: SELECT CASE(proj%code)
	543
	544	CASE(PROJ_PS)
	545	CALL set_ps(proj)
	546
	547	CASE(PROJ_PS_WGS84)
	548	CALL set_ps_wgs84(proj)
	549
	550	CASE(PROJ_ALBERS_NAD83)
	551	CALL set_albers_nad83(proj)
	552
	553	CASE(PROJ_LC)
	554	IF (ABS(proj%truelat2) .GT. 90.) THEN
	555	proj%truelat2=proj%truelat1
	556	ENDIF
	557	CALL set_lc(proj)
	558
	559	CASE (PROJ_MERC)
	560	CALL set_merc(proj)
	561
	562	CASE (PROJ_LATLON)
	563
	564	CASE (PROJ_GAUSS)
	565	CALL set_gauss(proj)
	566
	567	CASE (PROJ_CYL)
	568	CALL set_cyl(proj)
	569
	570	CASE (PROJ_CASSINI)
	571	CALL set_cassini(proj)
	572
	573	CASE (PROJ_ROTLL)
	574
	575	END SELECT pick_proj
	576	proj%init = .TRUE.
	577
	578	RETURN
	579
	580	END SUBROUTINE map_set
	581
	582
	583	SUBROUTINE latlon_to_ij(proj, lat, lon, i, j)
	584	! Converts input lat/lon values to the cartesian (i,j) value
	585	! for the given projection.
	586
	587	IMPLICIT NONE
	588	TYPE(proj_info), INTENT(IN) :: proj
	589	REAL, INTENT(IN) :: lat
	590	REAL, INTENT(IN) :: lon
	591	REAL, INTENT(OUT) :: i
	592	REAL, INTENT(OUT) :: j
	593
	594	IF (.NOT.proj%init) THEN
	595	PRINT '(A)', 'You have not called map_set for this projection'
	596	CALL wrf_error_fatal ( 'LATLON_TO_IJ' )
	597	ENDIF
	598
	599	SELECT CASE(proj%code)
	600
	601	CASE(PROJ_LATLON)
	602	CALL llij_latlon(lat,lon,proj,i,j)
	603
	604	CASE(PROJ_MERC)
	605	CALL llij_merc(lat,lon,proj,i,j)
	606
	607	CASE(PROJ_PS)
	608	CALL llij_ps(lat,lon,proj,i,j)
	609
	610	CASE(PROJ_PS_WGS84)
	611	CALL llij_ps_wgs84(lat,lon,proj,i,j)
	612
	613	CASE(PROJ_ALBERS_NAD83)
	614	CALL llij_albers_nad83(lat,lon,proj,i,j)
	615
	616	CASE(PROJ_LC)
	617	CALL llij_lc(lat,lon,proj,i,j)
	618
	619	CASE(PROJ_GAUSS)
	620	CALL llij_gauss(lat,lon,proj,i,j)
	621
	622	CASE(PROJ_CYL)
	623	CALL llij_cyl(lat,lon,proj,i,j)
	624
	625	CASE(PROJ_CASSINI)
	626	CALL llij_cassini(lat,lon,proj,i,j)
	627
	628	CASE(PROJ_ROTLL)
	629	CALL llij_rotlatlon(lat,lon,proj,i,j)
	630
	631	CASE DEFAULT
	632	PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
	633	CALL wrf_error_fatal ( 'LATLON_TO_IJ' )
	634
	635	END SELECT
	636
	637	RETURN
	638
	639	END SUBROUTINE latlon_to_ij
	640
	641
	642	SUBROUTINE ij_to_latlon(proj, i, j, lat, lon)
	643	! Computes geographical latitude and longitude for a given (i,j) point
	644	! in a grid with a projection of proj
	645
	646	IMPLICIT NONE
	647	TYPE(proj_info),INTENT(IN) :: proj
	648	REAL, INTENT(IN) :: i
	649	REAL, INTENT(IN) :: j
	650	REAL, INTENT(OUT) :: lat
	651	REAL, INTENT(OUT) :: lon
	652
	653	IF (.NOT.proj%init) THEN
	654	PRINT '(A)', 'You have not called map_set for this projection'
	655	CALL wrf_error_fatal ( 'IJ_TO_LATLON' )
	656	ENDIF
	657	SELECT CASE (proj%code)
	658
	659	CASE (PROJ_LATLON)
	660	CALL ijll_latlon(i, j, proj, lat, lon)
	661
	662	CASE (PROJ_MERC)
	663	CALL ijll_merc(i, j, proj, lat, lon)
	664
	665	CASE (PROJ_PS)
	666	CALL ijll_ps(i, j, proj, lat, lon)
	667
	668	CASE (PROJ_PS_WGS84)
	669	CALL ijll_ps_wgs84(i, j, proj, lat, lon)
	670
	671	CASE (PROJ_ALBERS_NAD83)
	672	CALL ijll_albers_nad83(i, j, proj, lat, lon)
	673
	674	CASE (PROJ_LC)
	675	CALL ijll_lc(i, j, proj, lat, lon)
	676
	677	CASE (PROJ_CYL)
	678	CALL ijll_cyl(i, j, proj, lat, lon)
	679
	680	CASE (PROJ_CASSINI)
	681	CALL ijll_cassini(i, j, proj, lat, lon)
	682
	683	CASE (PROJ_ROTLL)
	684	CALL ijll_rotlatlon(i, j, proj, lat, lon)
	685
	686	CASE DEFAULT
	687	PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
	688	CALL wrf_error_fatal ( 'IJ_TO_LATLON' )
	689
	690	END SELECT
	691	RETURN
	692	END SUBROUTINE ij_to_latlon
	693
	694
	695	SUBROUTINE set_ps(proj)
	696	! Initializes a polar-stereographic map projection from the partially
	697	! filled proj structure. This routine computes the radius to the
	698	! southwest corner and computes the i/j location of the pole for use
	699	! in llij_ps and ijll_ps.
	700	IMPLICIT NONE
	701
	702	! Declare args
	703	TYPE(proj_info), INTENT(INOUT) :: proj
	704
	705	! Local vars
	706	REAL :: ala1
	707	REAL :: alo1
	708	REAL :: reflon
	709	REAL :: scale_top
	710
	711	! Executable code
	712	reflon = proj%stdlon + 90.
	713
	714	! Compute numerator term of map scale factor
	715	scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
	716
	717	! Compute radius to lower-left (SW) corner
	718	ala1 = proj%lat1 * rad_per_deg
	719	proj%rsw = proj%rebydxCOS(ala1)scale_top/(1.+proj%hemi*SIN(ala1))
	720
	721	! Find the pole point
	722	alo1 = (proj%lon1 - reflon) * rad_per_deg
	723	proj%polei = proj%knowni - proj%rsw * COS(alo1)
	724	proj%polej = proj%knownj - proj%hemi * proj%rsw * SIN(alo1)
	725
	726	RETURN
	727
	728	END SUBROUTINE set_ps
	729
	730
	731	SUBROUTINE llij_ps(lat,lon,proj,i,j)
	732	! Given latitude (-90 to 90), longitude (-180 to 180), and the
	733	! standard polar-stereographic projection information via the
	734	! public proj structure, this routine returns the i/j indices which
	735	! if within the domain range from 1->nx and 1->ny, respectively.
	736
	737	IMPLICIT NONE
	738
	739	! Delcare input arguments
	740	REAL, INTENT(IN) :: lat
	741	REAL, INTENT(IN) :: lon
	742	TYPE(proj_info),INTENT(IN) :: proj
	743
	744	! Declare output arguments
	745	REAL, INTENT(OUT) :: i !(x-index)
	746	REAL, INTENT(OUT) :: j !(y-index)
	747
	748	! Declare local variables
	749
	750	REAL :: reflon
	751	REAL :: scale_top
	752	REAL :: ala
	753	REAL :: alo
	754	REAL :: rm
	755
	756	! BEGIN CODE
	757
	758	reflon = proj%stdlon + 90.
	759
	760	! Compute numerator term of map scale factor
	761
	762	scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
	763
	764	! Find radius to desired point
	765	ala = lat * rad_per_deg
	766	rm = proj%rebydx * COS(ala) * scale_top/(1. + proj%hemi *SIN(ala))
	767	alo = (lon - reflon) * rad_per_deg
	768	i = proj%polei + rm * COS(alo)
	769	j = proj%polej + proj%hemi * rm * SIN(alo)
	770
	771	RETURN
	772
	773	END SUBROUTINE llij_ps
	774
	775
	776	SUBROUTINE ijll_ps(i, j, proj, lat, lon)
	777
	778	! This is the inverse subroutine of llij_ps. It returns the
	779	! latitude and longitude of an i/j point given the projection info
	780	! structure.
	781
	782	IMPLICIT NONE
	783
	784	! Declare input arguments
	785	REAL, INTENT(IN) :: i ! Column
	786	REAL, INTENT(IN) :: j ! Row
	787	TYPE (proj_info), INTENT(IN) :: proj
	788
	789	! Declare output arguments
	790	REAL, INTENT(OUT) :: lat ! -90 -> 90 north
	791	REAL, INTENT(OUT) :: lon ! -180 -> 180 East
	792
	793	! Local variables
	794	REAL :: reflon
	795	REAL :: scale_top
	796	REAL :: xx,yy
	797	REAL :: gi2, r2
	798	REAL :: arccos
	799
	800	! Begin Code
	801
	802	! Compute the reference longitude by rotating 90 degrees to the east
	803	! to find the longitude line parallel to the positive x-axis.
	804	reflon = proj%stdlon + 90.
	805
	806	! Compute numerator term of map scale factor
	807	scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)
	808
	809	! Compute radius to point of interest
	810	xx = i - proj%polei
	811	yy = (j - proj%polej) * proj%hemi
	812	r2 = xx2 + yy2
	813
	814	! Now the magic code
	815	IF (r2 .EQ. 0.) THEN
	816	lat = proj%hemi * 90.
	817	lon = reflon
	818	ELSE
	819	gi2 = (proj%rebydx * scale_top)**2.
	820	lat = deg_per_rad * proj%hemi * ASIN((gi2-r2)/(gi2+r2))
	821	arccos = ACOS(xx/SQRT(r2))
	822	IF (yy .GT. 0) THEN
	823	lon = reflon + deg_per_rad * arccos
	824	ELSE
	825	lon = reflon - deg_per_rad * arccos
	826	ENDIF
	827	ENDIF
	828
	829	! Convert to a -180 -> 180 East convention
	830	IF (lon .GT. 180.) lon = lon - 360.
	831	IF (lon .LT. -180.) lon = lon + 360.
	832
	833	RETURN
	834
	835	END SUBROUTINE ijll_ps
	836
	837
	838	SUBROUTINE set_ps_wgs84(proj)
	839	! Initializes a polar-stereographic map projection (WGS84 ellipsoid)
	840	! from the partially filled proj structure. This routine computes the
	841	! radius to the southwest corner and computes the i/j location of the
	842	! pole for use in llij_ps and ijll_ps.
	843
	844	IMPLICIT NONE
	845
	846	! Arguments
	847	TYPE(proj_info), INTENT(INOUT) :: proj
	848
	849	! Local variables
	850	real :: h, mc, tc, t, rho
	851
	852	h = proj%hemi
	853
	854	mc = cos(hproj%truelat1rad_per_deg)/sqrt(1.0-(E_WGS84sin(hproj%truelat1rad_per_deg))*2.0)
	855	tc = sqrt(((1.0-sin(hproj%truelat1rad_per_deg))/(1.0+sin(hproj%truelat1rad_per_deg)))* &
	856	(((1.0+E_WGS84sin(hproj%truelat1rad_per_deg))/(1.0-E_WGS84sin(hproj%truelat1rad_per_deg)))**E_WGS84 ))
	857
	858	! Find the i/j location of reference lat/lon with respect to the pole of the projection
	859	t = sqrt(((1.0-sin(hproj%lat1rad_per_deg))/(1.0+sin(hproj%lat1rad_per_deg)))* &
	860	(((1.0+E_WGS84sin(hproj%lat1rad_per_deg))/(1.0-E_WGS84sin(hproj%lat1rad_per_deg)) )**E_WGS84 ) )
	861	rho = h * (A_WGS84 / proj%dx) * mc * t / tc
	862	proj%polei = rho * sin((hproj%lon1 - hproj%stdlon)*rad_per_deg)
	863	proj%polej = -rho * cos((hproj%lon1 - hproj%stdlon)*rad_per_deg)
	864
	865	RETURN
	866
	867	END SUBROUTINE set_ps_wgs84
	868
	869
	870	SUBROUTINE llij_ps_wgs84(lat,lon,proj,i,j)
	871	! Given latitude (-90 to 90), longitude (-180 to 180), and the
	872	! standard polar-stereographic projection information via the
	873	! public proj structure, this routine returns the i/j indices which
	874	! if within the domain range from 1->nx and 1->ny, respectively.
	875
	876	IMPLICIT NONE
	877
	878	! Arguments
	879	REAL, INTENT(IN) :: lat
	880	REAL, INTENT(IN) :: lon
	881	REAL, INTENT(OUT) :: i !(x-index)
	882	REAL, INTENT(OUT) :: j !(y-index)
	883	TYPE(proj_info),INTENT(IN) :: proj
	884
	885	! Local variables
	886	real :: h, mc, tc, t, rho
	887
	888	h = proj%hemi
	889
	890	mc = cos(hproj%truelat1rad_per_deg)/sqrt(1.0-(E_WGS84sin(hproj%truelat1rad_per_deg))*2.0)
	891	tc = sqrt(((1.0-sin(hproj%truelat1rad_per_deg))/(1.0+sin(hproj%truelat1rad_per_deg)))* &
	892	(((1.0+E_WGS84sin(hproj%truelat1rad_per_deg))/(1.0-E_WGS84sin(hproj%truelat1rad_per_deg)))**E_WGS84 ))
	893
	894	t = sqrt(((1.0-sin(hlatrad_per_deg))/(1.0+sin(hlatrad_per_deg))) * &
	895	(((1.0+E_WGS84sin(hlatrad_per_deg))/(1.0-E_WGS84sin(hlatrad_per_deg)))**E_WGS84))
	896
	897	! Find the x/y location of the requested lat/lon with respect to the pole of the projection
	898	rho = (A_WGS84 / proj%dx) * mc * t / tc
	899	i = h * rho * sin((hlon - hproj%stdlon)*rad_per_deg)
	900	j = h (-rho) cos((hlon - hproj%stdlon)*rad_per_deg)
	901
	902	! Get i/j relative to reference i/j
	903	i = proj%knowni + (i - proj%polei)
	904	j = proj%knownj + (j - proj%polej)
	905
	906	RETURN
	907
	908	END SUBROUTINE llij_ps_wgs84
	909
	910
	911	SUBROUTINE ijll_ps_wgs84(i, j, proj, lat, lon)
	912
	913	! This is the inverse subroutine of llij_ps. It returns the
	914	! latitude and longitude of an i/j point given the projection info
	915	! structure.
	916
	917	implicit none
	918
	919	! Arguments
	920	REAL, INTENT(IN) :: i ! Column
	921	REAL, INTENT(IN) :: j ! Row
	922	REAL, INTENT(OUT) :: lat ! -90 -> 90 north
	923	REAL, INTENT(OUT) :: lon ! -180 -> 180 East
	924	TYPE (proj_info), INTENT(IN) :: proj
	925
	926	! Local variables
	927	real :: h, mc, tc, t, rho, x, y
	928	real :: chi, a, b, c, d
	929
	930	h = proj%hemi
	931	x = (i - proj%knowni + proj%polei)
	932	y = (j - proj%knownj + proj%polej)
	933
	934	mc = cos(hproj%truelat1rad_per_deg)/sqrt(1.0-(E_WGS84sin(hproj%truelat1rad_per_deg))*2.0)
	935	tc = sqrt(((1.0-sin(hproj%truelat1rad_per_deg))/(1.0+sin(hproj%truelat1rad_per_deg))) * &
	936	(((1.0+E_WGS84sin(hproj%truelat1rad_per_deg))/(1.0-E_WGS84sin(hproj%truelat1rad_per_deg)))**E_WGS84 ))
	937
	938	rho = sqrt((xproj%dx)2.0 + (yproj%dx)**2.0)
	939	t = rho * tc / (A_WGS84 * mc)
	940
	941	lon = hproj%stdlon + hatan2(hx,h(-y))
	942
	943	chi = PI/2.0-2.0*atan(t)
	944	a = 1./2.E_WGS842. + 5./24.E_WGS84*4. + 1./40.E_WGS84*6. + 73./2016.E_WGS84**8.
	945	b = 7./24.E_WGS844. + 29./120.E_WGS84*6. + 54113./40320.E_WGS84**8.
	946	c = 7./30.E_WGS846. + 81./280.E_WGS84**8.
	947	d = 4279./20160.E_WGS84*8.
	948
	949	lat = chi + sin(2.chi)(a + cos(2.chi)(b + cos(2.chi)(c + dcos(2.chi))))
	950	lat = h * lat
	951
	952	lat = lat*deg_per_rad
	953	lon = lon*deg_per_rad
	954
	955	RETURN
	956
	957	END SUBROUTINE ijll_ps_wgs84
	958
	959
	960	SUBROUTINE set_albers_nad83(proj)
	961	! Initializes an Albers equal area map projection (NAD83 ellipsoid)
	962	! from the partially filled proj structure. This routine computes the
	963	! radius to the southwest corner and computes the i/j location of the
	964	! pole for use in llij_albers_nad83 and ijll_albers_nad83.
	965
	966	IMPLICIT NONE
	967
	968	! Arguments
	969	TYPE(proj_info), INTENT(INOUT) :: proj
	970
	971	! Local variables
	972	real :: h, m1, m2, q1, q2, theta, q, sinphi
	973
	974	h = proj%hemi
	975
	976	m1 = cos(hproj%truelat1rad_per_deg)/sqrt(1.0-(E_NAD83sin(hproj%truelat1rad_per_deg))*2.0)
	977	m2 = cos(hproj%truelat2rad_per_deg)/sqrt(1.0-(E_NAD83sin(hproj%truelat2rad_per_deg))*2.0)
	978
	979	sinphi = sin(proj%truelat1*rad_per_deg)
	980	q1 = (1.0-E_NAD83*2.0) &
	981	((sinphi/(1.0-(E_NAD83sinphi)2.0)) - 1.0/(2.0E_NAD83) * log((1.0-E_NAD83sinphi)/(1.0+E_NAD83sinphi)))
	982
	983	sinphi = sin(proj%truelat2*rad_per_deg)
	984	q2 = (1.0-E_NAD83*2.0) &
	985	((sinphi/(1.0-(E_NAD83sinphi)2.0)) - 1.0/(2.0E_NAD83) * log((1.0-E_NAD83sinphi)/(1.0+E_NAD83sinphi)))
	986
	987	if (proj%truelat1 == proj%truelat2) then
	988	proj%nc = sin(proj%truelat1*rad_per_deg)
	989	else
	990	proj%nc = (m12.0 - m22.0) / (q2 - q1)
	991	end if
	992
	993	proj%bigc = m1*2.0 + proj%ncq1
	994
	995	! Find the i/j location of reference lat/lon with respect to the pole of the projection
	996	sinphi = sin(proj%lat1*rad_per_deg)
	997	q = (1.0-E_NAD83*2.0) &
	998	((sinphi/(1.0-(E_NAD83sinphi)2.0)) - 1.0/(2.0E_NAD83) * log((1.0-E_NAD83sinphi)/(1.0+E_NAD83sinphi)))
	999
	1000	proj%rho0 = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc
	1001	theta = proj%nc(proj%lon1 - proj%stdlon)rad_per_deg
	1002
	1003	proj%polei = proj%rho0 * sin(h*theta)
	1004	proj%polej = proj%rho0 - proj%rho0 * cos(h*theta)
	1005
	1006	RETURN
	1007
	1008	END SUBROUTINE set_albers_nad83
	1009
	1010
	1011	SUBROUTINE llij_albers_nad83(lat,lon,proj,i,j)
	1012	! Given latitude (-90 to 90), longitude (-180 to 180), and the
	1013	! standard projection information via the
	1014	! public proj structure, this routine returns the i/j indices which
	1015	! if within the domain range from 1->nx and 1->ny, respectively.
	1016
	1017	IMPLICIT NONE
	1018
	1019	! Arguments
	1020	REAL, INTENT(IN) :: lat
	1021	REAL, INTENT(IN) :: lon
	1022	REAL, INTENT(OUT) :: i !(x-index)
	1023	REAL, INTENT(OUT) :: j !(y-index)
	1024	TYPE(proj_info),INTENT(IN) :: proj
	1025
	1026	! Local variables
	1027	real :: h, q, rho, theta, sinphi
	1028
	1029	h = proj%hemi
	1030
	1031	sinphi = sin(hlatrad_per_deg)
	1032
	1033	! Find the x/y location of the requested lat/lon with respect to the pole of the projection
	1034	q = (1.0-E_NAD83*2.0) &
	1035	((sinphi/(1.0-(E_NAD83sinphi)2.0)) - 1.0/(2.0E_NAD83) * log((1.0-E_NAD83sinphi)/(1.0+E_NAD83sinphi)))
	1036
	1037	rho = h * (A_NAD83 / proj%dx) * sqrt(proj%bigc - proj%nc * q) / proj%nc
	1038	theta = proj%nc * (hlon - hproj%stdlon)*rad_per_deg
	1039
	1040	i = hrhosin(theta)
	1041	j = hproj%rho0 - hrho*cos(theta)
	1042
	1043	! Get i/j relative to reference i/j
	1044	i = proj%knowni + (i - proj%polei)
	1045	j = proj%knownj + (j - proj%polej)
	1046
	1047	RETURN
	1048
	1049	END SUBROUTINE llij_albers_nad83
	1050
	1051
	1052	SUBROUTINE ijll_albers_nad83(i, j, proj, lat, lon)
	1053
	1054	! This is the inverse subroutine of llij_albers_nad83. It returns the
	1055	! latitude and longitude of an i/j point given the projection info
	1056	! structure.
	1057
	1058	implicit none
	1059
	1060	! Arguments
	1061	REAL, INTENT(IN) :: i ! Column
	1062	REAL, INTENT(IN) :: j ! Row
	1063	REAL, INTENT(OUT) :: lat ! -90 -> 90 north
	1064	REAL, INTENT(OUT) :: lon ! -180 -> 180 East
	1065	TYPE (proj_info), INTENT(IN) :: proj
	1066
	1067	! Local variables
	1068	real :: h, q, rho, theta, beta, x, y
	1069	real :: a, b, c
	1070
	1071	h = proj%hemi
	1072
	1073	x = (i - proj%knowni + proj%polei)
	1074	y = (j - proj%knownj + proj%polej)
	1075
	1076	rho = sqrt(x2.0 + (proj%rho0 - y)2.0)
	1077	theta = atan2(x, proj%rho0-y)
	1078
	1079	q = (proj%bigc - (rhoproj%ncproj%dx/A_NAD83)**2.0) / proj%nc
	1080
	1081	beta = asin(q/(1.0 - log((1.0-E_NAD83)/(1.0+E_NAD83))(1.0-E_NAD832.0)/(2.0E_NAD83)))
	1082	a = 1./3.E_NAD832. + 31./180.E_NAD83*4. + 517./5040.E_NAD83**6.
	1083	b = 23./360.E_NAD834. + 251./3780.E_NAD83**6.
	1084	c = 761./45360.E_NAD83*6.
	1085
	1086	lat = beta + asin(2.beta) + bsin(4.beta) + csin(6.beta)
	1087
	1088	lat = hlatdeg_per_rad
	1089	lon = proj%stdlon + theta*deg_per_rad/proj%nc
	1090
	1091	RETURN
	1092
	1093	END SUBROUTINE ijll_albers_nad83
	1094
	1095
	1096	SUBROUTINE set_lc(proj)
	1097	! Initialize the remaining items in the proj structure for a
	1098	! lambert conformal grid.
	1099
	1100	IMPLICIT NONE
	1101
	1102	TYPE(proj_info), INTENT(INOUT) :: proj
	1103
	1104	REAL :: arg
	1105	REAL :: deltalon1
	1106	REAL :: tl1r
	1107	REAL :: ctl1r
	1108
	1109	! Compute cone factor
	1110	CALL lc_cone(proj%truelat1, proj%truelat2, proj%cone)
	1111
	1112	! Compute longitude differences and ensure we stay out of the
	1113	! forbidden "cut zone"
	1114	deltalon1 = proj%lon1 - proj%stdlon
	1115	IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
	1116	IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.
	1117
	1118	! Convert truelat1 to radian and compute COS for later use
	1119	tl1r = proj%truelat1 * rad_per_deg
	1120	ctl1r = COS(tl1r)
	1121
	1122	! Compute the radius to our known lower-left (SW) corner
	1123	proj%rsw = proj%rebydx * ctl1r/proj%cone * &
	1124	(TAN((90.proj%hemi-proj%lat1)rad_per_deg/2.) / &
	1125	TAN((90.proj%hemi-proj%truelat1)rad_per_deg/2.))**proj%cone
	1126
	1127	! Find pole point
	1128	arg = proj%cone(deltalon1rad_per_deg)
	1129	proj%polei = proj%hemiproj%knowni - proj%hemi proj%rsw * SIN(arg)
	1130	proj%polej = proj%hemiproj%knownj + proj%rsw COS(arg)
	1131
	1132	RETURN
	1133
	1134	END SUBROUTINE set_lc
	1135
	1136
	1137	SUBROUTINE lc_cone(truelat1, truelat2, cone)
	1138
	1139	! Subroutine to compute the cone factor of a Lambert Conformal projection
	1140
	1141	IMPLICIT NONE
	1142
	1143	! Input Args
	1144	REAL, INTENT(IN) :: truelat1 ! (-90 -> 90 degrees N)
	1145	REAL, INTENT(IN) :: truelat2 ! " " " " "
	1146
	1147	! Output Args
	1148	REAL, INTENT(OUT) :: cone
	1149
	1150	! Locals
	1151
	1152	! BEGIN CODE
	1153
	1154	! First, see if this is a secant or tangent projection. For tangent
	1155	! projections, truelat1 = truelat2 and the cone is tangent to the
	1156	! Earth's surface at this latitude. For secant projections, the cone
	1157	! intersects the Earth's surface at each of the distinctly different
	1158	! latitudes
	1159	IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
	1160	cone = ALOG10(COS(truelat1*rad_per_deg)) - &
	1161	ALOG10(COS(truelat2*rad_per_deg))
	1162	cone = cone /(ALOG10(TAN((45.0 - ABS(truelat1)/2.0) * rad_per_deg)) - &
	1163	ALOG10(TAN((45.0 - ABS(truelat2)/2.0) * rad_per_deg)))
	1164	ELSE
	1165	cone = SIN(ABS(truelat1)*rad_per_deg )
	1166	ENDIF
	1167
	1168	RETURN
	1169
	1170	END SUBROUTINE lc_cone
	1171
	1172
	1173	SUBROUTINE ijll_lc( i, j, proj, lat, lon)
	1174
	1175	! Subroutine to convert from the (i,j) cartesian coordinate to the
	1176	! geographical latitude and longitude for a Lambert Conformal projection.
	1177
	1178	! History:
	1179	! 25 Jul 01: Corrected by B. Shaw, NOAA/FSL
	1180	!
	1181	IMPLICIT NONE
	1182
	1183	! Input Args
	1184	REAL, INTENT(IN) :: i ! Cartesian X coordinate
	1185	REAL, INTENT(IN) :: j ! Cartesian Y coordinate
	1186	TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
	1187
	1188	! Output Args
	1189	REAL, INTENT(OUT) :: lat ! Latitude (-90->90 deg N)
	1190	REAL, INTENT(OUT) :: lon ! Longitude (-180->180 E)
	1191
	1192	! Locals
	1193	REAL :: inew
	1194	REAL :: jnew
	1195	REAL :: r
	1196	REAL :: chi,chi1,chi2
	1197	REAL :: r2
	1198	REAL :: xx
	1199	REAL :: yy
	1200
	1201	! BEGIN CODE
	1202
	1203	chi1 = (90. - proj%hemiproj%truelat1)rad_per_deg
	1204	chi2 = (90. - proj%hemiproj%truelat2)rad_per_deg
	1205
	1206	! See if we are in the southern hemispere and flip the indices
	1207	! if we are.
	1208	inew = proj%hemi * i
	1209	jnew = proj%hemi * j
	1210
	1211	! Compute radius**2 to i/j location
	1212	xx = inew - proj%polei
	1213	yy = proj%polej - jnew
	1214	r2 = (xxxx + yyyy)
	1215	r = SQRT(r2)/proj%rebydx
	1216
	1217	! Convert to lat/lon
	1218	IF (r2 .EQ. 0.) THEN
	1219	lat = proj%hemi * 90.
	1220	lon = proj%stdlon
	1221	ELSE
	1222
	1223	! Longitude
	1224	lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
	1225	lon = AMOD(lon+360., 360.)
	1226
	1227	! Latitude. Latitude determined by solving an equation adapted
	1228	! from:
	1229	! Maling, D.H., 1973: Coordinate Systems and Map Projections
	1230	! Equations #20 in Appendix I.
	1231
	1232	IF (chi1 .EQ. chi2) THEN
	1233	chi = 2.0ATAN( ( r/TAN(chi1) )(1./proj%cone) TAN(chi1*0.5) )
	1234	ELSE
	1235	chi = 2.0ATAN( (rproj%cone/SIN(chi1))*(1./proj%cone) TAN(chi1*0.5))
	1236	ENDIF
	1237	lat = (90.0-chideg_per_rad)proj%hemi
	1238
	1239	ENDIF
	1240
	1241	IF (lon .GT. +180.) lon = lon - 360.
	1242	IF (lon .LT. -180.) lon = lon + 360.
	1243
	1244	RETURN
	1245
	1246	END SUBROUTINE ijll_lc
	1247
	1248
	1249	SUBROUTINE llij_lc( lat, lon, proj, i, j)
	1250
	1251	! Subroutine to compute the geographical latitude and longitude values
	1252	! to the cartesian x/y on a Lambert Conformal projection.
	1253
	1254	IMPLICIT NONE
	1255
	1256	! Input Args
	1257	REAL, INTENT(IN) :: lat ! Latitude (-90->90 deg N)
	1258	REAL, INTENT(IN) :: lon ! Longitude (-180->180 E)
	1259	TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure
	1260
	1261	! Output Args
	1262	REAL, INTENT(OUT) :: i ! Cartesian X coordinate
	1263	REAL, INTENT(OUT) :: j ! Cartesian Y coordinate
	1264
	1265	! Locals
	1266	REAL :: arg
	1267	REAL :: deltalon
	1268	REAL :: tl1r
	1269	REAL :: rm
	1270	REAL :: ctl1r
	1271
	1272
	1273	! BEGIN CODE
	1274
	1275	! Compute deltalon between known longitude and standard lon and ensure
	1276	! it is not in the cut zone
	1277	deltalon = lon - proj%stdlon
	1278	IF (deltalon .GT. +180.) deltalon = deltalon - 360.
	1279	IF (deltalon .LT. -180.) deltalon = deltalon + 360.
	1280
	1281	! Convert truelat1 to radian and compute COS for later use
	1282	tl1r = proj%truelat1 * rad_per_deg
	1283	ctl1r = COS(tl1r)
	1284
	1285	! Radius to desired point
	1286	rm = proj%rebydx * ctl1r/proj%cone * &
	1287	(TAN((90.proj%hemi-lat)rad_per_deg/2.) / &
	1288	TAN((90.proj%hemi-proj%truelat1)rad_per_deg/2.))**proj%cone
	1289
	1290	arg = proj%cone(deltalonrad_per_deg)
	1291	i = proj%polei + proj%hemi * rm * SIN(arg)
	1292	j = proj%polej - rm * COS(arg)
	1293
	1294	! Finally, if we are in the southern hemisphere, flip the i/j
	1295	! values to a coordinate system where (1,1) is the SW corner
	1296	! (what we assume) which is different than the original NCEP
	1297	! algorithms which used the NE corner as the origin in the
	1298	! southern hemisphere (left-hand vs. right-hand coordinate?)
	1299	i = proj%hemi * i
	1300	j = proj%hemi * j
	1301
	1302	RETURN
	1303	END SUBROUTINE llij_lc
	1304
	1305
	1306	SUBROUTINE set_merc(proj)
	1307
	1308	! Sets up the remaining basic elements for the mercator projection
	1309
	1310	IMPLICIT NONE
	1311	TYPE(proj_info), INTENT(INOUT) :: proj
	1312	REAL :: clain
	1313
	1314
	1315	! Preliminary variables
	1316
	1317	clain = COS(rad_per_deg*proj%truelat1)
	1318	proj%dlon = proj%dx / (proj%re_m * clain)
	1319
	1320	! Compute distance from equator to origin, and store in the
	1321	! proj%rsw tag.
	1322
	1323	proj%rsw = 0.
	1324	IF (proj%lat1 .NE. 0.) THEN
	1325	proj%rsw = (ALOG(TAN(0.5((proj%lat1+90.)rad_per_deg))))/proj%dlon
	1326	ENDIF
	1327
	1328	RETURN
	1329
	1330	END SUBROUTINE set_merc
	1331
	1332
	1333	SUBROUTINE llij_merc(lat, lon, proj, i, j)
	1334
	1335	! Compute i/j coordinate from lat lon for mercator projection
	1336
	1337	IMPLICIT NONE
	1338	REAL, INTENT(IN) :: lat
	1339	REAL, INTENT(IN) :: lon
	1340	TYPE(proj_info),INTENT(IN) :: proj
	1341	REAL,INTENT(OUT) :: i
	1342	REAL,INTENT(OUT) :: j
	1343	REAL :: deltalon
	1344
	1345	deltalon = lon - proj%lon1
	1346	IF (deltalon .LT. -180.) deltalon = deltalon + 360.
	1347	IF (deltalon .GT. 180.) deltalon = deltalon - 360.
	1348	i = proj%knowni + (deltalon/(proj%dlon*deg_per_rad))
	1349	j = proj%knownj + (ALOG(TAN(0.5((lat + 90.) rad_per_deg)))) / &
	1350	proj%dlon - proj%rsw
	1351
	1352	RETURN
	1353
	1354	END SUBROUTINE llij_merc
	1355
	1356
	1357	SUBROUTINE ijll_merc(i, j, proj, lat, lon)
	1358
	1359	! Compute the lat/lon from i/j for mercator projection
	1360
	1361	IMPLICIT NONE
	1362	REAL,INTENT(IN) :: i
	1363	REAL,INTENT(IN) :: j
	1364	TYPE(proj_info),INTENT(IN) :: proj
	1365	REAL, INTENT(OUT) :: lat
	1366	REAL, INTENT(OUT) :: lon
	1367
	1368
	1369	lat = 2.0ATAN(EXP(proj%dlon(proj%rsw + j-proj%knownj)))*deg_per_rad - 90.
	1370	lon = (i-proj%knowni)proj%dlondeg_per_rad + proj%lon1
	1371	IF (lon.GT.180.) lon = lon - 360.
	1372	IF (lon.LT.-180.) lon = lon + 360.
	1373	RETURN
	1374
	1375	END SUBROUTINE ijll_merc
	1376
	1377
	1378	SUBROUTINE llij_latlon(lat, lon, proj, i, j)
	1379
	1380	! Compute the i/j location of a lat/lon on a LATLON grid.
	1381	IMPLICIT NONE
	1382	REAL, INTENT(IN) :: lat
	1383	REAL, INTENT(IN) :: lon
	1384	TYPE(proj_info), INTENT(IN) :: proj
	1385	REAL, INTENT(OUT) :: i
	1386	REAL, INTENT(OUT) :: j
	1387
	1388	REAL :: deltalat
	1389	REAL :: deltalon
	1390
	1391	! Compute deltalat and deltalon as the difference between the input
	1392	! lat/lon and the origin lat/lon
	1393	deltalat = lat - proj%lat1
	1394	deltalon = lon - proj%lon1
	1395
	1396	! Compute i/j
	1397	i = deltalon/proj%loninc
	1398	j = deltalat/proj%latinc
	1399
	1400	i = i + proj%knowni
	1401	j = j + proj%knownj
	1402
	1403	RETURN
	1404
	1405	END SUBROUTINE llij_latlon
	1406
	1407
	1408	SUBROUTINE ijll_latlon(i, j, proj, lat, lon)
	1409
	1410	! Compute the lat/lon location of an i/j on a LATLON grid.
	1411	IMPLICIT NONE
	1412	REAL, INTENT(IN) :: i
	1413	REAL, INTENT(IN) :: j
	1414	TYPE(proj_info), INTENT(IN) :: proj
	1415	REAL, INTENT(OUT) :: lat
	1416	REAL, INTENT(OUT) :: lon
	1417
	1418	REAL :: i_work, j_work
	1419	REAL :: deltalat
	1420	REAL :: deltalon
	1421
	1422	i_work = i - proj%knowni
	1423	j_work = j - proj%knownj
	1424
	1425	! Compute deltalat and deltalon
	1426	deltalat = j_work*proj%latinc
	1427	deltalon = i_work*proj%loninc
	1428
	1429	lat = proj%lat1 + deltalat
	1430	lon = proj%lon1 + deltalon
	1431
	1432	RETURN
	1433
	1434	END SUBROUTINE ijll_latlon
	1435
	1436
	1437	SUBROUTINE set_cyl(proj)
	1438
	1439	implicit none
	1440
	1441	! Arguments
	1442	type(proj_info), intent(inout) :: proj
	1443
	1444	proj%hemi = 1.0
	1445
	1446	END SUBROUTINE set_cyl
	1447
	1448
	1449	SUBROUTINE llij_cyl(lat, lon, proj, i, j)
	1450
	1451	implicit none
	1452
	1453	! Arguments
	1454	real, intent(in) :: lat, lon
	1455	real, intent(out) :: i, j
	1456	type(proj_info), intent(in) :: proj
	1457
	1458	! Local variables
	1459	real :: deltalat
	1460	real :: deltalon
	1461
	1462	! Compute deltalat and deltalon as the difference between the input
	1463	! lat/lon and the origin lat/lon
	1464	deltalat = lat - proj%lat1
	1465	! deltalon = lon - proj%stdlon
	1466	deltalon = lon - proj%lon1
	1467
	1468	if (deltalon < 0.) deltalon = deltalon + 360.
	1469	if (deltalon > 360.) deltalon = deltalon - 360.
	1470
	1471	! Compute i/j
	1472	i = deltalon/proj%loninc
	1473	j = deltalat/proj%latinc
	1474
	1475	if (i <= 0.) i = i + 360./proj%loninc
	1476	if (i > 360./proj%loninc) i = i - 360./proj%loninc
	1477
	1478	i = i + proj%knowni
	1479	j = j + proj%knownj
	1480
	1481	END SUBROUTINE llij_cyl
	1482
	1483
	1484	SUBROUTINE ijll_cyl(i, j, proj, lat, lon)
	1485
	1486	implicit none
	1487
	1488	! Arguments
	1489	real, intent(in) :: i, j
	1490	real, intent(out) :: lat, lon
	1491	type(proj_info), intent(in) :: proj
	1492
	1493	! Local variables
	1494	real :: deltalat
	1495	real :: deltalon
	1496	real :: i_work, j_work
	1497
	1498	i_work = i - proj%knowni
	1499	j_work = j - proj%knownj
	1500
	1501	if (i_work < 0.) i_work = i_work + 360./proj%loninc
	1502	if (i_work >= 360./proj%loninc) i_work = i_work - 360./proj%loninc
	1503
	1504	! Compute deltalat and deltalon
	1505	deltalat = j_work*proj%latinc
	1506	deltalon = i_work*proj%loninc
	1507
	1508	lat = deltalat + proj%lat1
	1509	! lon = deltalon + proj%stdlon
	1510	lon = deltalon + proj%lon1
	1511
	1512	if (lon < -180.) lon = lon + 360.
	1513	if (lon > 180.) lon = lon - 360.
	1514
	1515	END SUBROUTINE ijll_cyl
	1516
	1517
	1518	SUBROUTINE set_cassini(proj)
	1519
	1520	implicit none
	1521
	1522	! Arguments
	1523	type(proj_info), intent(inout) :: proj
	1524
	1525	! Local variables
	1526	real :: comp_lat, comp_lon
	1527	logical :: global_domain
	1528
	1529	proj%hemi = 1.0
	1530
	1531	! Try to determine whether this domain has global coverage
	1532	if (abs(proj%lat1 - proj%latinc/2. + 90.) < 0.001 .and. &
	1533	abs(mod(proj%lon1 - proj%loninc/2. - proj%stdlon,360.)) < 0.001) then
	1534	global_domain = .true.
	1535	else
	1536	global_domain = .false.
	1537	end if
	1538
	1539	if (abs(proj%lat0) /= 90. .and. .not.global_domain) then
	1540	call rotate_coords(proj%lat1,proj%lon1,comp_lat,comp_lon,proj%lat0,proj%lon0,proj%stdlon,-1)
	1541	proj%lat1 = comp_lat
	1542	proj%lon1 = comp_lon
	1543	end if
	1544
	1545	END SUBROUTINE set_cassini
	1546
	1547
	1548	SUBROUTINE llij_cassini(lat, lon, proj, i, j)
	1549
	1550	implicit none
	1551
	1552	! Arguments
	1553	real, intent(in) :: lat, lon
	1554	real, intent(out) :: i, j
	1555	type(proj_info), intent(in) :: proj
	1556
	1557	! Local variables
	1558	real :: comp_lat, comp_lon
	1559
	1560	! Convert geographic to computational lat/lon
	1561	if (abs(proj%lat0) /= 90.) then
	1562	call rotate_coords(lat,lon,comp_lat,comp_lon,proj%lat0,proj%lon0,proj%stdlon,-1)
	1563	else
	1564	comp_lat = lat
	1565	comp_lon = lon
	1566	end if
	1567
	1568	! Convert computational lat/lon to i/j
	1569	call llij_cyl(comp_lat, comp_lon, proj, i, j)
	1570
	1571	END SUBROUTINE llij_cassini
	1572
	1573
	1574	SUBROUTINE ijll_cassini(i, j, proj, lat, lon)
	1575
	1576	implicit none
	1577
	1578	! Arguments
	1579	real, intent(in) :: i, j
	1580	real, intent(out) :: lat, lon
	1581	type(proj_info), intent(in) :: proj
	1582
	1583	! Local variables
	1584	real :: comp_lat, comp_lon
	1585
	1586	! Convert i/j to computational lat/lon
	1587	call ijll_cyl(i, j, proj, comp_lat, comp_lon)
	1588
	1589	! Convert computational to geographic lat/lon
	1590	if (abs(proj%lat0) /= 90.) then
	1591	call rotate_coords(comp_lat,comp_lon,lat,lon,proj%lat0,proj%lon0,proj%stdlon,1)
	1592	else
	1593	lat = comp_lat
	1594	lon = comp_lon
	1595	end if
	1596
	1597	END SUBROUTINE ijll_cassini
	1598
	1599
	1600	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1601	! Purpose: Converts between computational and geographic lat/lon for Cassini
	1602	!
	1603	! Notes: This routine was provided by Bill Skamarock, 2007-03-27
	1604	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1605	SUBROUTINE rotate_coords(ilat,ilon,olat,olon,lat_np,lon_np,lon_0,direction)
	1606
	1607	IMPLICIT NONE
	1608
	1609	REAL, INTENT(IN ) :: ilat, ilon
	1610	REAL, INTENT( OUT) :: olat, olon
	1611	REAL, INTENT(IN ) :: lat_np, lon_np, lon_0
	1612	INTEGER, INTENT(IN ), OPTIONAL :: direction
	1613	! >=0, default : computational -> geographical
	1614	! < 0 : geographical -> computational
	1615
	1616	REAL :: rlat, rlon
	1617	REAL :: phi_np, lam_np, lam_0, dlam
	1618	REAL :: sinphi, cosphi, coslam, sinlam
	1619
	1620	! Convert all angles to radians
	1621	phi_np = lat_np * rad_per_deg
	1622	lam_np = lon_np * rad_per_deg
	1623	lam_0 = lon_0 * rad_per_deg
	1624	rlat = ilat * rad_per_deg
	1625	rlon = ilon * rad_per_deg
	1626
	1627	IF (PRESENT(direction) .AND. (direction < 0)) THEN
	1628	! The equations are exactly the same except for one small difference
	1629	! with respect to longitude ...
	1630	dlam = PI - lam_0
	1631	ELSE
	1632	dlam = lam_np
	1633	END IF
	1634	sinphi = COS(phi_np)COS(rlat)COS(rlon-dlam) + SIN(phi_np)*SIN(rlat)
	1635	cosphi = SQRT(1.-sinphi*sinphi)
	1636	coslam = SIN(phi_np)COS(rlat)COS(rlon-dlam) - COS(phi_np)*SIN(rlat)
	1637	sinlam = COS(rlat)*SIN(rlon-dlam)
	1638	IF ( cosphi /= 0. ) THEN
	1639	coslam = coslam/cosphi
	1640	sinlam = sinlam/cosphi
	1641	END IF
	1642	olat = deg_per_rad*ASIN(sinphi)
	1643	olon = deg_per_rad*(ATAN2(sinlam,coslam)-dlam-lam_0+lam_np)
	1644	! Both of my F90 text books prefer the DO-EXIT form, and claim it is faster
	1645	! when optimization is turned on (as we will always do...)
	1646	DO
	1647	IF (olon >= -180.) EXIT
	1648	olon = olon + 360.
	1649	END DO
	1650	DO
	1651	IF (olon <= 180.) EXIT
	1652	olon = olon - 360.
	1653	END DO
	1654
	1655	END SUBROUTINE rotate_coords
	1656
	1657
	1658	SUBROUTINE llij_rotlatlon(lat, lon, proj, i_real, j_real)
	1659
	1660	IMPLICIT NONE
	1661
	1662	! Arguments
	1663	REAL, INTENT(IN) :: lat, lon
	1664	REAL :: i, j
	1665	REAL, INTENT(OUT) :: i_real, j_real
	1666	TYPE (proj_info), INTENT(IN) :: proj
	1667
	1668	! Local variables
	1669	INTEGER :: ii,imt,jj,jmt,k,krows,ncol,nrow,iri
	1670	REAL(KIND=HIGH) :: dphd,dlmd !Grid increments, degrees
	1671	REAL(KIND=HIGH) :: glatd !Geographic latitude, positive north
	1672	REAL(KIND=HIGH) :: glond !Geographic longitude, positive west
	1673	REAL(KIND=HIGH) :: col,d1,d2,d2r,dlm,dlm1,dlm2,dph,glat,glon, &
	1674	pi,r2d,row,tlat,tlat1,tlat2, &
	1675	tlon,tlon1,tlon2,tph0,tlm0,x,y,z
	1676
	1677	glatd = lat
	1678	glond = -lon
	1679
	1680	dphd = proj%phi/REAL((proj%jydim-1)/2)
	1681	dlmd = proj%lambda/REAL(proj%ixdim-1)
	1682
	1683	pi = ACOS(-1.0)
	1684	d2r = pi/180.
	1685	r2d = 1./d2r
	1686
	1687	imt = 2*proj%ixdim-1
	1688	jmt = proj%jydim/2+1
	1689
	1690	glat = glatd*d2r
	1691	glon = glond*d2r
	1692	dph = dphd*d2r
	1693	dlm = dlmd*d2r
	1694	tph0 = proj%lat1*d2r
	1695	tlm0 = -proj%lon1*d2r
	1696
	1697	x = COS(tph0)COS(glat)COS(glon-tlm0)+SIN(tph0)*SIN(glat)
	1698	y = -COS(glat)*SIN(glon-tlm0)
	1699	z = COS(tph0)SIN(glat)-SIN(tph0)COS(glat)*COS(glon-tlm0)
	1700	tlat = r2dATAN(z/SQRT(xx+y*y))
	1701	tlon = r2d*ATAN(y/x)
	1702
	1703	row = tlat/dphd+jmt
	1704	col = tlon/dlmd+proj%ixdim
	1705
	1706	if ( (row - INT(row)) .gt. 0.999) then
	1707	row = row + 0.0002
	1708	else if ( (col - INT(col)) .gt. 0.999) then
	1709	col = col + 0.0002
	1710	end if
	1711
	1712	nrow = INT(row)
	1713	ncol = INT(col)
	1714
	1715	! nrow = NINT(row)
	1716	! ncol = NINT(col)
	1717
	1718	tlat = tlat*d2r
	1719	tlon = tlon*d2r
	1720
	1721
	1722	IF (proj%stagger == HH) THEN
	1723
	1724	IF (mod(nrow,2) .eq. 0) then
	1725	i_real = col / 2.0
	1726	ELSE
	1727	i_real = col / 2.0 + 0.5
	1728	ENDIF
	1729	j_real=row
	1730
	1731
	1732	IF ((abs(MOD(nrow,2)) == 1 .AND. abs(MOD(ncol,2)) == 1) .OR. &
	1733	(MOD(nrow,2) == 0 .AND. MOD(ncol,2) == 0)) THEN
	1734
	1735	tlat1 = (nrow-jmt)*dph
	1736	tlat2 = tlat1+dph
	1737	tlon1 = (ncol-proj%ixdim)*dlm
	1738	tlon2 = tlon1+dlm
	1739
	1740	dlm1 = tlon-tlon1
	1741	dlm2 = tlon-tlon2
	1742	d1 = ACOS(COS(tlat)COS(tlat1)COS(dlm1)+SIN(tlat)*SIN(tlat1))
	1743	d2 = ACOS(COS(tlat)COS(tlat2)COS(dlm2)+SIN(tlat)*SIN(tlat2))
	1744
	1745	IF (d1 > d2) THEN
	1746	nrow = nrow+1
	1747	ncol = ncol+1
	1748	END IF
	1749
	1750	ELSE
	1751	tlat1 = (nrow+1-jmt)*dph
	1752	tlat2 = tlat1-dph
	1753	tlon1 = (ncol-proj%ixdim)*dlm
	1754	tlon2 = tlon1+dlm
	1755	dlm1 = tlon-tlon1
	1756	dlm2 = tlon-tlon2
	1757	d1 = ACOS(COS(tlat)COS(tlat1)COS(dlm1)+SIN(tlat)*SIN(tlat1))
	1758	d2 = ACOS(COS(tlat)COS(tlat2)COS(dlm2)+SIN(tlat)*SIN(tlat2))
	1759
	1760	IF (d1 < d2) THEN
	1761	nrow = nrow+1
	1762	ELSE
	1763	ncol = ncol+1
	1764	END IF
	1765	END IF
	1766
	1767	ELSE IF (proj%stagger == VV) THEN
	1768
	1769	IF (mod(nrow,2) .eq. 0) then
	1770	i_real = col / 2.0 + 0.5
	1771	ELSE
	1772	i_real = col / 2.0
	1773	ENDIF
	1774	j_real=row
	1775
	1776	IF ((MOD(nrow,2) == 0 .AND. abs(MOD(ncol,2)) == 1) .OR. &
	1777	(abs(MOD(nrow,2)) == 1 .AND. MOD(ncol,2) == 0)) THEN
	1778	tlat1 = (nrow-jmt)*dph
	1779	tlat2 = tlat1+dph
	1780	tlon1 = (ncol-proj%ixdim)*dlm
	1781	tlon2 = tlon1+dlm
	1782	dlm1 = tlon-tlon1
	1783	dlm2 = tlon-tlon2
	1784	d1 = ACOS(COS(tlat)COS(tlat1)COS(dlm1)+SIN(tlat)*SIN(tlat1))
	1785	d2 = ACOS(COS(tlat)COS(tlat2)COS(dlm2)+SIN(tlat)*SIN(tlat2))
	1786
	1787	IF (d1 > d2) THEN
	1788	nrow = nrow+1
	1789	ncol = ncol+1
	1790	END IF
	1791
	1792	ELSE
	1793	tlat1 = (nrow+1-jmt)*dph
	1794	tlat2 = tlat1-dph
	1795	tlon1 = (ncol-proj%ixdim)*dlm
	1796	tlon2 = tlon1+dlm
	1797	dlm1 = tlon-tlon1
	1798	dlm2 = tlon-tlon2
	1799	d1 = ACOS(COS(tlat)COS(tlat1)COS(dlm1)+SIN(tlat)*SIN(tlat1))
	1800	d2 = ACOS(COS(tlat)COS(tlat2)COS(dlm2)+SIN(tlat)*SIN(tlat2))
	1801
	1802	IF (d1 < d2) THEN
	1803	nrow = nrow+1
	1804	ELSE
	1805	ncol = ncol+1
	1806	END IF
	1807	END IF
	1808	END IF
	1809
	1810
	1811	!!! Added next line as a Kludge - not yet understood why needed
	1812	if (ncol .le. 0) ncol=ncol-1
	1813
	1814	jj = nrow
	1815	ii = ncol/2
	1816
	1817	IF (proj%stagger == HH) THEN
	1818	IF (abs(MOD(jj,2)) == 1) ii = ii+1
	1819	ELSE IF (proj%stagger == VV) THEN
	1820	IF (MOD(jj,2) == 0) ii=ii+1
	1821	END IF
	1822
	1823	i = REAL(ii)
	1824	j = REAL(jj)
	1825
	1826	END SUBROUTINE llij_rotlatlon
	1827
	1828
	1829	SUBROUTINE ijll_rotlatlon(i, j, proj, lat,lon)
	1830
	1831	IMPLICIT NONE
	1832
	1833	! Arguments
	1834	REAL, INTENT(IN) :: i, j
	1835	REAL, INTENT(OUT) :: lat, lon
	1836	TYPE (proj_info), INTENT(IN) :: proj
	1837
	1838	! Local variables
	1839	INTEGER :: ih,jh
	1840	INTEGER :: midcol,midrow,ncol,iadd1,iadd2,imt,jh2,knrow,krem,kv,nrow
	1841	REAL :: i_work, j_work
	1842	REAL :: dphd,dlmd !Grid increments, degrees
	1843	REAL(KIND=HIGH) :: arg1,arg2,d2r,fctr,glatr,glatd,glond,pi, &
	1844	r2d,tlatd,tlond,tlatr,tlonr,tlm0,tph0
	1845	REAL :: col
	1846
	1847	i_work = i
	1848	j_work = j
	1849
	1850	if ( (j - INT(j)) .gt. 0.999) then
	1851	j_work = j + 0.0002
	1852	endif
	1853
	1854	jh = INT(j_work)
	1855
	1856	dphd = proj%phi/REAL((proj%jydim-1)/2)
	1857	dlmd = proj%lambda/REAL(proj%ixdim-1)
	1858
	1859	pi = ACOS(-1.0)
	1860	d2r = pi/180.
	1861	r2d = 1./d2r
	1862	tph0 = proj%lat1*d2r
	1863	tlm0 = -proj%lon1*d2r
	1864
	1865	midrow = (proj%jydim+1)/2
	1866	midcol = proj%ixdim
	1867
	1868	col = 2*i_work-1+abs(MOD(jh+1,2))
	1869	tlatd = (j_work-midrow)*dphd
	1870	tlond = (col-midcol)*dlmd
	1871
	1872	IF (proj%stagger == VV) THEN
	1873	if (mod(jh,2) .eq. 0) then
	1874	tlond = tlond - DLMD
	1875	else
	1876	tlond = tlond + DLMD
	1877	end if
	1878	END IF
	1879
	1880	tlatr = tlatd*d2r
	1881	tlonr = tlond*d2r
	1882	arg1 = SIN(tlatr)COS(tph0)+COS(tlatr)SIN(tph0)*COS(tlonr)
	1883	glatr = ASIN(arg1)
	1884
	1885	glatd = glatr*r2d
	1886
	1887	arg2 = COS(tlatr)COS(tlonr)/(COS(glatr)COS(tph0))-TAN(glatr)*TAN(tph0)
	1888	IF (ABS(arg2) > 1.) arg2 = ABS(arg2)/arg2
	1889	fctr = 1.
	1890	IF (tlond > 0.) fctr = -1.
	1891
	1892	glond = tlm0r2d+fctrACOS(arg2)*r2d
	1893
	1894	lat = glatd
	1895	lon = -glond
	1896
	1897	IF (lon .GT. +180.) lon = lon - 360.
	1898	IF (lon .LT. -180.) lon = lon + 360.
	1899
	1900	END SUBROUTINE ijll_rotlatlon
	1901
	1902
	1903	SUBROUTINE set_gauss(proj)
	1904
	1905	IMPLICIT NONE
	1906
	1907	! Argument
	1908	type (proj_info), intent(inout) :: proj
	1909
	1910	! Initialize the array that will hold the Gaussian latitudes.
	1911
	1912	IF ( ASSOCIATED( proj%gauss_lat ) ) THEN
	1913	DEALLOCATE ( proj%gauss_lat )
	1914	END IF
	1915
	1916	! Get the needed space for our array.
	1917
	1918	ALLOCATE ( proj%gauss_lat(proj%nlat*2) )
	1919
	1920	! Compute the Gaussian latitudes.
	1921
	1922	CALL gausll( proj%nlat*2 , proj%gauss_lat )
	1923
	1924	! Now, these could be upside down from what we want, so let's check.
	1925	! We take advantage of the equatorial symmetry to remove any sort of
	1926	! array re-ordering.
	1927
	1928	IF ( ABS(proj%gauss_lat(1) - proj%lat1) .GT. 0.01 ) THEN
	1929	proj%gauss_lat = -1. * proj%gauss_lat
	1930	END IF
	1931
	1932	! Just a sanity check.
	1933
	1934	IF ( ABS(proj%gauss_lat(1) - proj%lat1) .GT. 0.01 ) THEN
	1935	PRINT '(A)','Oops, something is not right with the Gaussian latitude computation.'
	1936	PRINT '(A,F8.3,A)','The input data gave the starting latitude as ',proj%lat1,'.'
	1937	PRINT '(A,F8.3,A)','This routine computed the starting latitude as +-',ABS(proj%gauss_lat(1)),'.'
	1938	PRINT '(A,F8.3,A)','The difference is larger than 0.01 degrees, which is not expected.'
	1939	CALL wrf_error_fatal ( 'Gaussian_latitude_computation' )
	1940	END IF
	1941
	1942	END SUBROUTINE set_gauss
	1943
	1944
	1945	SUBROUTINE gausll ( nlat , lat_sp )
	1946
	1947	IMPLICIT NONE
	1948
	1949	INTEGER :: nlat , i
	1950	REAL (KIND=HIGH) , PARAMETER :: pi = 3.141592653589793
	1951	REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat , wos2 , lat
	1952	REAL , DIMENSION(nlat) :: lat_sp
	1953
	1954	CALL lggaus(nlat, cosc, gwt, sinc, colat, wos2)
	1955
	1956	DO i = 1, nlat
	1957	lat(i) = ACOS(sinc(i)) * 180._HIGH / pi
	1958	IF (i.gt.nlat/2) lat(i) = -lat(i)
	1959	END DO
	1960
	1961	lat_sp = REAL(lat)
	1962
	1963	END SUBROUTINE gausll
	1964
	1965
	1966	SUBROUTINE lggaus( nlat, cosc, gwt, sinc, colat, wos2 )
	1967
	1968	IMPLICIT NONE
	1969
	1970	! LGGAUS finds the Gaussian latitudes by finding the roots of the
	1971	! ordinary Legendre polynomial of degree NLAT using Newton's
	1972	! iteration method.
	1973
	1974	! On entry:
	1975	integer NLAT ! the number of latitudes (degree of the polynomial)
	1976
	1977	! On exit: for each Gaussian latitude
	1978	! COSC - cos(colatitude) or sin(latitude)
	1979	! GWT - the Gaussian weights
	1980	! SINC - sin(colatitude) or cos(latitude)
	1981	! COLAT - the colatitudes in radians
	1982	! WOS2 - Gaussian weight over sin**2(colatitude)
	1983
	1984	REAL (KIND=HIGH) , DIMENSION(nlat) :: cosc , gwt , sinc , colat , wos2
	1985	REAL (KIND=HIGH) , PARAMETER :: pi = 3.141592653589793
	1986
	1987	! Convergence criterion for iteration of cos latitude
	1988
	1989	REAL , PARAMETER :: xlim = 1.0E-14
	1990
	1991	INTEGER :: nzero, i, j
	1992	REAL (KIND=HIGH) :: fi, fi1, a, b, g, gm, gp, gt, delta, c, d
	1993
	1994	! The number of zeros between pole and equator
	1995
	1996	nzero = nlat/2
	1997
	1998	! Set first guess for cos(colat)
	1999
	2000	DO i=1,nzero
	2001	cosc(i) = SIN( (i-0.5)pi/nlat + pi0.5 )
	2002	END DO
	2003
	2004	! Constants for determining the derivative of the polynomial
	2005	fi = nlat
	2006	fi1 = fi+1.0
	2007	a = fifi1 / SQRT(4.0fi1*fi1-1.0)
	2008	b = fi1fi / SQRT(4.0fi*fi-1.0)
	2009
	2010	! Loop over latitudes, iterating the search for each root
	2011
	2012	DO i=1,nzero
	2013	j=0
	2014
	2015	! Determine the value of the ordinary Legendre polynomial for
	2016	! the current guess root
	2017
	2018	DO
	2019	CALL lgord( g, cosc(i), nlat )
	2020
	2021	! Determine the derivative of the polynomial at this point
	2022
	2023	CALL lgord( gm, cosc(i), nlat-1 )
	2024	CALL lgord( gp, cosc(i), nlat+1 )
	2025	gt = (cosc(i)cosc(i)-1.0) / (agp-b*gm)
	2026
	2027	! Update the estimate of the root
	2028
	2029	delta = g*gt
	2030	cosc(i) = cosc(i) - delta
	2031
	2032	! If convergence criterion has not been met, keep trying
	2033
	2034	j = j+1
	2035	IF( ABS(delta).GT.xlim ) CYCLE
	2036
	2037	! Determine the Gaussian weights
	2038
	2039	c = 2.0 ( 1.0-cosc(i)cosc(i) )
	2040	CALL lgord( d, cosc(i), nlat-1 )
	2041	d = ddfi*fi
	2042	gwt(i) = c *( fi-0.5 ) / d
	2043	EXIT
	2044
	2045	END DO
	2046
	2047	END DO
	2048
	2049	! Determine the colatitudes and sin(colat) and weights over sin**2
	2050
	2051	DO i=1,nzero
	2052	colat(i)= ACOS(cosc(i))
	2053	sinc(i) = SIN(colat(i))
	2054	wos2(i) = gwt(i) /( sinc(i)*sinc(i) )
	2055	END DO
	2056
	2057	! If NLAT is odd, set values at the equator
	2058
	2059	IF( MOD(nlat,2) .NE. 0 ) THEN
	2060	i = nzero+1
	2061	cosc(i) = 0.0
	2062	c = 2.0
	2063	CALL lgord( d, cosc(i), nlat-1 )
	2064	d = ddfi*fi
	2065	gwt(i) = c *( fi-0.5 ) / d
	2066	colat(i)= pi*0.5
	2067	sinc(i) = 1.0
	2068	wos2(i) = gwt(i)
	2069	END IF
	2070
	2071	! Determine the southern hemisphere values by symmetry
	2072
	2073	DO i=nlat-nzero+1,nlat
	2074	cosc(i) =-cosc(nlat+1-i)
	2075	gwt(i) = gwt(nlat+1-i)
	2076	colat(i)= pi-colat(nlat+1-i)
	2077	sinc(i) = sinc(nlat+1-i)
	2078	wos2(i) = wos2(nlat+1-i)
	2079	END DO
	2080
	2081	END SUBROUTINE lggaus
	2082
	2083
	2084	SUBROUTINE lgord( f, cosc, n )
	2085
	2086	IMPLICIT NONE
	2087
	2088	! LGORD calculates the value of an ordinary Legendre polynomial at a
	2089	! specific latitude.
	2090
	2091	! On entry:
	2092	! cosc - COS(colatitude)
	2093	! n - the degree of the polynomial
	2094
	2095	! On exit:
	2096	! f - the value of the Legendre polynomial of degree N at
	2097	! latitude ASIN(cosc)
	2098
	2099	REAL (KIND=HIGH) :: s1, c4, a, b, fk, f, cosc, colat, c1, fn, ang
	2100	INTEGER :: n, k
	2101
	2102	! Determine the colatitude
	2103
	2104	colat = ACOS(cosc)
	2105
	2106	c1 = SQRT(2.0_HIGH)
	2107	DO k=1,n
	2108	c1 = c1 * SQRT( 1.0 - 1.0/(4kk) )
	2109	END DO
	2110
	2111	fn = n
	2112	ang= fn * colat
	2113	s1 = 0.0
	2114	c4 = 1.0
	2115	a =-1.0
	2116	b = 0.0
	2117	DO k=0,n,2
	2118	IF (k.eq.n) c4 = 0.5 * c4
	2119	s1 = s1 + c4 * COS(ang)
	2120	a = a + 2.0
	2121	b = b + 1.0
	2122	fk = k
	2123	ang= colat * (fn-fk-2.0)
	2124	c4 = ( a * (fn-b+1.0) / ( b * (fn+fn-a) ) ) * c4
	2125	END DO
	2126
	2127	f = s1 * c1
	2128
	2129	END SUBROUTINE lgord
	2130
	2131
	2132	SUBROUTINE llij_gauss (lat, lon, proj, i, j)
	2133
	2134	IMPLICIT NONE
	2135
	2136	REAL , INTENT(IN) :: lat, lon
	2137	REAL , INTENT(OUT) :: i, j
	2138	TYPE (proj_info), INTENT(IN) :: proj
	2139
	2140	INTEGER :: n , n_low
	2141	LOGICAL :: found = .FALSE.
	2142	REAL :: diff_1 , diff_nlat
	2143
	2144	! The easy one first, get the x location. The calling routine has already made
	2145	! sure that the necessary assumptions concerning the sign of the deltalon and the
	2146	! relative east/west'ness of the longitude and the starting longitude are consistent
	2147	! to allow this easy computation.
	2148
	2149	i = ( lon - proj%lon1 ) / proj%loninc + 1.
	2150
	2151	! Since this is a global data set, we need to be concerned about wrapping the
	2152	! fields around the globe.
	2153
	2154	! IF ( ( proj%loninc .GT. 0 ) .AND. &
	2155	! ( FLOOR((lon-proj%lon1)/proj%loninc) + 1 .GE. proj%ixdim ) .AND. &
	2156	! ( lon + proj%loninc .GE. proj%lon1 + 360 ) ) THEN
	2157	!! BUG: We may need to set proj%wrap, but proj is intent(in)
	2158	!! WHAT IS THIS USED FOR?
	2159	!! proj%wrap = .TRUE.
	2160	! i = proj%ixdim
	2161	! ELSE IF ( ( proj%loninc .LT. 0 ) .AND. &
	2162	! ( FLOOR((lon-proj%lon1)/proj%loninc) + 1 .GE. proj%ixdim ) .AND. &
	2163	! ( lon + proj%loninc .LE. proj%lon1 - 360 ) ) THEN
	2164	! ! BUG: We may need to set proj%wrap, but proj is intent(in)
	2165	! ! WHAT IS THIS USED FOR?
	2166	! ! proj%wrap = .TRUE.
	2167	! i = proj%ixdim
	2168	! END IF
	2169
	2170	! Yet another quicky test, can we find bounding values? If not, then we may be
	2171	! dealing with putting data to a polar projection, so just give them them maximal
	2172	! value for the location. This is an OK assumption for the interpolation across the
	2173	! top of the pole, given how close the longitude lines are.
	2174
	2175	IF ( ABS(lat) .GT. ABS(proj%gauss_lat(1)) ) THEN
	2176
	2177	diff_1 = lat - proj%gauss_lat(1)
	2178	diff_nlat = lat - proj%gauss_lat(proj%nlat*2)
	2179
	2180	IF ( ABS(diff_1) .LT. ABS(diff_nlat) ) THEN
	2181	j = 1
	2182	ELSE
	2183	j = proj%nlat*2
	2184	END IF
	2185
	2186	! If the latitude is between the two bounding values, we have to search and interpolate.
	2187
	2188	ELSE
	2189
	2190	DO n = 1 , proj%nlat*2 -1
	2191	IF ( ( proj%gauss_lat(n) - lat ) * ( proj%gauss_lat(n+1) - lat ) .LE. 0 ) THEN
	2192	found = .TRUE.
	2193	n_low = n
	2194	EXIT
	2195	END IF
	2196	END DO
	2197
	2198	! Everything still OK?
	2199
	2200	IF ( .NOT. found ) THEN
	2201	PRINT '(A)','Troubles in river city. No bounding values of latitude found in the Gaussian routines.'
	2202	CALL wrf_error_fatal ( 'Gee_no_bounding_lats_Gaussian' )
	2203	END IF
	2204
	2205	j = ( ( proj%gauss_lat(n_low) - lat ) * ( n_low + 1 ) + &
	2206	( lat - proj%gauss_lat(n_low+1) ) * ( n_low ) ) / &
	2207	( proj%gauss_lat(n_low) - proj%gauss_lat(n_low+1) )
	2208
	2209	END IF
	2210
	2211	END SUBROUTINE llij_gauss
	2212
	2213	END MODULE module_llxy

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format