source: branches/jerome/src_flexwrf_v3.1/map_proj_wrf_subaa.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago

sources for flexwrf v3.1

File size: 35.8 KB
Line 
1! SUPPORTED PROJECTIONS
2! ---------------------
3! Cylindrical Lat/Lon (code = PROJ_LATLON)
4! Mercator (code = PROJ_MERC)
5! Lambert Conformal (code = PROJ_LC)
6! Polar Stereographic (code = PROJ_PS)
7!
8! REMARKS
9! -------
10! The routines contained within were adapted from routines
11! obtained from the NCEP w3 library.  The original NCEP routines were
12! less
13! flexible (e.g., polar-stereo routines only supported truelat of
14! 60N/60S)
15! than what we needed, so modifications based on equations in Hoke,
16! Hayes, and
17! Renninger (AFGWC/TN/79-003) were added to improve the flexibility. 
18! Additionally, coding was improved to F90 standards and the routines
19! were
20! combined into this module. 
21!-----------------------------------------------------------------------
22
23
24!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
25!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
26!dis   
27!dis    open source license/disclaimer, forecast systems laboratory
28!dis    noaa/oar/fsl, 325 broadway boulder, co 80305
29!dis   
30!dis    this software is distributed under the open source definition,
31!dis    which may be found at http://www.opensource.org/osd.html.
32!dis   
33!dis    in particular, redistribution and use in source and binary forms,
34!dis    with or without modification, are permitted provided that the
35!dis    following conditions are met:
36!dis   
37!dis    - redistributions of source code must retain this notice, this
38!dis    list of conditions and the following disclaimer.
39!dis   
40!dis    - redistributions in binary form must provide access to this
41!dis    notice, this list of conditions and the following disclaimer, and
42!dis    the underlying source code.
43!dis   
44!dis    - all modifications to this software must be clearly documented,
45!dis    and are solely the responsibility of the agent making the
46!dis    modifications.
47!dis   
48!dis    - if significant modifications or enhancements are made to this
49!dis    software, the fsl software policy manager
50!dis    (softwaremgr@fsl.noaa.gov) should be notified.
51!dis   
52!dis    this software and its documentation are in the public domain
53!dis    and are furnished "as is."  the authors, the united states
54!dis    government, its instrumentalities, officers, employees, and
55!dis    agents make no warranty, express or implied, as to the usefulness
56!dis    of the software and documentation for any purpose.  they assume
57!dis    no responsibility (1) for the use of the software and
58!dis    documentation; or (2) to provide technical support to users.
59!dis   
60!dis
61
62! module that defines constants, data structures, and
63! routines used to convert grid indices to lat/lon
64! and vice versa.   
65!
66! supported projections
67! ---------------------
68! cylindrical lat/lon (code = proj_latlon)
69! mercator (code = proj_merc)
70! lambert conformal (code = proj_lc)
71! polar stereographic (code = proj_ps)
72!
73! remarks
74! -------
75! the routines contained within were adapted from routines
76! obtained from the ncep w3 library.  the original ncep routines were less
77! flexible (e.g., polar-stereo routines only supported truelat of 60n/60s)
78! than what we needed, so modifications based on equations in hoke, hayes, and
79! renninger (afgwc/tn/79-003) were added to improve the flexibility. 
80! additionally, coding was improved to f90 standards and the routines were
81! combined into this module. 
82!
83! assumptions
84! -----------
85!  grid definition:
86!    for mercator, lambert conformal, and polar-stereographic projections,
87!    the routines within assume the following:
88!
89!       1.  grid is dimensioned (i,j) where i is the east-west direction,
90!           positive toward the east, and j is the north-south direction,
91!           positive toward the north. 
92!       2.  origin is at (1,1) and is located at the southwest corner,
93!           regardless of hemispere.
94!       3.  grid spacing (dx) is always positive.
95!       4.  values of true latitudes must be positive for nh domains
96!           and negative for sh domains.
97!
98!     for the latlon projection, the grid origin may be at any of the
99!     corners, and the deltalat and deltalon values can be signed to
100!     account for this using the following convention:
101!       origin location        deltalat sign      deltalon sign
102!       ---------------        -------------      -------------
103!        sw corner                  +                   +
104!        ne corner                  -                   -
105!        nw corner                  -                   +
106!        se corner                  +                   -
107!       
108!  data definitions:
109!       1. any arguments that are a latitude value are expressed in
110!          degrees north with a valid range of -90 -> 90
111!       2. any arguments that are a longitude value are expressed in
112!          degrees east with a valid range of -180 -> 180.
113!       3. distances are in meters and are always positive.
114!       4. the standard longitude (stdlon) is defined as the longitude
115!          line which is parallel to the y-axis (j-direction), along
116!          which latitude increases (not the absolute value of latitude, but
117!          the actual latitude, such that latitude increases continuously
118!          from the south pole to the north pole) as j increases. 
119!       5. one true latitude value is required for polar-stereographic and
120!          mercator projections, and defines at which latitude the
121!          grid spacing is true.  for lambert conformal, two true latitude
122!          values must be specified, but may be set equal to each other to
123!          specify a tangent projection instead of a secant projection. 
124!       
125! usage
126! -----
127! to use the routines in this module, the calling routines must have the
128! following statement at the beginning of its declaration block:
129!   use map_utils
130!
131! the use of the module not only provides access to the necessary routines,
132! but also defines a structure of type (proj_info) that can be used
133! to declare a variable of the same type to hold your map projection
134! information.  it also defines some integer parameters that contain
135! the projection codes so one only has to use those variable names rather
136! than remembering the acutal code when using them.  the basic steps are
137! as follows:
138
139!   1.  ensure the "use map_utils" is in your declarations.
140!   2.  declare the projection information structure as type(proj_info):
141!         type(proj_info) :: proj
142!   3.  populate your structure by calling the map_set routine:
143!         call map_set(code,lat1,lon1,dx,stdlon,truelat1,truelat2,nx,ny,proj)
144!       where:
145!         code (input) = one of proj_latlon, proj_merc, proj_lc, or proj_ps
146!         lat1 (input) = latitude of grid origin point (i,j)=(1,1)
147!                         (see assumptions!)
148!         lon1 (input) = longitude of grid origin
149!         dx (input) = grid spacing in meters (ignored for latlon projections)
150!         stdlon (input) = standard longitude for proj_ps and proj_lc,
151!               deltalon (see assumptions) for proj_latlon,
152!               ignored for proj_merc
153!         truelat1 (input) = 1st true latitude for proj_ps, proj_lc, and
154!                proj_merc, deltalat (see assumptions) for proj_latlon
155!         truelat2 (input) = 2nd true latitude for proj_lc,
156!                ignored for all others.
157!         nx = number of points in east-west direction
158!         ny = number of points in north-south direction
159!         proj (output) = the structure of type (proj_info) that will be fully
160!                populated after this call
161!
162!   4.  now that the proj structure is populated, you may call any
163!       of the following routines:
164!       
165!       latlon_to_ij(proj, lat, lon, i, j)
166!       ij_to_latlon(proj, i, j, lat, lon)
167!       truewind_to_gridwind(lon, proj, ugrid, vgrid, utrue, vtrue)
168!       gridwind_to_truewind(lon, proj, utrue, vtrue, ugrid, vgrid)
169!       compare_projections(proj1, proj2, same_proj)
170!
171!       it is incumbent upon the calling routine to determine whether or
172!       not the values returned are within your domain bounds.  all values
173!       of i, j, lat, and lon are real values.
174!
175!
176! references
177! ----------
178!  hoke, hayes, and renninger, "map preojections and grid systems for
179!       meteorological applications." afgwc/tn-79/003(rev), air weather
180!       service, 1985.
181!
182!  ncar mm5v3 modeling system, regridder program, module_first_guess_map.f
183!  ncep routines w3fb06, w3fb07, w3fb08, w3fb09, w3fb11, w3fb12
184!
185! history
186! -------
187! 27 mar 2001 - original version
188!               brent l. shaw, noaa/fsl (csu/cira)
189! 02 apr 2001 - added routines to rotate winds from true to grid
190!               and vice versa.
191!               brent l. shaw, noaa/fsl (csu/cira)
192! 09 apr 2001 - added compare_projections routine to compare two
193!               sets of projection parameters.
194!
195!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
196!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197
198
199
200   
201!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
202! ! define data structures to define various projections
203!
204! type proj_info
205!
206!   logical        proj_init     ! flag to indicate if this struct is
207!                                ! ready for use
208!   logical        proj_cyclic   ! flag indicating if this grid
209!                                ! is cyclic in the longitudinal
210!                                ! direction...happens with
211!                                ! global lat/lon grids like gfs/avn
212!   integer        proj_code     ! integer code for projection type
213!   integer        proj_nx
214!   integer        proj_ny
215!   real           proj_lat1    ! sw latitude (1,1) in degrees (-90->90n)
216!   real           proj_lon1    ! sw longitude (1,1) in degrees (-180->180e)
217!   real           proj_dx       ! grid spacing in meters at truelats, used
218!                                ! only for ps, lc, and merc projections
219!   real           proj_dlat     ! lat increment for lat/lon grids
220!   real           proj_dlon     ! lon increment for lat/lon grids
221!   real           proj_clat     ! center latitude of grid
222!   real           proj_clon     ! center longitude of grid
223!   real           proj_stdlon   ! longitude parallel to y-axis (-180->180e)
224!   real           proj_truelat1 ! first true latitude (all projections)
225!   real           proj_truelat2 ! second true lat (lc only)
226!   real           proj_hemi     ! 1 for nh, -1 for sh
227!   real           proj_cone     ! cone factor for lc projections
228!   real           proj_polei    ! computed i-location of pole point
229!   real           proj_polej    ! computed j-location of pole point
230!   real           proj_rsw      ! computed radius to sw corner
231!   real           proj_rebydx   ! earth radius divided by dx
232!
233! end type proj_info
234!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
235
236
237!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238      subroutine map_init
239    ! initializes the map projection structure to missing values
240
241   use wrf_map_utils_mod
242      implicit none
243
244!      include 'include_wrf_map_utils'
245      pi = 3.1415927
246      deg_per_rad = 180./pi
247      rad_per_deg = pi / 180.
248
249! mean earth radius in m.  the value below is consistent
250! with nceps routines and grids.
251      earth_radius_m = 6370000.
252!     earth_radius_m = 6371200.
253!     earth_radius_m = 6368750.
254
255      proj_latlon = 0
256      proj_merc = 1
257      proj_lc = 3
258      proj_ps = 5
259      proj_rotlat = 203
260
261      proj_lat1 =    -999.9
262      proj_lon1 =    -999.9
263      proj_dx    =    -999.9
264      proj_stdlon =   -999.9
265      proj_truelat1 = -999.9
266      proj_truelat2 = -999.9
267      proj_hemi     = 0.0
268      proj_cone     = -999.9
269      proj_polei    = -999.9
270      proj_polej    = -999.9
271      proj_rsw      = -999.9
272      proj_init     = .false.
273      proj_nx       = -99
274      proj_ny       = -99
275      proj_cyclic   = .false.
276
277      return
278      end subroutine map_init
279
280
281!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
282      subroutine map_set( &
283          proj_code_in, lat1, lon1, dx, &
284          stdlon, truelat1, truelat2, &
285          idim, jdim )
286      ! given a partially filled proj_info structure, this routine computes
287      ! polei, polej, rsw, and cone (if lc projection) to complete the
288      ! structure.  this allows us to eliminate redundant calculations when
289      ! calling the coordinate conversion routines multiple times for the
290      ! same map.
291      ! this will generally be the first routine called when a user wants
292      ! to be able to use the coordinate conversion routines, and it
293      ! will call the appropriate routines based on the
294      ! proj_code which indicates which projection type  this is.
295   use wrf_map_utils_mod
296
297      implicit none
298     
299      ! declare arguments
300      integer :: proj_code_in
301      real :: lat1
302      real :: lon1
303      real :: dx
304      real :: stdlon
305      real :: truelat1
306      real :: truelat2
307      integer :: idim
308      integer :: jdim
309
310      ! local variables
311      real :: center_i,center_j
312      real :: center_lat, center_lon
313
314!      include 'include_wrf_map_utils'
315
316      ! executable code
317
318      proj_code = proj_code_in
319
320      ! first, check for validity of mandatory variables in proj
321      if ( abs(lat1) .gt. 90.001 ) then
322        print '(a)', 'latitude of origin corner required as follows:'
323        print '(a)', '    -90n <= lat1 < = 90.n'
324        stop 'map_init'
325      endif
326      if ( abs(lon1) .gt. 180.) then
327        print '(a)', 'longitude of origin required as follows:'
328        print '(a)', '   -180e <= lon1 <= 180w'
329        stop 'map_init'
330      endif
331      if ((dx .le. 0.).and.(proj_code .ne. proj_latlon)) then
332        print '(a)', 'require grid spacing (dx) in meters be positive!'
333        stop 'map_init'
334      endif
335      if ((abs(stdlon) .gt. 180.).and.(proj_code .ne. proj_merc)) then
336        print '(a)', 'need orientation longitude (stdlon) as: '
337        print '(a)', '   -180e <= lon1 <= 180w'
338        stop 'map_init'
339      endif
340      if (abs(truelat1).gt.90.) then
341        print '(a)', 'set true latitude 1 for all projections!'
342        stop 'map_init'
343      endif
344     
345      call map_init       
346      proj_code  = proj_code_in
347      proj_lat1 = lat1
348      proj_lon1 = lon1
349      proj_dx    = dx
350      proj_stdlon = stdlon
351      proj_truelat1 = truelat1
352      proj_truelat2 = truelat2
353      proj_nx = idim
354      proj_ny = jdim
355      if (proj_code .ne. proj_latlon) then
356        proj_dx = dx
357        if (truelat1 .lt. 0.) then
358          proj_hemi = -1.0
359        else
360          proj_hemi = 1.0
361        endif
362        proj_rebydx = earth_radius_m / dx
363      endif
364
365
366      if (proj_code .eq. proj_ps) then
367          !print '(a)', 'setting up polar stereographic map...'
368          call set_ps
369
370      else if (proj_code .eq. proj_lc) then
371          !print '(a)', 'setting up lambert conformal map...'
372          if (abs(proj_truelat2) .gt. 90.) then
373            print '(a)',  &
374                'second true latitude not set, assuming a tangent'
375            print '(a,f10.3)', 'projection at truelat1: ', proj_truelat1
376            proj_truelat2=proj_truelat1
377          else
378            ! ensure truelat1 < truelat2
379            proj_truelat1 = min(truelat1,truelat2)
380            proj_truelat2 = max(truelat1,truelat2)
381          endif
382          call set_lc
383     
384      else if (proj_code .eq. proj_merc) then
385          !print '(a)', 'setting up mercator map...'
386          call set_merc
387     
388      else if (proj_code .eq. proj_latlon) then
389          !print '(a)', 'setting up cylindrical equidistant latlon map...'
390          ! convert lon1 to 0->360 notation
391          if (proj_lon1 .lt. 0.) proj_lon1 = proj_lon1 + 360.
392          proj_dlat = truelat1
393          proj_dlon = stdlon
394          if (nint(proj_dlon*real(proj_nx)) .eq. 360) proj_cyclic = .true.
395
396      else
397          print '(a,i2)', 'unknown projection code: ', proj_code
398          stop 'map_init'
399     
400      end if
401
402      proj_init = .true.
403
404      center_i = real(proj_nx+1)*0.5
405      center_j = real(proj_ny+1)*0.5
406      call ij_to_latlon(center_i,center_j,center_lat,center_lon)
407      proj_clat = center_lat
408      proj_clon = center_lon
409
410      return
411      end subroutine map_set
412
413
414!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
415      subroutine latlon_to_ij(lat, lon, i, j)
416      ! converts input lat/lon values to the cartesian (i,j) value
417      ! for the given projection.
418
419   use wrf_map_utils_mod
420      implicit none
421
422      ! arguments
423      real :: lat
424      real :: lon
425      real :: i
426      real :: j
427
428!      include 'include_wrf_map_utils'
429
430      if (.not.proj_init) then
431        print '(a)', 'you have not called map_set for this projection!'
432        stop 'latlon_to_ij'
433      endif
434
435      if (proj_code .eq. proj_ps) then
436          call llij_ps(lat,lon,i,j)
437       
438      else if (proj_code .eq. proj_lc) then
439          call llij_lc(lat,lon,i,j)
440 
441      else if (proj_code .eq. proj_latlon) then
442          call llij_latlon(lat,lon,i,j)
443!
444      else if (proj_code .eq. proj_merc) then
445          call llij_merc(lat,lon,i,j)
446!
447!     else if (proj_code .eq. proj_rotlat) then
448!       write(6,*)   'doing nothing in latlon_to_ij'
449
450      else
451          print '(a,i2)','unrecognized map projection code: ', proj_code
452          stop 'latlon_to_ij'
453   
454      end if
455
456      return
457      end subroutine latlon_to_ij
458
459
460!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
461      subroutine ij_to_latlon(i, j, lat, lon)
462      ! computes geographical latitude and longitude for a given (i,j) point
463      ! in a grid with a projection of proj
464
465   use wrf_map_utils_mod
466      implicit none
467
468      ! arguments
469      real :: i
470      real :: j
471      real :: lat
472      real :: lon
473
474!      include 'include_wrf_map_utils'
475
476      if (.not.proj_init) then
477        print '(a)', 'you have not called map_set for this projection!'
478        stop 'ij_to_latlon'
479      endif
480
481      if (proj_code .eq. proj_ps) then
482          call ijll_ps(i, j, lat, lon)
483
484      else if (proj_code .eq. proj_lc) then
485          call ijll_lc(i, j, lat, lon)
486 
487      else if (proj_code .eq. proj_latlon) then
488          call ijll_latlon(i, j, lat, lon)
489!
490      else if (proj_code .eq. proj_merc) then
491          call ijll_merc(i, j, lat, lon)
492!
493!     else if (proj_code .eq. proj_rotlat) then
494!       write(6,*)   'doing nothing in ij_to_latlon'
495
496      else
497          print '(a,i2)','unrecognized map projection code: ', proj_code
498          stop 'ij_to_latlon'
499
500      end if
501
502      return
503      end subroutine ij_to_latlon
504
505
506!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
507      subroutine set_ps
508      ! initializes a polar-stereographic map projection from the partially
509      ! filled proj structure. this routine computes the radius to the
510      ! southwest corner and computes the i/j location of the pole for use
511      ! in llij_ps and ijll_ps.
512
513   use wrf_map_utils_mod
514      implicit none
515   
516      ! local vars
517      real :: ala1
518      real :: alo1
519      real :: reflon
520      real :: scale_top
521      real :: dumlat, dumlon
522
523!      include 'include_wrf_map_utils'
524
525      ! executable code
526      reflon = proj_stdlon + 90.
527   
528      ! cone factor
529      proj_cone = 1.0
530
531      ! compute numerator term of map scale factor
532      scale_top = 1. + proj_hemi * sin(proj_truelat1 * rad_per_deg)
533
534      ! compute radius to lower-left (sw) corner
535      ala1 = proj_lat1 * rad_per_deg
536      proj_rsw =proj_rebydx*cos(ala1)*scale_top/(1.+proj_hemi*sin(ala1))
537
538      ! find the pole point
539      alo1 = (proj_lon1 - reflon) * rad_per_deg
540      proj_polei = 1. - proj_rsw * cos(alo1)
541      proj_polej = 1. - proj_hemi * proj_rsw * sin(alo1)
542!     print '(a,2f14.5)', 'set_ps - computed (i,j) of pole point: ',
543!    &    proj_polei,proj_polej
544
545      return
546      end subroutine set_ps
547
548
549!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
550      subroutine llij_ps( lat, lon, i, j )
551      ! given latitude (-90 to 90), longitude (-180 to 180), and the
552      ! standard polar-stereographic projection information via the
553      ! public proj structure, this routine returns the i/j indices which
554      ! if within the domain range from 1->nx and 1->ny, respectively.
555   use wrf_map_utils_mod
556
557      implicit none
558
559      ! delcare input arguments
560      real :: lat
561      real :: lon
562
563      ! declare output arguments     
564      real :: i !(x-index)
565      real :: j !(y-index)
566
567      ! declare local variables
568      real :: reflon
569      real :: scale_top
570      real :: ala
571      real :: alo
572      real :: rm
573
574!      include 'include_wrf_map_utils'
575
576      ! begin code
577   
578      reflon = proj_stdlon + 90.
579     
580      ! compute numerator term of map scale factor
581
582      scale_top = 1. + proj_hemi * sin(proj_truelat1 * rad_per_deg)
583
584      ! find radius to desired point
585      ala = lat * rad_per_deg
586      rm = proj_rebydx * cos(ala) * scale_top/(1. + proj_hemi *sin(ala))
587      alo = (lon - reflon) * rad_per_deg
588      i = proj_polei + rm * cos(alo)
589      j = proj_polej + proj_hemi * rm * sin(alo)
590   
591      return
592      end subroutine llij_ps
593
594
595!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
596      subroutine ijll_ps( i, j, lat, lon )
597
598      ! this is the inverse routine of llij_ps.  it returns the
599      ! latitude and longitude of an i/j point given the projection info
600      ! structure. 
601   use wrf_map_utils_mod
602
603      implicit none
604
605      ! declare input arguments
606      real ::  i    ! column
607      real :: j    ! row
608     
609      ! declare output arguments
610      real :: lat     ! -90 -> 90 north
611      real :: lon     ! -180 -> 180 east
612
613      ! local variables
614      real :: reflon
615      real :: scale_top
616      real :: xx,yy
617      real :: gi2, r2
618      real :: arccos
619
620!      include 'include_wrf_map_utils'
621
622      ! begin code
623
624      ! compute the reference longitude by rotating 90 degrees to the east
625      ! to find the longitude line parallel to the positive x-axis.
626      reflon = proj_stdlon + 90.
627     
628      ! compute numerator term of map scale factor
629      scale_top = 1. + proj_hemi * sin(proj_truelat1 * rad_per_deg)
630
631      ! compute radius to point of interest
632      xx = i - proj_polei
633      yy = (j - proj_polej) * proj_hemi
634      r2 = xx**2 + yy**2
635
636      ! now the magic code
637      if (r2 .eq. 0.) then
638        lat = proj_hemi * 90.
639        lon = reflon
640      else
641        gi2 = (proj_rebydx * scale_top)**2.
642        lat = deg_per_rad * proj_hemi * asin((gi2-r2)/(gi2+r2))
643        arccos = acos(xx/sqrt(r2))
644        if (yy .gt. 0) then
645          lon = reflon + deg_per_rad * arccos
646        else
647          lon = reflon - deg_per_rad * arccos
648        endif
649      endif
650   
651      ! convert to a -180 -> 180 east convention
652      if (lon .gt. 180.) lon = lon - 360.
653      if (lon .lt. -180.) lon = lon + 360.
654
655      return
656      end subroutine ijll_ps
657
658
659!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
660      subroutine set_lc     
661      ! initialize the remaining items in the proj structure for a
662      ! lambert conformal grid.
663
664   use wrf_map_utils_mod
665      implicit none
666
667!      include 'include_wrf_map_utils'
668     
669      real :: arg
670      real :: deltalon1
671      real :: tl1r
672      real :: ctl1r
673
674      ! compute cone factor
675      call lc_cone(proj_truelat1, proj_truelat2, proj_cone)
676      ! print '(a,f8.6)', 'computed cone factor: ', proj_cone
677      ! compute longitude differences and ensure we stay out of the
678      ! forbidden "cut zone"
679      deltalon1 = proj_lon1 - proj_stdlon
680      if (deltalon1 .gt. +180.) deltalon1 = deltalon1 - 360.
681      if (deltalon1 .lt. -180.) deltalon1 = deltalon1 + 360.
682
683      ! convert truelat1 to radian and compute cos for later use
684      tl1r = proj_truelat1 * rad_per_deg
685      ctl1r = cos(tl1r)
686
687      ! compute the radius to our known lower-left (sw) corner
688      proj_rsw = proj_rebydx * ctl1r/proj_cone *  &
689             (tan((90.*proj_hemi-proj_lat1)*rad_per_deg/2.) / &
690              tan((90.*proj_hemi-proj_truelat1)*rad_per_deg/2.)) &
691              **proj_cone
692
693      ! find pole point
694      arg = proj_cone*(deltalon1*rad_per_deg)
695      proj_polei = 1. - proj_hemi * proj_rsw * sin(arg)
696      proj_polej = 1. + proj_rsw * cos(arg) 
697      !print '(a,2f10.3)', 'computed pole i/j = ', proj_polei, proj_polej
698
699      return
700      end subroutine set_lc
701
702
703!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
704      subroutine lc_cone(truelat1, truelat2, cone)
705
706    ! routine to compute the cone factor of a lambert conformal projection
707
708   use wrf_map_utils_mod
709      implicit none
710
711!      include 'include_wrf_map_utils'
712     
713      ! input args
714      real :: truelat1  ! (-90 -> 90 degrees n)
715      real :: truelat2  !   "   "  "   "     "
716
717      ! output args
718      real :: cone
719
720      ! locals
721
722      ! begin code
723
724      ! first, see if this is a secant or tangent projection.  for tangent
725      ! projections, truelat1 = truelat2 and the cone is tangent to the
726      ! earth surface at this latitude.  for secant projections, the cone
727      ! intersects the earth surface at each of the distinctly different
728      ! latitudes
729      if (abs(truelat1-truelat2) .gt. 0.1) then
730
731        ! compute cone factor following:
732        cone=(alog(cos(truelat1*rad_per_deg)) &
733             -alog(cos(truelat2*rad_per_deg))) / &
734         (alog(tan((90.-abs(truelat1))*rad_per_deg*0.5 ))- &
735          alog(tan((90.-abs(truelat2))*rad_per_deg*0.5 )) )
736      else
737         cone = sin(abs(truelat1)*rad_per_deg ) 
738      endif
739
740      return
741      end subroutine lc_cone
742
743
744!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
745      subroutine ijll_lc( i, j, lat, lon)
746
747    ! routine to convert from the (i,j) cartesian coordinate to the
748    ! geographical latitude and longitude for a lambert conformal projection.
749
750    ! history:
751    ! 25 jul 01: corrected by b. shaw, noaa/fsl
752    !
753   use wrf_map_utils_mod
754      implicit none
755
756!      include 'include_wrf_map_utils'
757
758      ! input args
759      real :: i        ! cartesian x coordinate
760      real :: j        ! cartesian y coordinate
761
762      ! output args                 
763      real :: lat      ! latitude (-90->90 deg n)
764      real :: lon      ! longitude (-180->180 e)
765
766      ! locals
767      real :: inew
768      real :: jnew
769      real :: r
770      real :: chi,chi1,chi2
771      real :: r2
772      real :: xx
773      real :: yy
774
775      ! begin code
776
777      chi1 = (90. - proj_hemi*proj_truelat1)*rad_per_deg
778      chi2 = (90. - proj_hemi*proj_truelat2)*rad_per_deg
779
780      ! see if we are in the southern hemispere and flip the indices
781      ! if we are.
782      if (proj_hemi .eq. -1.) then
783        inew = -i + 2.
784        jnew = -j + 2.
785      else
786        inew = i
787        jnew = j
788      endif
789
790      ! compute radius**2 to i/j location
791      xx = inew - proj_polei
792      yy = proj_polej - jnew
793      r2 = (xx*xx + yy*yy)
794      r = sqrt(r2)/proj_rebydx
795     
796      ! convert to lat/lon
797      if (r2 .eq. 0.) then
798        lat = proj_hemi * 90.
799        lon = proj_stdlon
800      else
801         
802        ! longitude
803        lon = proj_stdlon + deg_per_rad * atan2(proj_hemi*xx,yy) &
804              /proj_cone
805        lon = amod(lon+360., 360.)
806
807        ! latitude.  latitude determined by solving an equation adapted
808        ! from:
809        !  maling, d.h., 1973: coordinate systems and map projections
810        ! equations #20 in appendix i. 
811         
812        if (chi1 .eq. chi2) then
813          chi = 2.0*atan( ( r/tan(chi1) )**(1./proj_cone)  &
814                * tan(chi1*0.5) )
815        else
816          chi = 2.0*atan( (r*proj_cone/sin(chi1))**(1./proj_cone)  &
817                * tan(chi1*0.5))
818        endif
819        lat = (90.0-chi*deg_per_rad)*proj_hemi
820
821      endif
822
823      if (lon .gt. +180.) lon = lon - 360.
824      if (lon .lt. -180.) lon = lon + 360.
825
826      return
827      end subroutine ijll_lc
828
829
830!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
831      subroutine llij_lc( lat, lon, i, j)
832
833    ! routine to compute the geographical latitude and longitude values
834    ! to the cartesian x/y on a lambert conformal projection.
835     
836    use wrf_map_utils_mod
837      implicit none
838!
839!      include 'include_wrf_map_utils'
840
841      ! input args
842      real :: lat      ! latitude (-90->90 deg n)
843      real :: lon      ! longitude (-180->180 e)
844
845      ! output args                 
846      real :: i        ! cartesian x coordinate
847      real :: j        ! cartesian y coordinate
848
849      ! locals
850      real :: arg
851      real :: deltalon
852      real :: tl1r
853      real :: rm
854      real :: ctl1r
855     
856
857      ! begin code
858!         print*,proj_clat, proj_clon,proj_stdlon,proj_truelat1,proj_code,proj_cone
859      ! compute deltalon between known longitude and standard lon and ensure
860      ! it is not in the cut zone
861      deltalon = lon - proj_stdlon
862      if (deltalon .gt. +180.) deltalon = deltalon - 360.
863      if (deltalon .lt. -180.) deltalon = deltalon + 360.
864     
865      ! convert truelat1 to radian and compute cos for later use
866      tl1r = proj_truelat1 * rad_per_deg
867      ctl1r = cos(tl1r)     
868     
869      ! radius to desired point
870      rm = proj_rebydx * ctl1r/proj_cone *  &
871          (tan((90.*proj_hemi-lat)*rad_per_deg/2.) / &
872           tan((90.*proj_hemi-proj_truelat1)*rad_per_deg/2.))**proj_cone
873
874      arg = proj_cone*(deltalon*rad_per_deg)
875      i = proj_polei + proj_hemi * rm * sin(arg)
876      j = proj_polej - rm * cos(arg)
877
878      ! finally, if we are in the southern hemisphere, flip the i/j
879      ! values to a coordinate system where (1,1) is the sw corner
880      ! (what we assume) which is different than the original ncep
881      ! algorithms which used the ne corner as the origin in the
882      ! southern hemisphere (left-hand vs. right-hand coordinate?)
883      if (proj_hemi .eq. -1.) then
884        i = 2. - i 
885        j = 2. - j
886      endif
887
888      return
889      end subroutine llij_lc
890
891!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
892  SUBROUTINE set_merc
893
894    ! Sets up the remaining basic elements for the mercator projection
895     use wrf_map_utils_mod
896
897    IMPLICIT NONE
898!   TYPE(proj_info), INTENT(INOUT)       :: proj
899    REAL                                 :: clain
900
901    !  Preliminary variables
902
903    clain = COS(rad_per_deg*proj_truelat1)
904    proj_dlon = proj_dx / (earth_radius_m * clain)
905
906    ! Compute distance from equator to origin, and store in the
907    ! proj%rsw tag.
908
909    proj_rsw = 0.
910    IF (proj_lat1 .NE. 0.) THEN
911      proj_rsw = (ALOG(TAN(0.5*((proj_lat1+90.)*rad_per_deg))))/proj_dlon
912    ENDIF
913    RETURN
914  END SUBROUTINE set_merc
915!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
916
917  SUBROUTINE llij_merc(lat, lon, i, j)
918
919    ! Compute i/j coordinate from lat lon for mercator projection
920     use wrf_map_utils_mod
921
922    IMPLICIT NONE
923    REAL, INTENT(IN)              :: lat
924    REAL, INTENT(IN)              :: lon
925!   TYPE(proj_info),INTENT(IN)    :: proj
926    REAL,INTENT(OUT)              :: i
927    REAL,INTENT(OUT)              :: j
928    REAL                          :: deltalon
929
930     deltalon = lon - proj_lon1
931!   deltalon = lon - proj_stdlon
932! JB modif
933!   IF (deltalon .LT. -180.) deltalon = deltalon + 360.
934!   IF (deltalon .GT. 180.) deltalon = deltalon - 360.
935    i = 1. + (deltalon/(proj_dlon*deg_per_rad))
936    j = 1. + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
937           proj_dlon - proj_rsw
938!       print*,'IN MERCATOR',lat,lon,i,j,proj_lon1,deltalon
939!      pause
940    RETURN
941  END SUBROUTINE llij_merc
942!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
943  SUBROUTINE ijll_merc(i, j, lat, lon)
944
945    ! Compute the lat/lon from i/j for mercator projection
946     use wrf_map_utils_mod
947
948    IMPLICIT NONE
949    REAL,INTENT(IN)               :: i
950    REAL,INTENT(IN)               :: j
951!   TYPE(proj_info),INTENT(IN)    :: proj
952    REAL, INTENT(OUT)             :: lat
953    REAL, INTENT(OUT)             :: lon
954
955
956    lat = 2.0*ATAN(EXP(proj_dlon*(proj_rsw + j-1.)))*deg_per_rad - 90.
957    lon = (i-1.)*proj_dlon*deg_per_rad + proj_lon1
958    IF (lon.GT.180.) lon = lon - 360.
959    IF (lon.LT.-180.) lon = lon + 360.
960    RETURN
961  END SUBROUTINE ijll_merc
962!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
963  SUBROUTINE llij_latlon(lat, lon, i, j)
964
965    ! Compute the i/j location of a lat/lon on a LATLON grid.
966     use wrf_map_utils_mod
967
968    IMPLICIT NONE
969    REAL, INTENT(IN)             :: lat
970    REAL, INTENT(IN)             :: lon
971!   TYPE(proj_info), INTENT(IN)  :: proj
972    REAL, INTENT(OUT)            :: i
973    REAL, INTENT(OUT)            :: j
974
975    REAL                         :: deltalat
976    REAL                         :: deltalon
977    REAL                         :: lon360
978    REAL                         :: latinc
979    REAL                         :: loninc
980
981    ! Extract the latitude and longitude increments for this grid
982    ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
983    ! loninc is saved in the stdlon tag and latinc is saved in truelat1
984
985    latinc = proj_truelat1
986    loninc = proj_stdlon
987
988    ! Compute deltalat and deltalon as the difference between the input
989    ! lat/lon and the origin lat/lon
990
991    deltalat = lat - proj_lat1
992    ! To account for issues around the dateline, convert the incoming
993    ! longitudes to be 0->360.
994    IF (lon .LT. 0) THEN
995      lon360 = lon + 360.
996    ELSE
997      lon360 = lon
998    ENDIF
999    deltalon = lon360 - proj_lon1
1000    IF (deltalon .LT. 0) deltalon = deltalon + 360.
1001
1002    ! Compute i/j
1003    i = deltalon/loninc + 1.
1004    j = deltalat/latinc + 1.
1005    RETURN
1006    END SUBROUTINE llij_latlon
1007!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1008  SUBROUTINE ijll_latlon(i, j, lat, lon)
1009
1010    ! Compute the lat/lon location of an i/j on a LATLON grid.
1011     use wrf_map_utils_mod
1012
1013    IMPLICIT NONE
1014    REAL, INTENT(IN)             :: i
1015    REAL, INTENT(IN)             :: j
1016!   TYPE(proj_info), INTENT(IN)  :: proj
1017    REAL, INTENT(OUT)            :: lat
1018    REAL, INTENT(OUT)            :: lon
1019
1020    REAL                         :: deltalat
1021    REAL                         :: deltalon
1022    REAL                         :: lon360
1023    REAL                         :: latinc
1024    REAL                         :: loninc
1025
1026    ! Extract the latitude and longitude increments for this grid
1027    ! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
1028    ! loninc is saved in the stdlon tag and latinc is saved in truelat1
1029
1030    latinc = proj_truelat1
1031    loninc = proj_stdlon
1032
1033    ! Compute deltalat and deltalon
1034
1035    deltalat = (j-1.)*latinc
1036    deltalon = (i-1.)*loninc
1037    lat = proj_lat1 + deltalat
1038    lon = proj_lon1 + deltalon
1039
1040    IF ((ABS(lat) .GT. 90.).OR.(ABS(deltalon) .GT.360.)) THEN
1041      ! Off the earth for this grid
1042      lat = -999.
1043      lon = -999.
1044    ELSE
1045      lon = lon + 360.
1046      lon = AMOD(lon,360.)
1047      IF (lon .GT. 180.) lon = lon -360.
1048    ENDIF
1049
1050    RETURN
1051  END SUBROUTINE ijll_latlon
1052!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1053  SUBROUTINE gridwind_to_truewind(lon,ugrid,vgrid,utrue,vtrue)
1054
1055    ! Subroutine to convert a wind from grid north to true north.
1056     use wrf_map_utils_mod
1057
1058    IMPLICIT NONE
1059
1060    ! Arguments
1061    REAL, INTENT(IN)          :: lon     ! Longitude of point in degrees
1062!   TYPE(proj_info),INTENT(IN):: proj    ! Projection info structure
1063    REAL, INTENT(IN)          :: ugrid   ! U-component, grid-relative
1064    REAL, INTENT(IN)          :: vgrid   ! V-component, grid-relative
1065    REAL, INTENT(OUT)         :: utrue   ! U-component, earth-relative
1066    REAL, INTENT(OUT)         :: vtrue   ! V-component, earth-relative
1067
1068    ! Locals
1069    REAL                      :: alpha
1070    REAL                      :: diff
1071
1072    IF ((proj_code .EQ. PROJ_PS).OR.(proj_code .EQ. PROJ_LC))THEN
1073      diff = lon - proj_stdlon
1074      IF (diff .GT. 180.) diff = diff - 360.
1075      IF (diff .LT.-180.) diff = diff + 360.
1076      alpha = diff * proj_cone * rad_per_deg * SIGN(1.,proj_truelat1)
1077      utrue = vgrid * SIN(alpha) + ugrid * COS(alpha)
1078      vtrue = vgrid * COS(alpha) - ugrid * SIN(alpha)
1079    ELSEIF ((proj_code .EQ. PROJ_MERC).OR.(proj_code .EQ.PROJ_LATLON))THEN
1080      utrue = ugrid
1081      vtrue = vgrid
1082    ELSE
1083      PRINT '(A)', 'Unrecognized projection.'
1084      STOP 'GRIDWIND_TO_TRUEWIND'
1085    ENDIF
1086
1087    RETURN
1088  END SUBROUTINE gridwind_to_truewind
1089!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1090
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG