NCEPLIBS-ip 4.0.0
ip_lambert_conf_grid_mod.f90
1module ip_lambert_conf_grid_mod
2 use ip_grid_descriptor_mod
3 use ip_grid_mod
4 use earth_radius_mod
5 use constants_mod
6 implicit none
7
8 private
9 public :: ip_lambert_conf_grid
10
11 type, extends(ip_grid) :: ip_lambert_conf_grid
12 real :: rlat1, rlon1, rlati1, rlati2, orient
13 real :: dxs, dys, h
14 integer :: irot
15 contains
16 procedure :: init_grib1
17 procedure :: init_grib2
18 procedure :: gdswzd => gdswzd_lambert_conf
19 end type ip_lambert_conf_grid
20
21
22 INTEGER :: IROT
23
24 REAL :: AN, DXS, DYS, H
25 REAL :: RERTH
26
27contains
28
29 subroutine init_grib1(self, g1_desc)
30 class(ip_lambert_conf_grid), intent(inout) :: self
31 type(grib1_descriptor), intent(in) :: g1_desc
32
33 real :: dx, dy, hi, hj
34 integer :: iproj, iscan, jscan
35
36 associate(kgds => g1_desc%gds)
37 self%rerth = 6.3712e6
38 self%eccen_squared = 0.0
39
40 self%IM=kgds(2)
41 self%JM=kgds(3)
42
43 self%RLAT1=kgds(4)*1.e-3
44 self%RLON1=kgds(5)*1.e-3
45
46 self%IROT=mod(kgds(6)/8,2)
47 self%ORIENT=kgds(7)*1.e-3
48
49 dx=kgds(8)
50 dy=kgds(9)
51
52 iproj=mod(kgds(10)/128,2)
53 iscan=mod(kgds(11)/128,2)
54 jscan=mod(kgds(11)/64,2)
55
56 self%RLATI1=kgds(12)*1.e-3
57 self%RLATI2=kgds(13)*1.e-3
58 self%H=(-1.)**iproj
59
60 hi=(-1.)**iscan
61 hj=(-1.)**(1-jscan)
62 self%DXS=dx*hi
63 self%DYS=dy*hj
64
65 self%iwrap = 0
66 self%jwrap1 = 0
67 self%jwrap2 = 0
68 self%nscan = mod(kgds(11) / 32, 2)
69 self%nscan_field_pos = self%nscan
70 self%kscan = 0
71 end associate
72
73 end subroutine init_grib1
74
75 subroutine init_grib2(self, g2_desc)
76 class(ip_lambert_conf_grid), intent(inout) :: self
77 type(grib2_descriptor), intent(in) :: g2_desc
78
79 real :: dx, dy, hi, hj
80 integer :: iproj, iscan, jscan
81
82
83 associate(igdtmpl => g2_desc%gdt_tmpl, igdtlen => g2_desc%gdt_len)
84 call earth_radius(igdtmpl, igdtlen, self%rerth, self%eccen_squared)
85
86 self%IM=igdtmpl(8)
87 self%JM=igdtmpl(9)
88
89 self%RLAT1=float(igdtmpl(10))*1.0e-6
90 self%RLON1=float(igdtmpl(11))*1.0e-6
91
92 self%IROT=mod(igdtmpl(12)/8,2)
93 self%ORIENT=float(igdtmpl(14))*1.0e-6
94
95 dx=float(igdtmpl(15))*1.0e-3
96 dy=float(igdtmpl(16))*1.0e-3
97
98 iproj=mod(igdtmpl(17)/128,2)
99 iscan=mod(igdtmpl(18)/128,2)
100 jscan=mod(igdtmpl(18)/64,2)
101
102 self%RLATI1=float(igdtmpl(19))*1.0e-6
103 self%RLATI2=float(igdtmpl(20))*1.0e-6
104
105 self%H=(-1.)**iproj
106 hi=(-1.)**iscan
107 hj=(-1.)**(1-jscan)
108 self%DXS=dx*hi
109 self%DYS=dy*hj
110
111 self%nscan = mod(igdtmpl(18) / 32, 2)
112 self%nscan_field_pos = self%nscan
113 self%iwrap = 0
114 self%jwrap1 = 0
115 self%jwrap2 = 0
116 self%kscan = 0
117 end associate
118 end subroutine init_grib2
119
120 SUBROUTINE gdswzd_lambert_conf(self,IOPT,NPTS,FILL, &
121 XPTS,YPTS,RLON,RLAT,NRET, &
122 CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
123 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
124 !
125 ! SUBPROGRAM: GDSWZD_LAMBERT_CONF GDS WIZARD FOR LAMBERT CONFORMAL CONICAL
126 ! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10
127 !
128 ! ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB 2 GRID DEFINITION
129 ! TEMPLATE (PASSED IN INTEGER FORM AS DECODED BY THE
130 ! NCEP G2 LIBRARY) AND RETURNS ONE OF THE FOLLOWING:
131 ! (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES
132 ! (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES
133 ! WORKS FOR LAMBERT CONFORMAL CONICAL PROJECTIONS.
134 ! IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT
135 ! BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT
136 ! OUTPUT ELEMENTS ARE SET TO FILL VALUES.
137 ! THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO.
138 ! OPTIONALLY, THE VECTOR ROTATIONS, MAP JACOBIANS AND
139 ! GRID BOX AREAS FOR THIS GRID MAY BE RETURNED AS WELL.
140 ! TO COMPUTE THE VECTOR ROTATIONS, THE OPTIONAL ARGUMENTS
141 ! 'SROT' AND 'CROT' MUST BE PRESENT. TO COMPUTE THE MAP
142 ! JACOBIANS, THE OPTIONAL ARGUMENTS 'XLON', 'XLAT',
143 ! 'YLON', 'YLAT' MUST BE PRESENT. TO COMPUTE THE GRID BOX
144 ! AREAS THE OPTIONAL ARGUMENT 'AREA' MUST BE PRESENT.
145 !
146 ! PROGRAM HISTORY LOG:
147 ! 96-04-10 IREDELL
148 ! 96-10-01 IREDELL PROTECTED AGAINST UNRESOLVABLE POINTS
149 ! 97-10-20 IREDELL INCLUDE MAP OPTIONS
150 ! 1999-04-27 GILBERT CORRECTED MINOR ERROR CALCULATING VARIABLE AN
151 ! FOR THE SECANT PROJECTION CASE (RLATI1.NE.RLATI2).
152 ! 2012-08-14 GAYNO FIX PROBLEM WITH SH GRIDS. ENSURE GRID BOX
153 ! AREA ALWAYS POSITIVE.
154 ! 2015-01-21 GAYNO MERGER OF GDSWIZ03 AND GDSWZD03. MAKE
155 ! CROT,SORT,XLON,XLAT,YLON,YLAT AND AREA
156 ! OPTIONAL ARGUMENTS. MAKE PART OF A MODULE.
157 ! MOVE VECTOR ROTATION, MAP JACOBIAN AND GRID
158 ! BOX AREA COMPUTATIONS TO SEPARATE SUBROUTINES.
159 ! 2015-07-13 GAYNO CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAY
160 ! WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAY.
161 ! RENAME ROUTINE AS "GDSWZD_LAMBERT_CONF".
162 ! 2018-07-20 WESLEY ADD THREADS.
163 !
164 ! USAGE: CALL GDSWZD_LAMBERT_CONF(IGDTNUM,IGDTMPL,IGDTLEN,IOPT,NPTS,
165 ! & FILL,XPTS,YPTS,RLON,RLAT,NRET,
166 ! & CROT,SROT,XLON,XLAT,YLON,YLAT,AREA)
167 !
168 ! INPUT ARGUMENT LIST:
169 ! IGDTNUM - INTEGER GRID DEFINITION TEMPLATE NUMBER.
170 ! CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
171 ! NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
172 ! MUST BE "30" FOR LAMBERT CONFORMAL GRIDS.
173 ! IGDTMPL - INTEGER (IGDTLEN) GRID DEFINITION TEMPLATE ARRAY.
174 ! CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT OF THE
175 ! NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE FOR SECTION
176 ! THREE:
177 ! (1): SHAPE OF EARTH, OCTET 15
178 ! (2): SCALE FACTOR OF SPHERICAL EARTH RADIUS,
179 ! OCTET 16
180 ! (3): SCALED VALUE OF RADIUS OF SPHERICAL EARTH,
181 ! OCTETS 17-20
182 ! (4): SCALE FACTOR OF MAJOR AXIS OF ELLIPTICAL EARTH,
183 ! OCTET 21
184 ! (5): SCALED VALUE OF MAJOR AXIS OF ELLIPTICAL EARTH,
185 ! OCTETS 22-25
186 ! (6): SCALE FACTOR OF MINOR AXIS OF ELLIPTICAL EARTH,
187 ! OCTET 26
188 ! (7): SCALED VALUE OF MINOR AXIS OF ELLIPTICAL EARTH,
189 ! OCTETS 27-30
190 ! (8): NUMBER OF POINTS ALONG X-AXIS, OCTS 31-34
191 ! (9): NUMBER OF POINTS ALONG Y-AXIS, OCTS 35-38
192 ! (10): LATITUDE OF FIRST POINT, OCTETS 39-42
193 ! (11): LONGITUDE OF FIRST POINT, OCTETS 43-46
194 ! (12): RESOLUTION OF COMPONENT FLAG, OCTET 47
195 ! (13): LATITUDE WHERE GRID LENGTHS SPECIFIED,
196 ! OCTETS 48-51
197 ! (14): LONGITUDE OF MERIDIAN THAT IS PARALLEL TO
198 ! Y-AXIS, OCTETS 52-55
199 ! (15): X-DIRECTION GRID LENGTH, OCTETS 56-59
200 ! (16): Y-DIRECTION GRID LENGTH, OCTETS 60-63
201 ! (17): PROJECTION CENTER FLAG, OCTET 64
202 ! (18): SCANNING MODE, OCTET 65
203 ! (19): FIRST TANGENT LATITUDE FROM POLE, OCTETS 66-69
204 ! (20): SECOND TANGENT LATITUDE FROM POLE, OCTETS 70-73
205 ! (21): LATITUDE OF SOUTH POLE OF PROJECTION,
206 ! OCTETS 74-77
207 ! (22): LONGITUDE OF SOUTH POLE OF PROJECTION,
208 ! OCTETS 78-81
209 ! IGDTLEN - INTEGER NUMBER OF ELEMENTS (22) OF THE GRID DEFINITION
210 ! TEMPLATE ARRAY. CORRESPONDS TO THE GFLD%IGDTLEN
211 ! COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
212 ! IOPT - INTEGER OPTION FLAG
213 ! (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
214 ! (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
215 ! NPTS - INTEGER MAXIMUM NUMBER OF COORDINATES
216 ! FILL - REAL FILL VALUE TO SET INVALID OUTPUT DATA
217 ! (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
218 ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
219 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
220 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
221 ! (ACCEPTABLE RANGE: -360. TO 360.)
222 ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
223 ! (ACCEPTABLE RANGE: -90. TO 90.)
224 !
225 ! OUTPUT ARGUMENT LIST:
226 ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0
227 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0
228 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0
229 ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0
230 ! NRET - INTEGER NUMBER OF VALID POINTS COMPUTED
231 ! CROT - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES
232 ! SROT - REAL, OPTIONAL (NPTS) CLOCKWISE VECTOR ROTATION SINES
233 ! (UGRID=CROT*UEARTH-SROT*VEARTH;
234 ! VGRID=SROT*UEARTH+CROT*VEARTH)
235 ! XLON - REAL, OPTIONAL (NPTS) DX/DLON IN 1/DEGREES
236 ! XLAT - REAL, OPTIONAL (NPTS) DX/DLAT IN 1/DEGREES
237 ! YLON - REAL, OPTIONAL (NPTS) DY/DLON IN 1/DEGREES
238 ! YLAT - REAL, OPTIONAL (NPTS) DY/DLAT IN 1/DEGREES
239 ! AREA - REAL, OPTIONAL (NPTS) AREA WEIGHTS IN M**2
240 ! (PROPORTIONAL TO THE SQUARE OF THE MAP FACTOR)
241 !
242 ! ATTRIBUTES:
243 ! LANGUAGE: FORTRAN 90
244 !
245 !$$$
246 IMPLICIT NONE
247 !
248 class(ip_lambert_conf_grid), intent(in) :: self
249 INTEGER, INTENT(IN ) :: IOPT, NPTS
250 INTEGER, INTENT( OUT) :: NRET
251 !
252 REAL, INTENT(IN ) :: FILL
253 REAL, INTENT(INOUT) :: RLON(NPTS),RLAT(NPTS)
254 REAL, INTENT(INOUT) :: XPTS(NPTS),YPTS(NPTS)
255 REAL, OPTIONAL, INTENT( OUT) :: CROT(NPTS),SROT(NPTS)
256 REAL, OPTIONAL, INTENT( OUT) :: XLON(NPTS),XLAT(NPTS)
257 REAL, OPTIONAL, INTENT( OUT) :: YLON(NPTS),YLAT(NPTS),AREA(NPTS)
258 !
259 INTEGER :: IM, JM, N
260 !
261 LOGICAL :: LROT, LMAP, LAREA
262 !
263 REAL :: ANTR, DI, DJ
264 REAL :: DLON1
265 REAL :: DE, DE2, DR2
266 REAL :: ORIENT, RLAT1, RLON1
267 REAL :: RLATI1, RLATI2
268 REAL :: XMAX, XMIN, YMAX, YMIN, XP, YP
269 REAL :: DLON, DR
270 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271 IF(PRESENT(crot)) crot=fill
272 IF(PRESENT(srot)) srot=fill
273 IF(PRESENT(xlon)) xlon=fill
274 IF(PRESENT(xlat)) xlat=fill
275 IF(PRESENT(ylon)) ylon=fill
276 IF(PRESENT(ylat)) ylat=fill
277 IF(PRESENT(area)) area=fill
278
279 im=self%im
280 jm=self%jm
281
282 rlat1=self%rlat1
283 rlon1=self%rlon1
284
285 irot=self%irot
286 orient=self%orient
287
288 rlati1=self%rlati1
289 rlati2=self%rlati2
290
291 h=self%h
292 dxs=self%dxs
293 dys=self%dys
294
295 rerth = self%rerth
296
297 IF(rlati1.EQ.rlati2) THEN
298 an=sin(rlati1/dpr)
299 ELSE
300 an=log(cos(rlati1/dpr)/cos(rlati2/dpr))/ &
301 log(tan((90-rlati1)/2/dpr)/tan((90-rlati2)/2/dpr))
302 ENDIF
303 de=rerth*cos(rlati1/dpr)*tan((rlati1+90)/2/dpr)**an/an
304 IF(h*rlat1.EQ.90) THEN
305 xp=1
306 yp=1
307 ELSE
308 dr=de/tan((rlat1+90)/2/dpr)**an
309 dlon1=mod(rlon1-orient+180+3600,360.)-180
310 xp=1-sin(an*dlon1/dpr)*dr/dxs
311 yp=1+cos(an*dlon1/dpr)*dr/dys
312 ENDIF
313 antr=1/(2*an)
314 de2=de**2
315 xmin=0
316 xmax=im+1
317 ymin=0
318 ymax=jm+1
319 nret=0
320 IF(PRESENT(crot).AND.PRESENT(srot))THEN
321 lrot=.true.
322 ELSE
323 lrot=.false.
324 ENDIF
325 IF(PRESENT(xlon).AND.PRESENT(xlat).AND.PRESENT(ylon).AND.PRESENT(ylat))THEN
326 lmap=.true.
327 ELSE
328 lmap=.false.
329 ENDIF
330 IF(PRESENT(area))THEN
331 larea=.true.
332 ELSE
333 larea=.false.
334 ENDIF
335 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
336 ! TRANSLATE GRID COORDINATES TO EARTH COORDINATES
337 IF(iopt.EQ.0.OR.iopt.EQ.1) THEN
338 !$OMP PARALLEL DO PRIVATE(N,DI,DJ,DR2,DR,DLON) REDUCTION(+:NRET) SCHEDULE(STATIC)
339 DO n=1,npts
340 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
341 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
342 di=h*(xpts(n)-xp)*dxs
343 dj=h*(ypts(n)-yp)*dys
344 dr2=di**2+dj**2
345 dr=sqrt(dr2)
346 IF(dr2.LT.de2*1.e-6) THEN
347 rlon(n)=0.
348 rlat(n)=h*90.
349 ELSE
350 rlon(n)=mod(orient+1./an*dpr*atan2(di,-dj)+3600,360.)
351 rlat(n)=(2*dpr*atan((de2/dr2)**antr)-90)
352 ENDIF
353 nret=nret+1
354 dlon=mod(rlon(n)-orient+180+3600,360.)-180
355 IF(lrot) CALL lambert_conf_vect_rot(dlon,crot(n),srot(n))
356 IF(lmap) CALL lambert_conf_map_jacob(rlat(n),fill, dlon, dr, &
357 xlon(n),xlat(n),ylon(n),ylat(n))
358 IF(larea) CALL lambert_conf_grid_area(rlat(n),fill,dr,area(n))
359 ELSE
360 rlon(n)=fill
361 rlat(n)=fill
362 ENDIF
363 ENDDO
364 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
365 ! TRANSLATE EARTH COORDINATES TO GRID COORDINATES
366 ELSEIF(iopt.EQ.-1) THEN
367 !$OMP PARALLEL DO PRIVATE(N,DR,DLON) REDUCTION(+:NRET) SCHEDULE(STATIC)
368 DO n=1,npts
369 IF(abs(rlon(n)).LE.360.AND.abs(rlat(n)).LE.90.AND. &
370 h*rlat(n).NE.-90) THEN
371 dr=h*de*tan((90-rlat(n))/2/dpr)**an
372 dlon=mod(rlon(n)-orient+180+3600,360.)-180
373 xpts(n)=xp+h*sin(an*dlon/dpr)*dr/dxs
374 ypts(n)=yp-h*cos(an*dlon/dpr)*dr/dys
375 IF(xpts(n).GE.xmin.AND.xpts(n).LE.xmax.AND. &
376 ypts(n).GE.ymin.AND.ypts(n).LE.ymax) THEN
377 nret=nret+1
378 IF(lrot) CALL lambert_conf_vect_rot(dlon,crot(n),srot(n))
379 IF(lmap) CALL lambert_conf_map_jacob(rlat(n),fill,dlon,dr, &
380 xlon(n),xlat(n),ylon(n),ylat(n))
381 IF(larea) CALL lambert_conf_grid_area(rlat(n),fill,dr,area(n))
382 ELSE
383 xpts(n)=fill
384 ypts(n)=fill
385 ENDIF
386 ELSE
387 xpts(n)=fill
388 ypts(n)=fill
389 ENDIF
390 ENDDO
391 !$OMP END PARALLEL DO
392 ENDIF
393 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
394 END SUBROUTINE gdswzd_lambert_conf
395 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396 SUBROUTINE lambert_conf_error(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
397 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
398 !
399 ! SUBPROGRAM: LAMBERT_CONF_ERROR ERROR HANDLER
400 ! PRGMMR: GAYNO ORG: W/NMC23 DATE: 2015-07-13
401 !
402 ! ABSTRACT: UPON AN ERROR, THIS SUBPROGRAM ASSIGNS
403 ! A "FILL" VALUE TO THE OUTPUT FIELDS.
404
405 ! PROGRAM HISTORY LOG:
406 ! 2015-07-13 GAYNO INITIAL VERSION
407 !
408 ! USAGE: CALL LAMBERT_CONF_ERROR(IOPT,FILL,RLAT,RLON,XPTS,YPTS,NPTS)
409 !
410 ! INPUT ARGUMENT LIST:
411 ! IOPT - INTEGER OPTION FLAG
412 ! (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS)
413 ! (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS)
414 ! NPTS - INTEGER MAXIMUM NUMBER OF COORDINATES
415 ! FILL - REAL FILL VALUE TO SET INVALID OUTPUT DATA
416 ! (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.)
417 ! OUTPUT ARGUMENT LIST:
418 ! RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0
419 ! RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0
420 ! XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0
421 ! YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0
422 !
423 ! ATTRIBUTES:
424 ! LANGUAGE: FORTRAN 90
425 !
426 !$$$
427 IMPLICIT NONE
428 !
429 INTEGER, INTENT(IN ) :: IOPT, NPTS
430 !
431 REAL, INTENT(IN ) :: FILL
432 REAL, INTENT( OUT) :: RLAT(NPTS),RLON(NPTS)
433 REAL, INTENT( OUT) :: XPTS(NPTS),YPTS(NPTS)
434 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
435 IF(iopt>=0) THEN
436 rlon=fill
437 rlat=fill
438 ENDIF
439 IF(iopt<=0) THEN
440 xpts=fill
441 ypts=fill
442 ENDIF
443 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
444 END SUBROUTINE lambert_conf_error
445 !
446 SUBROUTINE lambert_conf_vect_rot(DLON,CROT,SROT)
447 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
448 !
449 ! SUBPROGRAM: LAMBERT_CONF_VECT_ROT VECTOR ROTATION FIELDS FOR
450 ! LAMBERT CONFORMAL CONICAL
451 !
452 ! PRGMMR: GAYNO ORG: W/NMC23 DATE: 2015-01-21
453 !
454 ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE VECTOR ROTATION SINES AND
455 ! COSINES FOR A LAMBERT CONFORMAL CONICAL GRID
456 !
457 ! PROGRAM HISTORY LOG:
458 ! 2015-01-21 GAYNO INITIAL VERSION
459 ! 2015-09-17 GAYNO RENAME AS "LAMBERT_CONF_VECT_ROT"
460 ! 2018-07-20 WESLEY PASS IN DLON FOR THREADING.
461 !
462 ! USAGE: CALL LAMBERT_CONF_VECT_ROT(DLON,CROT,SROT)
463 !
464 ! INPUT ARGUMENT LIST:
465 ! DLON - DISTANCE FROM ORIENTATION LONGITUDE (REAL)
466 !
467 ! OUTPUT ARGUMENT LIST:
468 ! CROT - CLOCKWISE VECTOR ROTATION COSINES (REAL)
469 ! SROT - CLOCKWISE VECTOR ROTATION SINES (REAL)
470 ! (UGRID=CROT*UEARTH-SROT*VEARTH;
471 ! VGRID=SROT*UEARTH+CROT*VEARTH)
472 !
473 ! ATTRIBUTES:
474 ! LANGUAGE: FORTRAN 90
475 !
476 !$$$
477 IMPLICIT NONE
478 REAL, INTENT( IN) :: DLON
479 REAL, INTENT( OUT) :: CROT, SROT
480
481 IF(irot.EQ.1) THEN
482 crot=cos(an*dlon/dpr)
483 srot=sin(an*dlon/dpr)
484 ELSE
485 crot=1.
486 srot=0.
487 ENDIF
488
489 END SUBROUTINE lambert_conf_vect_rot
490 !
491 SUBROUTINE lambert_conf_map_jacob(RLAT,FILL,DLON,DR,XLON,XLAT,YLON,YLAT)
492 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
493 !
494 ! SUBPROGRAM: LAMBERT_CONF_MAP_JACOB MAP JACOBIANS FOR
495 ! LAMBERT CONFORMAL CONICAL
496 !
497 ! PRGMMR: GAYNO ORG: W/NMC23 DATE: 2015-01-21
498 !
499 ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE MAP JACOBIANS FOR
500 ! A LAMBERT CONFORMAL CONICAL GRID.
501 !
502 ! PROGRAM HISTORY LOG:
503 ! 2015-01-21 GAYNO INITIAL VERSION
504 ! 2015-09-17 GAYNO RENAME AS "LAMBERT_CONF_MAP_JACOB"
505 ! 2018-07-20 WESLEY PASS DLON AND DR FOR THREADING.
506 !
507 ! USAGE: CALL LAMBERT_CONF_MAP_JACOB(RLAT,FILL,DLON,DR,XLON,XLAT,YLON,YLAT)
508 !
509 ! INPUT ARGUMENT LIST:
510 ! RLAT - GRID POINT LATITUDE IN DEGREES (REAL)
511 ! FILL - FILL VALUE FOR UNDEFINED POINTS (REAL)
512 ! DLON - DISTANCE FROM ORIENTATION LONGITUDE (REAL)
513 ! DR - DISTANCE FROM POLE POINT (REAL)
514 !
515 ! OUTPUT ARGUMENT LIST:
516 ! XLON - DX/DLON IN 1/DEGREES (REAL)
517 ! XLAT - DX/DLAT IN 1/DEGREES (REAL)
518 ! YLON - DY/DLON IN 1/DEGREES (REAL)
519 ! YLAT - DY/DLAT IN 1/DEGREES (REAL)
520 !
521 ! ATTRIBUTES:
522 ! LANGUAGE: FORTRAN 90
523 !
524 !$$$
525
526 IMPLICIT NONE
527
528 REAL, INTENT(IN ) :: RLAT, FILL, DLON, DR
529 REAL, INTENT( OUT) :: XLON, XLAT, YLON, YLAT
530
531 REAL :: CLAT
532
533 clat=cos(rlat/dpr)
534 IF(clat.LE.0.OR.dr.LE.0) THEN
535 xlon=fill
536 xlat=fill
537 ylon=fill
538 ylat=fill
539 ELSE
540 xlon=h*cos(an*dlon/dpr)*an/dpr*dr/dxs
541 xlat=-h*sin(an*dlon/dpr)*an/dpr*dr/dxs/clat
542 ylon=h*sin(an*dlon/dpr)*an/dpr*dr/dys
543 ylat=h*cos(an*dlon/dpr)*an/dpr*dr/dys/clat
544 ENDIF
545
546 END SUBROUTINE lambert_conf_map_jacob
547 !
548 SUBROUTINE lambert_conf_grid_area(RLAT,FILL,DR,AREA)
549 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
550 !
551 ! SUBPROGRAM: LAMBERT_CONF_GRID_AREA GRID BOX AREA FOR
552 ! LAMBERT CONFORMAL CONICAL
553 !
554 ! PRGMMR: GAYNO ORG: W/NMC23 DATE: 2015-01-21
555 !
556 ! ABSTRACT: THIS SUBPROGRAM COMPUTES THE GRID BOX AREA FOR
557 ! A LAMBERT CONFORMAL CONICAL GRID.
558 !
559 ! PROGRAM HISTORY LOG:
560 ! 2015-01-21 GAYNO INITIAL VERSION
561 ! 2015-09-17 GAYNO RENAME AS "LAMBERT_CONF_GRID_AREA"
562 ! 2018-07-20 WESLEY PASS IN DR FOR THREADING.
563 !
564 ! USAGE: CALL LAMBERT_CONF_GRID_AREA(RLAT,FILL,DR,AREA)
565 !
566 ! INPUT ARGUMENT LIST:
567 ! RLAT - LATITUDE OF GRID POINT IN DEGREES (REAL)
568 ! FILL - FILL VALUE FOR UNDEFINED POINTS (REAL)
569 ! DR - DISTANCE FROM POLE POINT (REAL)
570 !
571 ! OUTPUT ARGUMENT LIST:
572 ! AREA - AREA WEIGHTS IN M**2 (REAL)
573 !
574 ! ATTRIBUTES:
575 ! LANGUAGE: FORTRAN 90
576 !
577 !$$$
578
579 IMPLICIT NONE
580
581 REAL, INTENT(IN ) :: RLAT
582 REAL, INTENT(IN ) :: FILL
583 REAL, INTENT(IN ) :: DR
584 REAL, INTENT( OUT) :: AREA
585
586 REAL :: CLAT
587
588 clat=cos(rlat/dpr)
589 IF(clat.LE.0.OR.dr.LE.0) THEN
590 area=fill
591 ELSE
592 area=rerth**2*clat**2*abs(dxs)*abs(dys)/(an*dr)**2
593 ENDIF
594
595 END SUBROUTINE lambert_conf_grid_area
596
597end module ip_lambert_conf_grid_mod
598
constants_mod
Module containing common constants.
Definition: constants_mod.f90:9
ip_grid_descriptor_mod
Uses derived type grid descriptor objects to abstract away the raw Grib-1 and Grib-2 grid definitions...
Definition: ip_grid_descriptor_mod.f90:15
ip_grid_mod
Abstract ip_grid type.
Definition: ip_grid_mod.f90:8
ip_grid_mod::ip_grid
Abstract grid that holds fields and methods common to all grids.
Definition: ip_grid_mod.f90:45
ip_lambert_conf_grid_mod::ip_lambert_conf_grid
Definition: ip_lambert_conf_grid_mod.f90:11