Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: regrid_a2a_mod.F90
7 : !
8 : ! !DESCRIPTION: Module HCO\_REGRID\_A2A\_MOD uses an algorithm adapted from
9 : ! MAP\_A2A code to regrid from one horizontal grid to another.
10 : !\\
11 : !\\
12 : ! !INTERFACE:
13 : !
14 : MODULE HCO_Regrid_A2A_Mod
15 : !
16 : ! !USES:
17 : !
18 : USE HCO_PRECISION_MOD ! For GEOS-Chem Precision (fp)
19 :
20 : IMPLICIT NONE
21 : PRIVATE
22 : !
23 : ! !PUBLIC MEMBER FUNCTIONS:
24 : !
25 : PUBLIC :: Map_A2A
26 : PUBLIC :: Init_Map_A2A
27 : PUBLIC :: Cleanup_Map_A2A
28 :
29 : ! Map_A2A overloads these routines
30 : INTERFACE Map_A2A
31 : MODULE PROCEDURE Map_A2A_R8R8
32 : MODULE PROCEDURE Map_A2A_R4R8
33 : MODULE PROCEDURE Map_A2A_R4R4
34 : MODULE PROCEDURE Map_A2A_R8R4
35 : END INTERFACE Map_A2A
36 : !
37 : ! !PRIVATE MEMBER FUNCTIONS:
38 : !
39 : PRIVATE :: Map_A2A_R8R8
40 : PRIVATE :: Map_A2A_R4R4
41 : PRIVATE :: Map_A2A_R4R8
42 : PRIVATE :: Map_A2A_R8R4
43 : PRIVATE :: Ymap_R8R8
44 : PRIVATE :: Ymap_R4R8
45 : PRIVATE :: Ymap_R4R4
46 : PRIVATE :: Ymap_R8R4
47 : PRIVATE :: Xmap_R8R8
48 : PRIVATE :: Xmap_R4R4
49 : PRIVATE :: Xmap_R4R8
50 : PRIVATE :: Xmap_R8R4
51 : !
52 : ! !REVISION HISTORY:
53 : ! 13 Mar 2012 - M. Cooper - Initial version
54 : ! See https://github.com/geoschem/hemco for complete history
55 : !EOP
56 : !------------------------------------------------------------------------------
57 : !BOC
58 : !
59 : ! !PRIVATE TYPES:
60 : !
61 :
62 : !---------------------------------------------------------------------------
63 : ! These are now kept locally, to "shadow" variables from other parts of
64 : ! GEOS-Chem. This avoids depending on GEOS-Chem code within the core
65 : ! HEMCO modules. (bmy, 7/14/14)
66 : !---------------------------------------------------------------------------
67 : CHARACTER(LEN=255) :: NC_DIR ! Directory w/ netCDF files
68 : INTEGER :: OUTNX ! # of longitudes (x-dimension) in grid
69 : INTEGER :: OUTNY ! # of latitudes (y-dimension) in grid
70 : REAL(fp), ALLOCATABLE :: OUTLON (: ) ! Longitude on output grid
71 : REAL(fp), ALLOCATABLE :: OUTSIN (: ) ! Sines of latitudes on output grid
72 : REAL(fp), ALLOCATABLE :: OUTAREA(:,:) ! Surface areas on output grid
73 : !
74 : ! !DEFINED PARAMETERS:
75 : !
76 : !---------------------------------------------------------------------------
77 : ! Physical constants taken from the GEOS-Chem "physconstants.F90" module,
78 : ! which uses values from NIST 2014. (ewl, bmy, 03 Mar 2022)
79 : !---------------------------------------------------------------------------
80 : REAL(fp), PARAMETER :: PI = 3.14159265358979323_fp ! Pi
81 : REAL(fp), PARAMETER :: Re = 6.3710072e+6_fp ! Earth radius [m]
82 :
83 : !---------------------------------------------------------------------------
84 : ! Tiny numbers for single and double precision. These are being used for
85 : ! skipping missing values. miss_r4 and miss_r8 are the default missing values
86 : ! for single and double precision, respectively. (ckeller, 4/8/2017)
87 : !---------------------------------------------------------------------------
88 : REAL*4, PARAMETER :: tiny_r4 = 1.0e-30 !1.0e-20
89 : REAL*4, PARAMETER :: miss_r4 = 0.0e0
90 : REAL*8, PARAMETER :: tiny_r8 = 1.0d-40
91 : REAL*8, PARAMETER :: miss_r8 = 0.0d0
92 :
93 : CONTAINS
94 : !EOC
95 : !------------------------------------------------------------------------------
96 : ! Harmonized Emissions Component (HEMCO) !
97 : !------------------------------------------------------------------------------
98 : !BOP
99 : !
100 : ! !IROUTINE: Map_A2A_r8r8
101 : !
102 : ! !DESCRIPTION: Subroutine MAP\_A2A\_R8R8 is a horizontal arbitrary grid to
103 : ! arbitrary grid conservative high-order mapping regridding routine by S-J
104 : ! Lin. Both the input data and output data have REAL(fp) precision.
105 : !\\
106 : !\\
107 : ! !INTERFACE:
108 : !
109 0 : SUBROUTINE Map_A2A_r8r8( im, jm, lon1, sin1, q1, &
110 0 : in, jn, lon2, sin2, q2, ig, iv, missval)
111 : !
112 : ! !INPUT PARAMETERS:
113 : !
114 : ! Longitude and Latitude dimensions of INPUT grid
115 : INTEGER, INTENT(IN) :: im, jm
116 :
117 : ! Longitude and Latitude dimensions of OUTPUT grid
118 : INTEGER, INTENT(IN) :: in, jn
119 :
120 : ! IG=0: pole to pole;
121 : ! IG=1 J=1 is half-dy north of south pole
122 : INTEGER, INTENT(IN) :: ig
123 :
124 : ! IV=0: Regrid scalar quantity
125 : ! IV=1: Regrid vector quantity
126 : INTEGER, INTENT(IN) :: iv
127 :
128 : ! Longitude edges (degrees) of INPUT and OUTPUT grids
129 : REAL*8, INTENT(IN) :: lon1(im+1), lon2(in+1)
130 :
131 : ! Sine of Latitude Edges (radians) of INPUT and OUTPUT grids
132 : REAL*8, INTENT(IN) :: sin1(jm+1), sin2(jn+1)
133 :
134 : ! Quantity on INPUT grid
135 : REAL*8, INTENT(IN) :: q1(im,jm)
136 :
137 : !
138 : ! !OUTPUT PARAMETERS:
139 : !
140 : ! Regridded quantity on OUTPUT grid
141 : REAL*8, INTENT(OUT) :: q2(in,jn)
142 : !
143 : ! !OPTIONAL ARGUMENTS
144 : !
145 : REAL*8, INTENT(IN), OPTIONAL :: missval
146 : !
147 : ! !REMARKS:
148 : ! This routine is overloaded by the MAP_A2A interface.
149 : !
150 : ! !REVISION HISTORY:
151 : ! (1) Original subroutine by S-J Lin. Converted to F90 freeform format
152 : ! and inserted into "Geos3RegridModule" by Bob Yantosca (9/21/00)
153 : ! See https://github.com/geoschem/hemco for complete history
154 : !EOP
155 : !------------------------------------------------------------------------------
156 : !BOC
157 : !
158 : ! !LOCAL VARIABLES:
159 : !
160 : INTEGER :: i,j,k
161 0 : REAL*8 :: qtmp(in,jm)
162 :
163 : ! Init
164 0 : IF ( PRESENT(missval) ) THEN
165 0 : qtmp = missval
166 0 : q2 = missval
167 : ELSE
168 0 : qtmp = miss_r8
169 0 : q2 = miss_r8
170 : ENDIF
171 :
172 : !===================================================================
173 : ! E-W regridding
174 : !===================================================================
175 0 : IF ( im == in .and. &
176 0 : lon1(1) == lon2(1) .and. &
177 : lon1(im+1) == lon2(in+1) ) THEN
178 :
179 : ! Don't call XMAP if both grids have the same # of longitudes
180 : ! but save the input data in the QTMP array
181 : !$OMP PARALLEL DO &
182 : !$OMP DEFAULT( SHARED ) &
183 : !$OMP PRIVATE( I, J )
184 0 : DO j=1,jm-ig
185 0 : DO i=1,im
186 0 : qtmp(i,j+ig) = q1(i,j+ig)
187 : ENDDO
188 : ENDDO
189 : !$OMP END PARALLEL DO
190 :
191 : ELSE
192 :
193 : ! Otherwise, call XMAP to regrid in the E-W direction
194 0 : CALL xmap_r8r8(im, jm-ig, lon1, q1(1,1+ig),in, lon2, qtmp(1,1+ig), &
195 0 : missval=missval )
196 :
197 : ENDIF
198 :
199 : !===================================================================
200 : ! N-S regridding
201 : !===================================================================
202 0 : IF ( jm == jn .and. &
203 0 : sin1(1) == sin2(1) .and. &
204 : sin1(jm+1) == sin2(jn+1) ) THEN
205 :
206 : ! Don't call XMAP if both grids have the same # of longitudes,
207 : ! but assign the value of QTMP to the output Q2 array
208 : !$OMP PARALLEL DO &
209 : !$OMP DEFAULT( SHARED ) &
210 : !$OMP PRIVATE( I, J )
211 0 : DO j=1,jm-ig
212 0 : DO i=1,in
213 0 : q2(i,j+ig) = qtmp(i,j+ig)
214 : ENDDO
215 : ENDDO
216 : !$OMP END PARALLEL DO
217 :
218 : ELSE
219 :
220 : ! Otherwise, call YMAP to regrid in the N-S direction
221 0 : CALL ymap_r8r8(in, jm, sin1, qtmp(1,1+ig), jn, sin2, q2(1,1+ig), ig, iv, &
222 0 : missval=missval )
223 :
224 : ENDIF
225 :
226 0 : END SUBROUTINE Map_A2A_r8r8
227 : !EOC
228 : !------------------------------------------------------------------------------
229 : ! Harmonized Emissions Component (HEMCO) !
230 : !------------------------------------------------------------------------------
231 : !BOP
232 : !
233 : ! !IROUTINE: Map_A2A_r4r4
234 : !
235 : ! !DESCRIPTION: Subroutine MAP\_A2A\_R4R4 is a horizontal arbitrary grid
236 : ! to arbitrary grid conservative high-order mapping regridding routine
237 : ! by S-J Lin. Both the input and output data have REAL*4 precision.
238 : !\\
239 : !\\
240 : ! !INTERFACE:
241 : !
242 0 : SUBROUTINE Map_A2A_r4r4( im, jm, lon1, sin1, q1, &
243 0 : in, jn, lon2, sin2, q2, ig, iv, missval)
244 : !
245 : ! !INPUT PARAMETERS:
246 : !
247 : ! Longitude and Latitude dimensions of INPUT grid
248 : INTEGER, INTENT(IN) :: im, jm
249 :
250 : ! Longitude and Latitude dimensions of OUTPUT grid
251 : INTEGER, INTENT(IN) :: in, jn
252 :
253 : ! IG=0: pole to pole;
254 : ! IG=1 J=1 is half-dy north of south pole
255 : INTEGER, INTENT(IN) :: ig
256 :
257 : ! IV=0: Regrid scalar quantity
258 : ! IV=1: Regrid vector quantity
259 : INTEGER, INTENT(IN) :: iv
260 :
261 : ! Longitude edges (degrees) of INPUT and OUTPUT grids
262 : REAL*4, INTENT(IN) :: lon1(im+1), lon2(in+1)
263 :
264 : ! Sine of Latitude Edges (radians) of INPUT and OUTPUT grids
265 : REAL*4, INTENT(IN) :: sin1(jm+1), sin2(jn+1)
266 :
267 : ! Quantity on INPUT grid
268 : REAL*4, INTENT(IN) :: q1(im,jm)
269 : !
270 : ! !OUTPUT PARAMETERS:
271 : !
272 : ! Regridded quantity on OUTPUT grid
273 : REAL*4, INTENT(OUT) :: q2(in,jn)
274 : !
275 : ! !OPTIONAL ARGUMENTS
276 : !
277 : REAL*4, INTENT(IN), OPTIONAL :: missval
278 : !
279 : ! !REMARKS:
280 : ! This routine is overloaded by the MAP_A2A interface.
281 : !
282 : ! !REVISION HISTORY:
283 : ! (1) Original subroutine by S-J Lin. Converted to F90 freeform format
284 : ! and inserted into "Geos3RegridModule" by Bob Yantosca (9/21/00)
285 : ! See https://github.com/geoschem/hemco for complete history
286 : !EOP
287 : !------------------------------------------------------------------------------
288 : !BOC
289 : !
290 : ! !LOCAL VARIABLES:
291 : !
292 : INTEGER :: i,j,k
293 0 : REAL*4 :: qtmp(in,jm)
294 :
295 : ! Init
296 0 : IF ( PRESENT(missval) ) THEN
297 0 : qtmp = missval
298 0 : q2 = missval
299 : ELSE
300 0 : qtmp = miss_r4
301 0 : q2 = miss_r4
302 : ENDIF
303 :
304 : !===================================================================
305 : ! E-W regridding
306 : !===================================================================
307 0 : IF ( im == in .and. &
308 0 : lon1(1) == lon2(1) .and. &
309 : lon1(im+1) == lon2(in+1) ) THEN
310 :
311 : ! Don't call XMAP if both grids have the same # of longitudes
312 : ! but save the input data in the QTMP array
313 : !$OMP PARALLEL DO &
314 : !$OMP DEFAULT( SHARED ) &
315 : !$OMP PRIVATE( I, J )
316 0 : DO j=1,jm-ig
317 0 : DO i=1,im
318 0 : qtmp(i,j+ig) = q1(i,j+ig)
319 : ENDDO
320 : ENDDO
321 : !$OMP END PARALLEL DO
322 :
323 : ELSE
324 :
325 : ! Otherwise, call XMAP to regrid in the E-W direction
326 0 : CALL xmap_r4r4(im, jm-ig, lon1, q1(1,1+ig),in, lon2, qtmp(1,1+ig), &
327 0 : missval=missval )
328 :
329 : ENDIF
330 :
331 : !===================================================================
332 : ! N-S regridding
333 : !===================================================================
334 0 : IF ( jm == jn .and. &
335 0 : sin1(1) == sin2(1) .and. &
336 : sin1(jm+1) == sin2(jn+1) ) THEN
337 :
338 : ! Don't call XMAP if both grids have the same # of longitudes,
339 : ! but assign the value of QTMP to the output Q2 array
340 : !$OMP PARALLEL DO &
341 : !$OMP DEFAULT( SHARED ) &
342 : !$OMP PRIVATE( I, J )
343 0 : DO j=1,jm-ig
344 0 : DO i=1,in
345 0 : q2(i,j+ig) = qtmp(i,j+ig)
346 : ENDDO
347 : ENDDO
348 : !$OMP END PARALLEL DO
349 :
350 : ELSE
351 :
352 : ! Otherwise, call YMAP to regrid in the N-S direction
353 0 : CALL ymap_r4r4(in, jm, sin1, qtmp(1,1+ig), jn, sin2, q2(1,1+ig), ig, iv, &
354 0 : missval=missval)
355 :
356 : ENDIF
357 :
358 0 : END SUBROUTINE Map_A2A_r4r4
359 : !EOC
360 : !------------------------------------------------------------------------------
361 : ! Harmonized Emissions Component (HEMCO) !
362 : !------------------------------------------------------------------------------
363 : !BOP
364 : !
365 : ! !IROUTINE: Map_A2A_r4r8
366 : !
367 : ! !DESCRIPTION: Subroutine MAP\_A2A\_R4R8 is a horizontal arbitrary grid to
368 : ! arbitrary grid conservative high-order mapping regridding routine by
369 : ! S-J Lin. The input data has REAL*4 precision, but the output argument
370 : ! has REAL(fp) precision.
371 : !\\
372 : !\\
373 : ! !INTERFACE:
374 : !
375 0 : SUBROUTINE Map_A2A_r4r8( im, jm, lon1, sin1, q1, &
376 0 : in, jn, lon2, sin2, q2, ig, iv, missval)
377 : !
378 : ! !INPUT PARAMETERS:
379 : !
380 : ! Longitude and Latitude dimensions of INPUT grid
381 : INTEGER, INTENT(IN) :: im, jm
382 :
383 : ! Longitude and Latitude dimensions of OUTPUT grid
384 : INTEGER, INTENT(IN) :: in, jn
385 :
386 : ! IG=0: pole to pole;
387 : ! IG=1 J=1 is half-dy north of south pole
388 : INTEGER, INTENT(IN) :: ig
389 :
390 : ! IV=0: Regrid scalar quantity
391 : ! IV=1: Regrid vector quantity
392 : INTEGER, INTENT(IN) :: iv
393 :
394 : ! Longitude edges (degrees) of INPUT and OUTPUT grids
395 : REAL*4, INTENT(IN) :: lon1(im+1), lon2(in+1)
396 :
397 : ! Sine of Latitude Edges (radians) of INPUT and OUTPUT grids
398 : REAL*4, INTENT(IN) :: sin1(jm+1), sin2(jn+1)
399 :
400 : ! Quantity on INPUT grid
401 : REAL*4, INTENT(IN) :: q1(im,jm)
402 : !
403 : ! !OUTPUT PARAMETERS:
404 : !
405 : ! Regridded quantity on OUTPUT grid
406 : REAL*8, INTENT(OUT) :: q2(in,jn)
407 : !
408 : ! !OPTIONAL ARGUMENTS
409 : !
410 : REAL*4, INTENT(IN), OPTIONAL :: missval
411 : !
412 : ! !REMARKS:
413 : ! This routine is overloaded by the MAP_A2A interface.
414 : !
415 : ! !REVISION HISTORY:
416 : ! (1) Original subroutine by S-J Lin. Converted to F90 freeform format
417 : ! and inserted into "Geos3RegridModule" by Bob Yantosca (9/21/00)
418 : ! See https://github.com/geoschem/hemco for complete history
419 : !EOP
420 : !------------------------------------------------------------------------------
421 : !BOC
422 : !
423 : ! !LOCAL VARIABLES:
424 : !
425 : INTEGER :: i,j,k
426 0 : REAL*8 :: qtmp(in,jm)
427 :
428 : ! Init
429 0 : IF ( PRESENT(missval) ) THEN
430 0 : qtmp = real(missval,kind=f8)
431 0 : q2 = real(missval,kind=f8)
432 : ELSE
433 0 : qtmp = miss_r8
434 0 : q2 = miss_r8
435 : ENDIF
436 :
437 : !===================================================================
438 : ! E-W regridding
439 : !===================================================================
440 0 : IF ( im == in .and. &
441 0 : lon1(1) == lon2(1) .and. &
442 : lon1(im+1) == lon2(in+1) ) THEN
443 :
444 : ! Don't call XMAP if both grids have the same # of longitudes
445 : ! but save the input data in the QTMP array
446 : !$OMP PARALLEL DO &
447 : !$OMP DEFAULT( SHARED ) &
448 : !$OMP PRIVATE( I, J )
449 0 : DO j=1,jm-ig
450 0 : DO i=1,im
451 0 : qtmp(i,j+ig) = q1(i,j+ig)
452 : ENDDO
453 : ENDDO
454 : !$OMP END PARALLEL DO
455 :
456 : ELSE
457 :
458 : ! Otherwise, call XMAP to regrid in the E-W direction
459 0 : CALL xmap_r4r8(im, jm-ig, lon1, q1(1,1+ig),in, lon2, qtmp(1,1+ig), &
460 0 : missval=missval )
461 :
462 : ENDIF
463 :
464 : !===================================================================
465 : ! N-S regridding
466 : !===================================================================
467 0 : IF ( jm == jn .and. &
468 0 : sin1(1) == sin2(1) .and. &
469 : sin1(jm+1) == sin2(jn+1) ) THEN
470 :
471 : ! Don't call XMAP if both grids have the same # of longitudes,
472 : ! but assign the value of QTMP to the output Q2 array
473 : !$OMP PARALLEL DO &
474 : !$OMP DEFAULT( SHARED ) &
475 : !$OMP PRIVATE( I, J )
476 0 : DO j=1,jm-ig
477 0 : DO i=1,in
478 0 : q2(i,j+ig) = qtmp(i,j+ig)
479 : ENDDO
480 : ENDDO
481 : !$OMP END PARALLEL DO
482 :
483 : ELSE
484 :
485 : ! Otherwise, call YMAP to regrid in the N-S direction
486 0 : CALL ymap_r4r8(in, jm, sin1, qtmp(1,1+ig), jn, sin2, q2(1,1+ig), ig, iv, &
487 0 : missval=missval )
488 :
489 : ENDIF
490 :
491 0 : END SUBROUTINE Map_A2A_r4r8
492 : !EOC
493 : !------------------------------------------------------------------------------
494 : ! Harmonized Emissions Component (HEMCO) !
495 : !------------------------------------------------------------------------------
496 : !BOP
497 : !
498 : ! !IROUTINE: Map_A2A_r8r4
499 : !
500 : ! !DESCRIPTION: Subroutine MAP\_A2A\_R8R4 is a horizontal arbitrary grid to
501 : ! arbitrary grid conservative high-order mapping regridding routine by
502 : ! S-J Lin. The input data has REAL*8 precision, but the output argument
503 : ! has REAL*4 precision.
504 : !\\
505 : !\\
506 : ! !INTERFACE:
507 : !
508 0 : SUBROUTINE Map_A2A_r8r4( im, jm, lon1, sin1, q1, &
509 0 : in, jn, lon2, sin2, q2, ig, iv, missval)
510 : !
511 : ! !INPUT PARAMETERS:
512 : !
513 : ! Longitude and Latitude dimensions of INPUT grid
514 : INTEGER, INTENT(IN) :: im, jm
515 :
516 : ! Longitude and Latitude dimensions of OUTPUT grid
517 : INTEGER, INTENT(IN) :: in, jn
518 :
519 : ! IG=0: pole to pole;
520 : ! IG=1 J=1 is half-dy north of south pole
521 : INTEGER, INTENT(IN) :: ig
522 :
523 : ! IV=0: Regrid scalar quantity
524 : ! IV=1: Regrid vector quantity
525 : INTEGER, INTENT(IN) :: iv
526 :
527 : ! Longitude edges (degrees) of INPUT and OUTPUT grids
528 : REAL*4, INTENT(IN) :: lon1(im+1), lon2(in+1)
529 :
530 : ! Sine of Latitude Edges (radians) of INPUT and OUTPUT grids
531 : REAL*4, INTENT(IN) :: sin1(jm+1), sin2(jn+1)
532 :
533 : ! Quantity on INPUT grid
534 : REAL*8, INTENT(IN) :: q1(im,jm)
535 : !
536 : ! !OUTPUT PARAMETERS:
537 : !
538 : ! Regridded quantity on OUTPUT grid
539 : REAL*4, INTENT(OUT) :: q2(in,jn)
540 : !
541 : ! !OPTIONAL ARGUMENTS
542 : !
543 : REAL*8, INTENT(IN), OPTIONAL :: missval
544 : !
545 : ! !REMARKS:
546 : ! This routine is overloaded by the MAP_A2A interface.
547 : !
548 : ! !REVISION HISTORY:
549 : ! (1) Original subroutine by S-J Lin. Converted to F90 freeform format
550 : ! and inserted into "Geos3RegridModule" by Bob Yantosca (9/21/00)
551 : ! See https://github.com/geoschem/hemco for complete history
552 : !EOP
553 : !------------------------------------------------------------------------------
554 : !BOC
555 : !
556 : ! !LOCAL VARIABLES:
557 : !
558 : INTEGER :: i,j,k
559 0 : REAL*4 :: qtmp(in,jm)
560 :
561 : ! Init
562 0 : IF ( PRESENT(missval) ) THEN
563 0 : qtmp = real(missval,kind=f4)
564 0 : q2 = real(missval,kind=f4)
565 : ELSE
566 0 : qtmp = miss_r4
567 0 : q2 = miss_r4
568 : ENDIF
569 :
570 : !===================================================================
571 : ! E-W regridding
572 : !===================================================================
573 0 : IF ( im == in .and. &
574 0 : lon1(1) == lon2(1) .and. &
575 : lon1(im+1) == lon2(in+1) ) THEN
576 :
577 : ! Don't call XMAP if both grids have the same # of longitudes
578 : ! but save the input data in the QTMP array
579 : !$OMP PARALLEL DO &
580 : !$OMP DEFAULT( SHARED ) &
581 : !$OMP PRIVATE( I, J )
582 0 : DO j=1,jm-ig
583 0 : DO i=1,im
584 0 : qtmp(i,j+ig) = q1(i,j+ig)
585 : ENDDO
586 : ENDDO
587 : !$OMP END PARALLEL DO
588 :
589 : ELSE
590 :
591 : ! Otherwise, call XMAP to regrid in the E-W direction
592 0 : CALL xmap_r8r4(im, jm-ig, lon1, q1(1,1+ig),in, lon2, qtmp(1,1+ig), &
593 0 : missval=missval )
594 :
595 : ENDIF
596 :
597 : !===================================================================
598 : ! N-S regridding
599 : !===================================================================
600 0 : IF ( jm == jn .and. &
601 0 : sin1(1) == sin2(1) .and. &
602 : sin1(jm+1) == sin2(jn+1) ) THEN
603 :
604 : ! Don't call XMAP if both grids have the same # of longitudes,
605 : ! but assign the value of QTMP to the output Q2 array
606 : !$OMP PARALLEL DO &
607 : !$OMP DEFAULT( SHARED ) &
608 : !$OMP PRIVATE( I, J )
609 0 : DO j=1,jm-ig
610 0 : DO i=1,in
611 0 : q2(i,j+ig) = qtmp(i,j+ig)
612 : ENDDO
613 : ENDDO
614 : !$OMP END PARALLEL DO
615 :
616 : ELSE
617 :
618 : ! Otherwise, call YMAP to regrid in the N-S direction
619 0 : CALL ymap_r4r4(in, jm, sin1, qtmp(1,1+ig), jn, sin2, q2(1,1+ig), ig, iv, &
620 0 : missval=real(missval,kind=f4))
621 :
622 : ENDIF
623 :
624 0 : END SUBROUTINE Map_A2A_r8r4
625 : !EOC
626 : !------------------------------------------------------------------------------
627 : ! Harmonized Emissions Component (HEMCO) !
628 : !------------------------------------------------------------------------------
629 : !BOP
630 : !
631 : ! !IROUTINE: Ymap_r8r8
632 : !
633 : ! !DESCRIPTION: Routine to perform area preserving mapping in N-S from an
634 : ! arbitrary resolution to another. Both the input and output arguments
635 : ! have REAL(fp) precision.
636 : !\\
637 : !\\
638 : ! !INTERFACE:
639 : !
640 0 : SUBROUTINE ymap_r8r8(im, jm, sin1, q1, jn, sin2, q2, ig, iv, missval )
641 : !
642 : ! !INPUT PARAMETERS:
643 : !
644 : ! original E-W dimension
645 : INTEGER, INTENT(IN) :: im
646 :
647 : ! original N-S dimension
648 : INTEGER, INTENT(IN) :: jm
649 :
650 : ! Target N-S dimension
651 : INTEGER, INTENT(IN) :: jn
652 :
653 : ! IG=0: scalars from SP to NP (D-grid v-wind is also IG=0)
654 : ! IG=1: D-grid u-wind
655 : INTEGER, INTENT(IN) :: ig
656 :
657 : ! IV=0: scalar;
658 : ! IV=1: vector
659 : INTEGER, INTENT(IN) :: iv
660 :
661 : ! Original southern edge of the cell sin(lat1)
662 : REAL*8, INTENT(IN) :: sin1(jm+1-ig)
663 :
664 : ! Original data at center of the cell
665 : REAL*8, INTENT(IN) :: q1(im,jm)
666 :
667 : ! Target cell's southern edge sin(lat2)
668 : REAL*8, INTENT(IN) :: sin2(jn+1-ig)
669 : !
670 : ! !OPTIONAL INPUT PARAMETERS:
671 : !
672 : REAL*8, INTENT(IN), OPTIONAL :: missval
673 : !
674 : ! !OUTPUT PARAMETERS:
675 : !
676 : ! Mapped data at the target resolution
677 : REAL*8, INTENT(OUT) :: q2(im,jn)
678 : !
679 : ! !REMARKS:
680 : !
681 : ! sin1 (1) = -1 must be south pole; sin1(jm+1)=1 must be N pole.
682 : !
683 : ! sin1(1) < sin1(2) < sin1(3) < ... < sin1(jm) < sin1(jm+1)
684 : ! sin2(1) < sin2(2) < sin2(3) < ... < sin2(jn) < sin2(jn+1)!
685 : !
686 : ! !AUTHOR:
687 : ! Developer: Prasad Kasibhatla
688 : ! March 6, 2012
689 : !
690 : ! !REVISION HISTORY
691 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
692 : ! See https://github.com/geoschem/hemco for complete history
693 : !EOP
694 : !------------------------------------------------------------------------------
695 : !BOC
696 : !
697 : ! !LOCAL VARIABLES:
698 : !
699 : INTEGER :: i, j0, m, mm, j
700 0 : REAL*8 :: dy1(jm)
701 : REAL*8 :: dy
702 : REAL*8 :: qsum, sum
703 : REAL*8 :: dlat, nlon, miss
704 :
705 : ! YMAP begins here!
706 0 : do j=1,jm-ig
707 0 : dy1(j) = sin1(j+1) - sin1(j)
708 : enddo
709 :
710 : ! Missing value
711 0 : miss = miss_r8
712 0 : if ( present(missval) ) miss=missval
713 :
714 : !===============================================================
715 : ! Area preserving mapping
716 : !===============================================================
717 :
718 : !$OMP PARALLEL DO &
719 : !$OMP DEFAULT( SHARED ) &
720 : !$OMP PRIVATE( I, J0, J, M, QSUM, DLAT, MM, DY )
721 0 : do 1000 i=1,im
722 0 : qsum = 0.0d0
723 0 : dlat = 0.0d0
724 0 : j0 = 1
725 0 : do 555 j=1,jn-ig
726 0 : do 100 m=j0,jm-ig
727 :
728 : !=========================================================
729 : ! locate the southern edge: sin2(i)
730 : !=========================================================
731 0 : if(sin2(j) .ge. sin1(m) .and. sin2(j) .le. sin1(m+1)) then
732 :
733 0 : if(sin2(j+1) .le. sin1(m+1)) then
734 :
735 : ! entire new cell is within the original cell
736 0 : if( abs(q1(i,m)-miss)>tiny_r8 ) q2(i,j)=q1(i,m)
737 : j0 = m
738 : goto 555
739 : else
740 :
741 : ! South most fractional area
742 0 : if( abs(q1(i,m)-miss)>tiny_r8 ) then
743 0 : dlat= sin1(m+1)-sin2(j)
744 0 : qsum=(sin1(m+1)-sin2(j))*q1(i,m)
745 : endif
746 :
747 0 : do mm=m+1,jm-ig
748 :
749 : ! locate the northern edge: sin2(j+1)
750 0 : if(sin2(j+1) .gt. sin1(mm+1) ) then
751 :
752 : ! Whole layer
753 0 : if( abs(q1(i,mm)-miss)>tiny_r8 ) then
754 0 : dlat = dlat + dy1(mm)
755 0 : qsum = qsum + dy1(mm)*q1(i,mm)
756 : endif
757 : else
758 :
759 : ! North most fractional area
760 0 : dy = sin2(j+1)-sin1(mm)
761 0 : if ( abs(q1(i,mm)-miss)>tiny_r8 ) then
762 0 : qsum=qsum+dy*q1(i,mm)
763 0 : dlat=dlat+dy
764 : endif
765 : j0 = mm
766 : goto 123
767 : endif
768 : enddo
769 : goto 123
770 : endif
771 : endif
772 0 : 100 continue
773 : !123 q2(i,j) = qsum / ( sin2(j+1) - sin2(j) )
774 0 : 123 if ( ABS( dlat ) > 0.0d0 ) q2(i,j) = qsum / dlat
775 0 : 555 continue
776 0 : 1000 continue
777 : !$OMP END PARALLEL DO
778 :
779 : !===================================================================
780 : ! Final processing for poles
781 : !===================================================================
782 0 : if ( ig .eq. 0 .and. iv .eq. 0 ) then
783 :
784 : ! South pole
785 0 : if ( sin2(1) .eq. -1.0_fp ) then
786 : sum = 0.e+0_fp
787 : nlon= 0.0d0
788 0 : do i=1,im
789 0 : if(abs(q2(i,1)-miss)>tiny_r8 ) then
790 0 : sum = sum + q2(i,1)
791 0 : nlon= nlon + 1.0d0
792 : endif
793 : enddo
794 :
795 0 : if ( nlon > 0.0d0 ) sum = sum / nlon
796 0 : do i=1,im
797 0 : q2(i,1) = sum
798 : enddo
799 : endif
800 :
801 : ! North pole:
802 0 : if( sin2(jn+1) .eq. 1.0_fp ) then
803 : sum = 0.e+0_fp
804 : nlon= 0.0d0
805 0 : do i=1,im
806 0 : if( abs(q2(i,jn)-miss)>tiny_r8 ) then
807 0 : sum = sum + q2(i,jn)
808 0 : nlon= nlon + 1.0d0
809 : endif
810 : enddo
811 :
812 0 : if ( nlon > 0.0d0 ) sum = sum / DBLE( im )
813 0 : do i=1,im
814 0 : q2(i,jn) = sum
815 : enddo
816 : endif
817 :
818 : endif
819 :
820 0 : END SUBROUTINE ymap_r8r8
821 : !EOC
822 : !------------------------------------------------------------------------------
823 : ! Harmonized Emissions Component (HEMCO) !
824 : !------------------------------------------------------------------------------
825 : !BOP
826 : !
827 : ! !IROUTINE: Ymap_r4r8
828 : !
829 : ! !DESCRIPTION: Routine to perform area preserving mapping in N-S from an
830 : ! arbitrary resolution to another. The input argument has REAL*4 precision
831 : ! but the output argument has REAL(fp) precision.
832 : !\\
833 : !\\
834 : ! !INTERFACE:
835 : !
836 0 : SUBROUTINE ymap_r4r8(im, jm, sin1, q1, jn, sin2, q2, ig, iv, missval)
837 : !
838 : ! !INPUT PARAMETERS:
839 : !
840 :
841 : ! original E-W dimension
842 : INTEGER, INTENT(IN) :: im
843 :
844 : ! original N-S dimension
845 : INTEGER, INTENT(IN) :: jm
846 :
847 : ! Target N-S dimension
848 : INTEGER, INTENT(IN) :: jn
849 :
850 : ! IG=0: scalars from SP to NP (D-grid v-wind is also IG=0)
851 : ! IG=1: D-grid u-wind
852 : INTEGER, INTENT(IN) :: ig
853 :
854 : ! IV=0: scalar;
855 : ! IV=1: vector
856 : INTEGER, INTENT(IN) :: iv
857 :
858 : ! Original southern edge of the cell sin(lat1)
859 : REAL*4, INTENT(IN) :: sin1(jm+1-ig)
860 :
861 : ! Original data at center of the cell
862 : REAL*8, INTENT(IN) :: q1(im,jm)
863 :
864 : ! Target cell's southern edge sin(lat2)
865 : REAL*4, INTENT(IN) :: sin2(jn+1-ig)
866 : !
867 : ! !OPTIONAL INPUT PARAMETERS:
868 : !
869 : REAL*4, INTENT(IN), OPTIONAL :: missval
870 : !
871 : !
872 : ! !OUTPUT PARAMETERS:
873 : !
874 : ! Mapped data at the target resolution
875 : REAL*8, INTENT(OUT) :: q2(im,jn)
876 : !
877 : ! !REMARKS:
878 : !
879 : ! sin1 (1) = -1 must be south pole; sin1(jm+1)=1 must be N pole.
880 : !
881 : ! sin1(1) < sin1(2) < sin1(3) < ... < sin1(jm) < sin1(jm+1)
882 : ! sin2(1) < sin2(2) < sin2(3) < ... < sin2(jn) < sin2(jn+1)!
883 : !
884 : ! !AUTHOR:
885 : ! Developer: Prasad Kasibhatla
886 : ! March 6, 2012
887 : !
888 : ! !REVISION HISTORY
889 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
890 : ! See https://github.com/geoschem/hemco for complete history
891 : !EOP
892 : !------------------------------------------------------------------------------
893 : !BOC
894 : !
895 : ! !LOCAL VARIABLES:
896 : !
897 : INTEGER :: i, j0, m, mm, j
898 0 : REAL*8 :: dy1(jm)
899 : REAL*8 :: dy
900 : REAL*8 :: qsum, dlat, nlon, sum
901 : REAL*4 :: miss
902 :
903 : ! YMAP begins here!
904 0 : do j=1,jm-ig
905 0 : dy1(j) = sin1(j+1) - sin1(j)
906 : enddo
907 :
908 : ! Missing value
909 0 : miss = miss_r4
910 0 : if ( present(missval) ) miss=missval
911 :
912 : !===============================================================
913 : ! Area preserving mapping
914 : !===============================================================
915 :
916 : !$OMP PARALLEL DO &
917 : !$OMP DEFAULT( SHARED ) &
918 : !$OMP PRIVATE( I, J0, J, M, QSUM, DLAT, MM, DY )
919 0 : do 1000 i=1,im
920 0 : qsum = 0.0d0
921 0 : dlat = 0.0d0
922 0 : j0 = 1
923 0 : do 555 j=1,jn-ig
924 0 : do 100 m=j0,jm-ig
925 :
926 : !=========================================================
927 : ! locate the southern edge: sin2(i)
928 : !=========================================================
929 0 : if(sin2(j) .ge. sin1(m) .and. sin2(j) .le. sin1(m+1)) then
930 :
931 0 : if(sin2(j+1) .le. sin1(m+1)) then
932 :
933 : ! entire new cell is within the original cell
934 0 : if ( abs(q1(i,m)-miss)>tiny_r4 ) q2(i,j)=q1(i,m)
935 : j0 = m
936 : goto 555
937 : else
938 :
939 : ! South most fractional area
940 0 : if( abs(q1(i,m)-miss)>tiny_r4 ) then
941 0 : dlat= sin1(m+1)-sin2(j)
942 0 : qsum=(sin1(m+1)-sin2(j))*q1(i,m)
943 : endif
944 :
945 0 : do mm=m+1,jm-ig
946 :
947 : ! locate the northern edge: sin2(j+1)
948 0 : if(sin2(j+1) .gt. sin1(mm+1) ) then
949 :
950 : ! Whole layer
951 0 : if( abs(q1(i,mm)-miss)>tiny_r4 ) then
952 0 : qsum = qsum + dy1(mm)*q1(i,mm)
953 0 : dlat = dlat + dy1(mm)
954 : endif
955 : else
956 :
957 : ! North most fractional area
958 0 : if( abs(q1(i,mm)-miss)>tiny_r4 ) then
959 0 : dy = sin2(j+1)-sin1(mm)
960 0 : qsum=qsum+dy*q1(i,mm)
961 0 : dlat=dlat+dy
962 : endif
963 : j0 = mm
964 : goto 123
965 : endif
966 : enddo
967 : goto 123
968 : endif
969 : endif
970 0 : 100 continue
971 0 : 123 if ( ABS( dlat ) > 0.0d0 ) q2(i,j) = qsum / dlat
972 0 : 555 continue
973 0 : 1000 continue
974 : !$OMP END PARALLEL DO
975 :
976 : !===================================================================
977 : ! Final processing for poles
978 : !===================================================================
979 0 : if ( ig .eq. 0 .and. iv .eq. 0 ) then
980 :
981 : ! South pole:
982 0 : if ( sin2(1) .eq. -1.0_fp ) then
983 : sum = 0.e+0_fp
984 : nlon= 0.0d0
985 0 : do i=1,im
986 0 : if( abs(q2(i,1)-miss)>tiny_r4 ) then
987 0 : sum = sum + q2(i,1)
988 0 : nlon = nlon + 1.0d0
989 : endif
990 : enddo
991 :
992 0 : if ( nlon > 0.0d0 ) sum = sum / nlon
993 0 : do i=1,im
994 0 : q2(i,1) = sum
995 : enddo
996 : endif
997 :
998 : ! North pole:
999 0 : if( sin2(jn+1) .eq. 1.0_fp ) then
1000 : sum = 0.e+0_fp
1001 : nlon = 0.0d0
1002 0 : do i=1,im
1003 0 : if( abs(q2(i,jn)-miss)>tiny_r4 ) then
1004 0 : sum = sum + q2(i,jn)
1005 0 : nlon = nlon + 1.0d0
1006 : endif
1007 : enddo
1008 :
1009 0 : if ( nlon > 0.0d0 ) sum = sum / nlon
1010 0 : do i=1,im
1011 0 : q2(i,jn) = sum
1012 : enddo
1013 : endif
1014 :
1015 : endif
1016 :
1017 0 : END SUBROUTINE ymap_r4r8
1018 : !EOC
1019 : !------------------------------------------------------------------------------
1020 : ! Harmonized Emissions Component (HEMCO) !
1021 : !------------------------------------------------------------------------------
1022 : !BOP
1023 : !
1024 : ! !IROUTINE: Ymap_r8r4
1025 : !
1026 : ! !DESCRIPTION: Routine to perform area preserving mapping in N-S from an
1027 : ! arbitrary resolution to another. The input argument has REAL*8 precision
1028 : ! but the output argument has REAL*4 precision.
1029 : !\\
1030 : !\\
1031 : ! !INTERFACE:
1032 : !
1033 : SUBROUTINE ymap_r8r4(im, jm, sin1, q1, jn, sin2, q2, ig, iv, missval)
1034 : !
1035 : ! !INPUT PARAMETERS:
1036 : !
1037 :
1038 : ! original E-W dimension
1039 : INTEGER, INTENT(IN) :: im
1040 :
1041 : ! original N-S dimension
1042 : INTEGER, INTENT(IN) :: jm
1043 :
1044 : ! Target N-S dimension
1045 : INTEGER, INTENT(IN) :: jn
1046 :
1047 : ! IG=0: scalars from SP to NP (D-grid v-wind is also IG=0)
1048 : ! IG=1: D-grid u-wind
1049 : INTEGER, INTENT(IN) :: ig
1050 :
1051 : ! IV=0: scalar;
1052 : ! IV=1: vector
1053 : INTEGER, INTENT(IN) :: iv
1054 :
1055 : ! Original southern edge of the cell sin(lat1)
1056 : REAL*4, INTENT(IN) :: sin1(jm+1-ig)
1057 :
1058 : ! Original data at center of the cell
1059 : REAL*8, INTENT(IN) :: q1(im,jm)
1060 :
1061 : ! Target cell's southern edge sin(lat2)
1062 : REAL*4, INTENT(IN) :: sin2(jn+1-ig)
1063 : !
1064 : ! !OPTIONAL INPUT PARAMETERS:
1065 : !
1066 : REAL*8, INTENT(IN), OPTIONAL :: missval
1067 : !
1068 : !
1069 : ! !OUTPUT PARAMETERS:
1070 : !
1071 : ! Mapped data at the target resolution
1072 : REAL*4, INTENT(OUT) :: q2(im,jn)
1073 : !
1074 : ! !REMARKS:
1075 : !
1076 : ! sin1 (1) = -1 must be south pole; sin1(jm+1)=1 must be N pole.
1077 : !
1078 : ! sin1(1) < sin1(2) < sin1(3) < ... < sin1(jm) < sin1(jm+1)
1079 : ! sin2(1) < sin2(2) < sin2(3) < ... < sin2(jn) < sin2(jn+1)!
1080 : !
1081 : ! !AUTHOR:
1082 : ! Developer: Prasad Kasibhatla
1083 : ! March 6, 2012
1084 : !
1085 : ! !REVISION HISTORY
1086 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
1087 : ! See https://github.com/geoschem/hemco for complete history
1088 : !EOP
1089 : !------------------------------------------------------------------------------
1090 : !BOC
1091 : !
1092 : ! !LOCAL VARIABLES:
1093 : !
1094 : INTEGER :: i, j0, m, mm, j
1095 : REAL*8 :: dy1(jm)
1096 : REAL*8 :: dy
1097 : REAL*8 :: qsum, sum, dlat
1098 : REAL*8 :: miss
1099 : REAL*4 :: nlon
1100 :
1101 : ! YMAP begins here!
1102 : do j=1,jm-ig
1103 : dy1(j) = sin1(j+1) - sin1(j)
1104 : enddo
1105 :
1106 : ! Missing value
1107 : miss = miss_r8
1108 : if ( present(missval) ) miss=missval
1109 :
1110 : !===============================================================
1111 : ! Area preserving mapping
1112 : !===============================================================
1113 :
1114 : !$OMP PARALLEL DO &
1115 : !$OMP DEFAULT( SHARED ) &
1116 : !$OMP PRIVATE( I, J0, J, M, QSUM, DLAT, MM, DY )
1117 : do 1000 i=1,im
1118 : qsum = 0.0d0
1119 : dlat = 0.0d0
1120 : j0 = 1
1121 : do 555 j=1,jn-ig
1122 : do 100 m=j0,jm-ig
1123 :
1124 : !=========================================================
1125 : ! locate the southern edge: sin2(i)
1126 : !=========================================================
1127 : if(sin2(j) .ge. sin1(m) .and. sin2(j) .le. sin1(m+1)) then
1128 :
1129 : if(sin2(j+1) .le. sin1(m+1)) then
1130 :
1131 : ! entire new cell is within the original cell
1132 : if( abs(q1(i,m)-miss)>tiny_r8 ) q2(i,j)=q1(i,m)
1133 : j0 = m
1134 : goto 555
1135 : else
1136 :
1137 : ! South most fractional area
1138 : if( abs(q1(i,m)-miss)>tiny_r8 ) then
1139 : dlat= sin1(m+1)-sin2(j)
1140 : qsum=(sin1(m+1)-sin2(j))*q1(i,m)
1141 : endif
1142 :
1143 : do mm=m+1,jm-ig
1144 :
1145 : ! locate the northern edge: sin2(j+1)
1146 : if(sin2(j+1) .gt. sin1(mm+1) ) then
1147 :
1148 : ! Whole layer
1149 : if( abs(q1(i,mm)-miss)>tiny_r8 ) then
1150 : qsum = qsum + dy1(mm)*q1(i,mm)
1151 : dlat = dlat + dy1(mm)
1152 : endif
1153 : else
1154 :
1155 : ! North most fractional area
1156 : dy = sin2(j+1)-sin1(mm)
1157 : if( abs(q1(i,mm)-miss)>tiny_r8 ) then
1158 : qsum=qsum+dy*q1(i,mm)
1159 : dlat=dlat+dy
1160 : endif
1161 : j0 = mm
1162 : goto 123
1163 : endif
1164 : enddo
1165 : goto 123
1166 : endif
1167 : endif
1168 : 100 continue
1169 : 123 if ( ABS( dlat ) > 0.0d0 ) q2(i,j) = qsum / dlat
1170 : 555 continue
1171 : 1000 continue
1172 : !$OMP END PARALLEL DO
1173 :
1174 : !===================================================================
1175 : ! Final processing for poles
1176 : !===================================================================
1177 : if ( ig .eq. 0 .and. iv .eq. 0 ) then
1178 :
1179 : ! South pole
1180 : if ( sin2(1) .eq. -1.0_fp ) then
1181 : sum = 0.0_f4
1182 : nlon= 0.0
1183 : do i=1,im
1184 : if( abs(q2(i,1)-miss)>tiny_r8 ) then
1185 : sum = sum + q2(i,1)
1186 : nlon= nlon + 1.0
1187 : endif
1188 : enddo
1189 :
1190 : if ( nlon > 0.0 ) sum = sum / nlon
1191 : do i=1,im
1192 : q2(i,1) = sum
1193 : enddo
1194 : endif
1195 :
1196 : ! North pole:
1197 : if( sin2(jn+1) .eq. 1.0_fp ) then
1198 : sum = 0.0_f4
1199 : nlon= 0.
1200 : do i=1,im
1201 : if( abs(q2(i,jn)-miss)>tiny_r8 ) then
1202 : sum = sum + q2(i,jn)
1203 : nlon= nlon + 1.0
1204 : endif
1205 : enddo
1206 :
1207 : if ( nlon > 0.0 ) sum = sum / nlon
1208 : do i=1,im
1209 : q2(i,jn) = sum
1210 : enddo
1211 : endif
1212 :
1213 : endif
1214 :
1215 : END SUBROUTINE ymap_r8r4
1216 : !EOC
1217 : !------------------------------------------------------------------------------
1218 : ! Harmonized Emissions Component (HEMCO) !
1219 : !------------------------------------------------------------------------------
1220 : !BOP
1221 : !
1222 : ! !IROUTINE: Ymap_r4r4
1223 : !
1224 : ! !DESCRIPTION: Routine to perform area preserving mapping in N-S from an
1225 : ! arbitrary resolution to another. Both the input and output arguments
1226 : ! have REAL(fp) precision.
1227 : !\\
1228 : !\\
1229 : ! !INTERFACE:
1230 : !
1231 0 : SUBROUTINE ymap_r4r4(im, jm, sin1, q1, jn, sin2, q2, ig, iv, missval)
1232 : !
1233 : ! !INPUT PARAMETERS:
1234 : !
1235 :
1236 : ! original E-W dimension
1237 : INTEGER, INTENT(IN) :: im
1238 :
1239 : ! original N-S dimension
1240 : INTEGER, INTENT(IN) :: jm
1241 :
1242 : ! Target N-S dimension
1243 : INTEGER, INTENT(IN) :: jn
1244 :
1245 : ! IG=0: scalars from SP to NP (D-grid v-wind is also IG=0)
1246 : ! IG=1: D-grid u-wind
1247 : INTEGER, INTENT(IN) :: ig
1248 :
1249 : ! IV=0: scalar;
1250 : ! IV=1: vector
1251 : INTEGER, INTENT(IN) :: iv
1252 :
1253 : ! Original southern edge of the cell sin(lat1)
1254 : REAL*4, INTENT(IN) :: sin1(jm+1-ig)
1255 :
1256 : ! Original data at center of the cell
1257 : REAL*4, INTENT(IN) :: q1(im,jm)
1258 :
1259 : ! Target cell's southern edge sin(lat2)
1260 : REAL*4, INTENT(IN) :: sin2(jn+1-ig)
1261 : !
1262 : ! !OPTIONAL INPUT PARAMETERS:
1263 : !
1264 : ! Missing value
1265 : REAL*4, INTENT(IN), OPTIONAL :: missval
1266 : !
1267 : ! !OUTPUT PARAMETERS:
1268 : !
1269 : ! Mapped data at the target resolution
1270 : REAL*4, INTENT(OUT) :: q2(im,jn)
1271 : !
1272 : ! !REMARKS:
1273 : !
1274 : ! sin1 (1) = -1 must be south pole; sin1(jm+1)=1 must be N pole.
1275 : !
1276 : ! sin1(1) < sin1(2) < sin1(3) < ... < sin1(jm) < sin1(jm+1)
1277 : ! sin2(1) < sin2(2) < sin2(3) < ... < sin2(jn) < sin2(jn+1)!
1278 : !
1279 : ! !AUTHOR:
1280 : ! Developer: Prasad Kasibhatla
1281 : ! March 6, 2012
1282 : !
1283 : ! !REVISION HISTORY
1284 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
1285 : ! See https://github.com/geoschem/hemco for complete history
1286 : !EOP
1287 : !------------------------------------------------------------------------------
1288 : !BOC
1289 : !
1290 : ! !LOCAL VARIABLES:
1291 : !
1292 : INTEGER :: i, j0, m, mm, j
1293 0 : REAL*4 :: dy1(jm)
1294 : REAL*4 :: dy
1295 : REAL*4 :: qsum, sum
1296 : REAL*4 :: dlat, nlon, miss
1297 :
1298 : ! YMAP begins here!
1299 0 : do j=1,jm-ig
1300 0 : dy1(j) = sin1(j+1) - sin1(j)
1301 : enddo
1302 :
1303 : ! missing value
1304 0 : miss = miss_r4
1305 0 : if ( present(missval) ) miss = missval
1306 :
1307 : !===============================================================
1308 : ! Area preserving mapping
1309 : !===============================================================
1310 :
1311 : !$OMP PARALLEL DO &
1312 : !$OMP DEFAULT( SHARED ) &
1313 : !$OMP PRIVATE( I, J0, J, M, QSUM, DLAT, MM, DY )
1314 0 : do 1000 i=1,im
1315 0 : qsum = 0.0
1316 0 : dlat = 0.0
1317 0 : j0 = 1
1318 0 : do 555 j=1,jn-ig
1319 0 : do 100 m=j0,jm-ig
1320 :
1321 : !=========================================================
1322 : ! locate the southern edge: sin2(i)
1323 : !=========================================================
1324 0 : if(sin2(j) .ge. sin1(m) .and. sin2(j) .le. sin1(m+1)) then
1325 :
1326 0 : if(sin2(j+1) .le. sin1(m+1)) then
1327 :
1328 : ! entire new cell is within the original cell
1329 0 : if( abs(q1(i,m)-miss)>tiny_r4 ) q2(i,j)=q1(i,m)
1330 : j0 = m
1331 : goto 555
1332 : else
1333 :
1334 : ! South most fractional area
1335 0 : if( abs(q1(i,m)-miss)>tiny_r4 ) then
1336 0 : dlat=sin1(m+1)-sin2(j)
1337 0 : qsum=dlat*q1(i,m)
1338 : endif
1339 :
1340 0 : do mm=m+1,jm-ig
1341 :
1342 : ! locate the northern edge: sin2(j+1)
1343 0 : if(sin2(j+1) .gt. sin1(mm+1) ) then
1344 :
1345 : ! Whole layer
1346 0 : if( abs(q1(i,mm)-miss)>tiny_r4 ) then
1347 0 : qsum = qsum + dy1(mm)*q1(i,mm)
1348 0 : dlat = dlat + dy1(mm)
1349 : endif
1350 : else
1351 :
1352 : ! North most fractional area
1353 0 : if( abs(q1(i,mm)-miss)>tiny_r4 ) then
1354 0 : dy = sin2(j+1)-sin1(mm)
1355 0 : qsum=qsum+dy*q1(i,mm)
1356 0 : dlat=dlat+dy
1357 : endif
1358 : j0 = mm
1359 : goto 123
1360 : endif
1361 : enddo
1362 : goto 123
1363 : endif
1364 : endif
1365 0 : 100 continue
1366 : !123 q2(i,j) = qsum / ( sin2(j+1) - sin2(j) )
1367 0 : 123 if ( ABS( dlat ) > 0.0e0 ) q2(i,j) = qsum / dlat
1368 0 : 555 continue
1369 0 : 1000 continue
1370 : !$OMP END PARALLEL DO
1371 :
1372 : !===================================================================
1373 : ! Final processing for poles
1374 : !===================================================================
1375 0 : if ( ig .eq. 0 .and. iv .eq. 0 ) then
1376 :
1377 : ! South pole
1378 0 : if ( sin2(1) .eq. -1.0_fp ) then
1379 : sum = 0.e+0_fp
1380 : nlon = 0.0
1381 0 : do i=1,im
1382 0 : if( abs(q2(i,1)-miss)>tiny_r4 ) then
1383 0 : sum = sum + q2(i,1)
1384 0 : nlon = nlon + 1.0
1385 : endif
1386 : enddo
1387 :
1388 0 : if ( nlon > 0.0 ) sum = sum / nlon
1389 : !sum = sum / REAL( im, 4 )
1390 0 : do i=1,im
1391 0 : q2(i,1) = sum
1392 : enddo
1393 : endif
1394 :
1395 : ! North pole:
1396 0 : if( sin2(jn+1) .eq. 1.0_fp ) then
1397 : sum = 0.e+0_fp
1398 : nlon= 0.0
1399 0 : do i=1,im
1400 0 : if( abs(q2(i,jn)-miss)>tiny_r4 ) then
1401 0 : sum = sum + q2(i,jn)
1402 0 : nlon = nlon + 1.0
1403 : endif
1404 : enddo
1405 :
1406 : !sum = sum / REAL( im, 4 )
1407 0 : if ( nlon > 0.0 ) sum = sum / nlon
1408 0 : do i=1,im
1409 0 : q2(i,jn) = sum
1410 : enddo
1411 : endif
1412 : endif
1413 :
1414 0 : END SUBROUTINE ymap_r4r4
1415 : !EOC
1416 : !------------------------------------------------------------------------------
1417 : ! Harmonized Emissions Component (HEMCO) !
1418 : !------------------------------------------------------------------------------
1419 : !BOP
1420 : !
1421 : ! !IROUTINE: Xmap_r8r8
1422 : !
1423 : ! !DESCRIPTION: Routine to perform area preserving mapping in E-W from an
1424 : ! arbitrary resolution to another. Both the input and output arguments
1425 : ! have REAL(fp) precision.
1426 : !\\
1427 : !\\
1428 : ! Periodic domain will be assumed, i.e., the eastern wall bounding cell
1429 : ! im is lon1(im+1) = lon1(1); Note the equal sign is true geographysically.
1430 : !\\
1431 : !\\
1432 : ! !INTERFACE:
1433 : !
1434 0 : SUBROUTINE xmap_r8r8(im, jm, lon1, q1, iin, ilon2, iq2, missval)
1435 : !
1436 : ! !INPUT PARAMETERS:
1437 : !
1438 : ! Original E-W dimension
1439 : INTEGER, INTENT(IN) :: im
1440 :
1441 : ! Target E-W dimension
1442 : INTEGER, INTENT(IN) :: iin
1443 :
1444 : ! Original N-S dimension
1445 : INTEGER, INTENT(IN) :: jm
1446 :
1447 : ! Original western edge of the cell
1448 : REAL*8, INTENT(IN) :: lon1(im+1)
1449 :
1450 : ! Original data at center of the cell
1451 : REAL*8, INTENT(IN) :: q1(im,jm)
1452 :
1453 : ! Target cell's western edge
1454 : REAL*8, INTENT(IN), TARGET :: ilon2(iin+1)
1455 : !
1456 : ! !OPTIONAL INPUT PARAMETERS:
1457 : !
1458 : ! Missing value
1459 : REAL*8, INTENT(IN), OPTIONAL :: missval
1460 : !
1461 : ! !OUTPUT PARAMETERS:
1462 : !
1463 : ! Mapped data at the target resolution
1464 : REAL*8, INTENT(OUT), TARGET :: iq2(iin,jm)
1465 : !
1466 : ! !REMARKS:
1467 : ! lon1(1) < lon1(2) < lon1(3) < ... < lon1(im) < lon1(im+1)
1468 : ! lon2(1) < lon2(2) < lon2(3) < ... < lon2(in) < lon2(in+1)
1469 : !
1470 : ! !AUTHOR:
1471 : ! Developer: Prasad Kasibhatla
1472 : ! March 6, 2012
1473 : !
1474 : ! !REVISION HISTORY
1475 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
1476 : ! See https://github.com/geoschem/hemco for complete history
1477 : !EOP
1478 : !------------------------------------------------------------------------------
1479 : !BOC
1480 : !
1481 : ! !LOCAL VARIABLES:
1482 : !
1483 : INTEGER :: i1, i2, i, i0, m, mm, j
1484 0 : REAL*8 :: qtmp(-im:im+im)
1485 0 : REAL*8 :: x1(-im:im+im+1)
1486 0 : REAL*8 :: dx1(-im:im+im)
1487 : REAL*8 :: dx
1488 : REAL*8 :: qsum, dlon
1489 : LOGICAL :: found
1490 :
1491 : ! Update
1492 : INTEGER :: n1, n2
1493 : INTEGER :: in
1494 0 : REAL*8, POINTER :: lon2(:)
1495 0 : REAL*8, POINTER :: q2(:,:)
1496 : REAL*8 :: minlon, maxlon
1497 0 : REAL*8 :: lon1s(im+1)
1498 :
1499 : ! Ghost correction
1500 : Logical :: isGlobal
1501 : Real*8 :: xSpan
1502 :
1503 : ! Missing value
1504 : Real*8 :: miss
1505 :
1506 : ! Initialize pointers
1507 0 : lon2 => NULL()
1508 0 : q2 => NULL()
1509 :
1510 : ! missing value
1511 0 : miss = miss_r8
1512 0 : if ( present(missval) ) miss = missval
1513 :
1514 : ! XMAP begins here!
1515 0 : do i=1,im+1
1516 0 : x1(i) = lon1(i)
1517 : enddo
1518 :
1519 0 : do i=1,im
1520 0 : dx1(i) = x1(i+1) - x1(i)
1521 : enddo
1522 :
1523 : !===================================================================
1524 : ! define minimum and maximum longitude on output grid
1525 : ! to be used. Remapping will be restricted to this
1526 : ! domain. This procedure allows remapping of nested
1527 : ! domains onto larger (e.g. global) domains.
1528 : ! ckeller, 2/11/15).
1529 : !===================================================================
1530 0 : minlon = minval(lon1)
1531 0 : maxlon = maxval(lon1)
1532 :
1533 : ! check for values > 180.0
1534 0 : if(maxlon > 180.0d0) then
1535 0 : lon1s = lon1
1536 0 : do while(maxlon > 180.0d0)
1537 0 : WHERE(lon1s > 180.0d0) lon1s = lon1s - 360.0d0
1538 0 : minlon = minval(lon1s)
1539 0 : maxlon = maxval(lon1s)
1540 : enddo
1541 : endif
1542 :
1543 : ! maxlon must represent the easter edge of the grid:
1544 0 : maxlon = maxlon + ( lon1(im+1)-lon1(im) )
1545 :
1546 : ! Reduce input grid
1547 0 : n1 = 1
1548 0 : n2 = iin+1
1549 0 : do i=1,iin+1
1550 0 : if ( ilon2(i) < minlon ) n1 = i
1551 0 : if ( ilon2(iin+2-i) > maxlon ) n2 = iin+2-i
1552 : enddo
1553 0 : in = n2 - n1
1554 0 : lon2 => ilon2(n1:n2)
1555 0 : q2 => iq2(n1:(n2-1),:)
1556 :
1557 : ! if there is no overlap between original grid and output grid
1558 : ! reduced will be zero and missing values should be returned
1559 0 : if ( in .eq. 0 ) then
1560 0 : iq2 = missval
1561 0 : lon2 => NULL()
1562 0 : q2 => NULL()
1563 : return
1564 : endif
1565 :
1566 : ! Periodic BC only valid if the variable is "global"
1567 0 : xSpan = x1(im+1)-x1(1)
1568 0 : isGlobal = ((xSpan.ge.355.0).and.(xSpan.le.365.0))
1569 :
1570 : !===================================================================
1571 : ! check to see if ghosting is necessary
1572 : ! Western edge:
1573 : !===================================================================
1574 0 : found = .false.
1575 0 : i1 = 1
1576 : do while ( .not. found )
1577 0 : if( lon2(1) .ge. x1(i1) ) then
1578 : found = .true.
1579 : else
1580 0 : i1 = i1 - 1
1581 0 : if (i1 .lt. -im) then
1582 0 : write(6,*) 'Failed in Xmap_R8R8 (regrid_a2a_mod.F90)'
1583 0 : stop
1584 : else
1585 0 : x1(i1) = x1(i1+1) - dx1(im+i1)
1586 0 : dx1(i1) = dx1(im+i1)
1587 : endif
1588 : endif
1589 : enddo
1590 :
1591 : !===================================================================
1592 : ! Eastern edge:
1593 : !===================================================================
1594 : found = .false.
1595 : i2 = im+1
1596 : do while ( .not. found )
1597 0 : if( lon2(in+1) .le. x1(i2) ) then
1598 : found = .true.
1599 : else
1600 0 : i2 = i2 + 1
1601 0 : if (i2 .gt. 2*im) then
1602 0 : write(6,*) 'Failed in Xmap_R8R8 (regrid_a2a_mod.F90)'
1603 0 : stop
1604 : else
1605 0 : dx1(i2-1) = dx1(i2-1-im)
1606 0 : x1(i2) = x1(i2-1) + dx1(i2-1)
1607 : endif
1608 : endif
1609 : enddo
1610 :
1611 : !$OMP PARALLEL DO &
1612 : !$OMP DEFAULT( SHARED ) &
1613 : !$OMP PRIVATE( J, QTMP, I, I0, M, QSUM, DLON, MM, DX )
1614 0 : do 1000 j=1,jm
1615 :
1616 : !=================================================================
1617 : ! Area preserving mapping
1618 : !================================================================
1619 :
1620 0 : qtmp(:) = 0.0d0
1621 0 : do i=1,im
1622 0 : qtmp(i)=q1(i,j)
1623 : enddo
1624 :
1625 : ! SDE 2017-01-07
1626 : ! Only have shadow regions if we are on a global grid. Otherwise, we
1627 : ! should keep the zero boundary conditions.
1628 0 : If (isGlobal) Then
1629 0 : qtmp(0)=q1(im,j)
1630 0 : qtmp(im+1)=q1(1,j)
1631 :
1632 : ! check to see if ghosting is necessary
1633 : ! Western edge
1634 0 : if ( i1 .le. 0 ) then
1635 0 : do i=i1,0
1636 0 : qtmp(i) = qtmp(im+i)
1637 : enddo
1638 : endif
1639 :
1640 : ! Eastern edge:
1641 0 : if ( i2 .gt. im+1 ) then
1642 0 : do i=im+1,i2-1
1643 0 : qtmp(i) = qtmp(i-im)
1644 : enddo
1645 : endif
1646 : End If
1647 :
1648 : i0 = i1
1649 :
1650 0 : do 555 i=1,in
1651 0 : do 100 m=i0,i2-1
1652 :
1653 : !=============================================================
1654 : ! locate the western edge: lon2(i)
1655 : !=============================================================
1656 0 : if(lon2(i) .ge. x1(m) .and. lon2(i) .le. x1(m+1)) then
1657 :
1658 0 : if(lon2(i+1) .le. x1(m+1)) then
1659 :
1660 : ! entire new grid is within the original grid
1661 0 : if( abs(qtmp(m)-miss)>tiny_r8 ) q2(i,j)=qtmp(m)
1662 : i0 = m
1663 : goto 555
1664 : else
1665 :
1666 : ! Left most fractional area
1667 0 : if( abs(qtmp(m)-miss)>tiny_r8 ) then
1668 0 : qsum=(x1(m+1)-lon2(i))*qtmp(m)
1669 : dlon= x1(m+1)-lon2(i)
1670 : else
1671 : qsum = 0.0d0
1672 : dlon = 0.0d0
1673 : endif
1674 0 : do mm=m+1,i2-1
1675 :
1676 : ! locate the eastern edge: lon2(i+1)
1677 0 : if(lon2(i+1) .gt. x1(mm+1) ) then
1678 :
1679 : ! Whole layer
1680 0 : if( abs(qtmp(mm)-miss)>tiny_r8 ) then
1681 0 : qsum = qsum + dx1(mm)*qtmp(mm)
1682 0 : dlon = dlon + dx1(mm)
1683 : endif
1684 : else
1685 : ! Right most fractional area
1686 0 : if( abs(qtmp(mm)-miss)>tiny_r8 ) then
1687 0 : dx = lon2(i+1)-x1(mm)
1688 0 : qsum=qsum+dx*qtmp(mm)
1689 0 : dlon=dlon+dx
1690 : endif
1691 : i0 = mm
1692 : goto 123
1693 : endif
1694 : enddo
1695 : goto 123
1696 : endif
1697 : endif
1698 0 : 100 continue
1699 0 : 123 if ( ABS( dlon ) > 0.0d0 ) q2(i,j) = qsum / dlon
1700 0 : 555 continue
1701 0 : 1000 continue
1702 : !$OMP END PARALLEL DO
1703 :
1704 : ! Cleanup
1705 0 : lon2 => NULL()
1706 0 : q2 => NULL()
1707 :
1708 0 : END SUBROUTINE xmap_r8r8
1709 : !EOC
1710 : !------------------------------------------------------------------------------
1711 : ! Harmonized Emissions Component (HEMCO) !
1712 : !------------------------------------------------------------------------------
1713 : !BOP
1714 : !
1715 : ! !IROUTINE: Xmap_r4r4
1716 : !
1717 : ! !DESCRIPTION: Routine to perform area preserving mapping in E-W from an
1718 : ! arbitrary resolution to another. Both the input and output arguments
1719 : ! have REAL*4 precision.
1720 : !\\
1721 : !\\
1722 : ! Periodic domain will be assumed, i.e., the eastern wall bounding cell
1723 : ! im is lon1(im+1) = lon1(1); Note the equal sign is true geographysically.
1724 : !\\
1725 : !\\
1726 : ! !INTERFACE:
1727 : !
1728 0 : SUBROUTINE xmap_r4r4(im, jm, lon1, q1, iin, ilon2, iq2, missval)
1729 : !
1730 : ! !INPUT PARAMETERS:
1731 : !
1732 : ! Original E-W dimension
1733 : INTEGER, INTENT(IN) :: im
1734 :
1735 : ! Target E-W dimension
1736 : INTEGER, INTENT(IN) :: iin
1737 :
1738 : ! Original N-S dimension
1739 : INTEGER, INTENT(IN) :: jm
1740 :
1741 : ! Original western edge of the cell
1742 : REAL*4, INTENT(IN) :: lon1(im+1)
1743 :
1744 : ! Original data at center of the cell
1745 : REAL*4, INTENT(IN) :: q1(im,jm)
1746 :
1747 : ! Target cell's western edge
1748 : REAL*4, INTENT(IN), TARGET :: ilon2(iin+1)
1749 : !
1750 : ! !OPTIONAL INPUT PARAMETERS:
1751 : !
1752 : ! Missing value
1753 : REAL*4, INTENT(IN), OPTIONAL :: missval
1754 : !
1755 : ! !OUTPUT PARAMETERS:
1756 : !
1757 : ! Mapped data at the target resolution
1758 : REAL*4, INTENT(OUT), TARGET :: iq2(iin,jm)
1759 : !
1760 : ! !REMARKS:
1761 : ! lon1(1) < lon1(2) < lon1(3) < ... < lon1(im) < lon1(im+1)
1762 : ! lon2(1) < lon2(2) < lon2(3) < ... < lon2(in) < lon2(in+1)
1763 : !
1764 : ! !AUTHOR:
1765 : ! Developer: Prasad Kasibhatla
1766 : ! March 6, 2012
1767 : !
1768 : ! !REVISION HISTORY
1769 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
1770 : ! See https://github.com/geoschem/hemco for complete history
1771 : !EOP
1772 : !------------------------------------------------------------------------------
1773 : !BOC
1774 : !
1775 : ! !LOCAL VARIABLES:
1776 : !
1777 : INTEGER :: i1, i2, i, i0, m, mm, j
1778 0 : REAL*4 :: qtmp(-im:im+im)
1779 0 : REAL*4 :: x1(-im:im+im+1)
1780 0 : REAL*4 :: dx1(-im:im+im)
1781 : REAL*4 :: dx
1782 : REAL*4 :: qsum, dlon
1783 : LOGICAL :: found
1784 :
1785 : ! Update
1786 : INTEGER :: n1, n2
1787 : INTEGER :: in
1788 0 : REAL*4, POINTER :: lon2(:)
1789 0 : REAL*4, POINTER :: q2(:,:)
1790 : REAL*4 :: minlon, maxlon
1791 0 : REAL*4 :: lon1s(im+1)
1792 :
1793 : ! Ghost correction
1794 : Logical :: isGlobal
1795 : Real*4 :: xSpan
1796 :
1797 : ! Missing value
1798 : REAL*4 :: miss
1799 :
1800 : ! Initialize
1801 0 : lon2 => NULL()
1802 0 : q2 => NULL()
1803 :
1804 : ! Missing value
1805 0 : miss = miss_r4
1806 0 : if ( present(missval) ) miss = missval
1807 :
1808 : ! XMAP begins here!
1809 0 : do i=1,im+1
1810 0 : x1(i) = lon1(i)
1811 : enddo
1812 :
1813 0 : do i=1,im
1814 0 : dx1(i) = x1(i+1) - x1(i)
1815 : enddo
1816 :
1817 : !===================================================================
1818 : ! define minimum and maximum longitude on output grid
1819 : ! to be used. Remapping will be restricted to this
1820 : ! domain. This procedure allows remapping of nested
1821 : ! domains onto larger (e.g. global) domains.
1822 : ! ckeller, (2/11/15).
1823 : !===================================================================
1824 0 : minlon = minval(lon1)
1825 0 : maxlon = maxval(lon1)
1826 :
1827 : ! check for values > 180.0
1828 0 : if(maxlon > 180.0) then
1829 0 : lon1s = lon1
1830 0 : do while(maxlon > 180.0)
1831 0 : WHERE(lon1s > 180.0) lon1s = lon1s - 360.0
1832 0 : minlon = minval(lon1s)
1833 0 : maxlon = maxval(lon1s)
1834 : enddo
1835 : endif
1836 :
1837 : ! maxlon must represent the easter edge of the grid:
1838 0 : maxlon = maxlon + ( lon1(im+1)-lon1(im) )
1839 :
1840 : ! Reduce output grid
1841 0 : n1 = 1
1842 0 : n2 = iin+1
1843 0 : do i=1,iin+1
1844 0 : if ( ilon2(i) < minlon ) n1 = i
1845 0 : if ( ilon2(iin+2-i) > maxlon ) n2 = iin+2-i
1846 : enddo
1847 0 : in = n2 - n1
1848 0 : lon2 => ilon2(n1:n2)
1849 0 : q2 => iq2(n1:(n2-1),:)
1850 :
1851 : ! if there is no overlap between original grid and output grid
1852 : ! reduced will be zero and missing values should be returned
1853 0 : if ( in .eq. 0 ) then
1854 0 : iq2 = missval
1855 0 : lon2 => NULL()
1856 0 : q2 => NULL()
1857 : return
1858 : endif
1859 :
1860 : ! shadow variables to selected range
1861 : ! Periodic BC only valid if the variable is "global"
1862 0 : xSpan = x1(im+1)-x1(1)
1863 0 : isGlobal = ((xSpan.ge.355.0).and.(xSpan.le.365.0))
1864 :
1865 : !===================================================================
1866 : ! check to see if ghosting is necessary
1867 : ! Western edge:
1868 : !===================================================================
1869 0 : found = .false.
1870 0 : i1 = 1
1871 : do while ( .not. found )
1872 0 : if( lon2(1) .ge. x1(i1) ) then
1873 : found = .true.
1874 : else
1875 0 : i1 = i1 - 1
1876 0 : if (i1 .lt. -im) then
1877 0 : write(6,*) 'Failed in Xmap_R4R4 (regrid_a2a_mod.F90)'
1878 0 : stop
1879 : else
1880 0 : x1(i1) = x1(i1+1) - dx1(im+i1)
1881 0 : dx1(i1) = dx1(im+i1)
1882 : endif
1883 : endif
1884 : enddo
1885 :
1886 : !===================================================================
1887 : ! Eastern edge:
1888 : !===================================================================
1889 : found = .false.
1890 : i2 = im+1
1891 : do while ( .not. found )
1892 0 : if( lon2(in+1) .le. x1(i2) ) then
1893 : found = .true.
1894 : else
1895 0 : i2 = i2 + 1
1896 0 : if (i2 .gt. 2*im) then
1897 0 : write(6,*) 'Failed in Xmap_R4R4 (regrid_a2a_mod.F90)'
1898 0 : stop
1899 : else
1900 0 : dx1(i2-1) = dx1(i2-1-im)
1901 0 : x1(i2) = x1(i2-1) + dx1(i2-1)
1902 : endif
1903 : endif
1904 : enddo
1905 :
1906 : !$OMP PARALLEL DO &
1907 : !$OMP DEFAULT( SHARED ) &
1908 : !$OMP PRIVATE( J, QTMP, I, I0, M, QSUM, DLON, MM, DX )
1909 0 : do 1000 j=1,jm
1910 :
1911 : !=================================================================
1912 : ! Area preserving mapping
1913 : !================================================================
1914 :
1915 0 : qtmp(:) = 0.0
1916 0 : do i=1,im
1917 0 : qtmp(i)=q1(i,j)
1918 : enddo
1919 :
1920 : ! SDE 2017-01-07
1921 : ! Only have shadow regions if we are on a global grid. Otherwise, we
1922 : ! should keep the zero boundary conditions.
1923 0 : If (isGlobal) Then
1924 0 : qtmp(0)=q1(im,j)
1925 0 : qtmp(im+1)=q1(1,j)
1926 :
1927 : ! check to see if ghosting is necessary
1928 : ! Western edge
1929 0 : if ( i1 .le. 0 ) then
1930 0 : do i=i1,0
1931 0 : qtmp(i) = qtmp(im+i)
1932 : enddo
1933 : endif
1934 :
1935 : ! Eastern edge:
1936 0 : if ( i2 .gt. im+1 ) then
1937 0 : do i=im+1,i2-1
1938 0 : qtmp(i) = qtmp(i-im)
1939 : enddo
1940 : endif
1941 : End If
1942 :
1943 : i0 = i1
1944 :
1945 0 : do 555 i=1,in
1946 0 : do 100 m=i0,i2-1
1947 :
1948 : !=============================================================
1949 : ! locate the western edge: lon2(i)
1950 : !=============================================================
1951 0 : if(lon2(i) .ge. x1(m) .and. lon2(i) .le. x1(m+1)) then
1952 :
1953 0 : if(lon2(i+1) .le. x1(m+1)) then
1954 :
1955 : ! entire new grid is within the original grid
1956 0 : if ( abs(qtmp(m)-miss)>tiny_r4 ) q2(i,j)=qtmp(m)
1957 : i0 = m
1958 : goto 555
1959 : else
1960 :
1961 : ! Left most fractional area
1962 0 : if( abs(qtmp(m)-miss)>tiny_r4 ) then
1963 0 : dlon=x1(m+1)-lon2(i)
1964 0 : qsum=dlon*qtmp(m)
1965 : else
1966 : dlon=0.0
1967 : qsum=0.0
1968 : endif
1969 0 : do mm=m+1,i2-1
1970 :
1971 : ! locate the eastern edge: lon2(i+1)
1972 0 : if(lon2(i+1) .gt. x1(mm+1) ) then
1973 :
1974 : ! Whole layer
1975 0 : if( abs(qtmp(mm)-miss)>tiny_r4 ) then
1976 0 : qsum = qsum + dx1(mm)*qtmp(mm)
1977 0 : dlon = dlon + dx1(mm)
1978 : endif
1979 :
1980 : else
1981 : ! Right most fractional area
1982 0 : if( abs(qtmp(mm)-miss)>tiny_r4 ) then
1983 0 : dx = lon2(i+1)-x1(mm)
1984 0 : qsum=qsum+dx*qtmp(mm)
1985 0 : dlon=dlon+dx
1986 : endif
1987 : i0 = mm
1988 : goto 123
1989 : endif
1990 : enddo
1991 : goto 123
1992 : endif
1993 : endif
1994 0 : 100 continue
1995 0 : 123 if( ABS( dlon ) > 0.0e0 ) q2(i,j) = qsum / dlon
1996 0 : 555 continue
1997 0 : 1000 continue
1998 : !$OMP END PARALLEL DO
1999 :
2000 : ! Cleanup
2001 0 : lon2 => NULL()
2002 0 : q2 => NULL()
2003 :
2004 0 : END SUBROUTINE xmap_r4r4
2005 : !EOC
2006 : !------------------------------------------------------------------------------
2007 : ! Harmonized Emissions Component (HEMCO) !
2008 : !------------------------------------------------------------------------------
2009 : !BOP
2010 : !
2011 : ! !IROUTINE: Xmap_r4r8
2012 : !
2013 : ! !DESCRIPTION: Routine to perform area preserving mapping in E-W from an
2014 : ! arbitrary resolution to another. The input argument has REAL*4 precision
2015 : ! but the output argument has REAL(fp) precision.
2016 : !\\
2017 : !\\
2018 : ! Periodic domain will be assumed, i.e., the eastern wall bounding cell
2019 : ! im is lon1(im+1) = lon1(1); Note the equal sign is true geographysically.
2020 : !\\
2021 : !\\
2022 : ! !INTERFACE:
2023 : !
2024 0 : SUBROUTINE xmap_r4r8(im, jm, lon1, q1, iin, ilon2, iq2, missval)
2025 : !
2026 : ! !INPUT PARAMETERS:
2027 : !
2028 : ! Original E-W dimension
2029 : INTEGER, INTENT(IN) :: im
2030 :
2031 : ! Target E-W dimension
2032 : INTEGER, INTENT(IN) :: iin
2033 :
2034 : ! Original N-S dimension
2035 : INTEGER, INTENT(IN) :: jm
2036 :
2037 : ! Original western edge of the cell
2038 : REAL*4, INTENT(IN) :: lon1(im+1)
2039 :
2040 : ! Original data at center of the cell
2041 : REAL*4, INTENT(IN) :: q1(im,jm)
2042 :
2043 : ! Target cell's western edge
2044 : REAL*4, INTENT(IN), TARGET :: ilon2(iin+1)
2045 : !
2046 : ! !OPTIONAL INPUT PARAMETERS:
2047 : !
2048 : REAL*4, INTENT(IN), OPTIONAL :: missval
2049 : !
2050 : ! !OUTPUT PARAMETERS:
2051 : !
2052 : ! Mapped data at the target resolution
2053 : REAL*8, INTENT(OUT), TARGET :: iq2(iin,jm)
2054 : !
2055 : ! !REMARKS:
2056 : ! lon1(1) < lon1(2) < lon1(3) < ... < lon1(im) < lon1(im+1)
2057 : ! lon2(1) < lon2(2) < lon2(3) < ... < lon2(in) < lon2(in+1)
2058 : !
2059 : ! !AUTHOR:
2060 : ! Developer: Prasad Kasibhatla
2061 : ! March 6, 2012
2062 : !
2063 : ! !REVISION HISTORY
2064 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
2065 : ! See https://github.com/geoschem/hemco for complete history
2066 : !EOP
2067 : !------------------------------------------------------------------------------
2068 : !BOC
2069 : !
2070 : ! !LOCAL VARIABLES:
2071 : !
2072 : INTEGER :: i1, i2, i, i0, m, mm, j
2073 0 : REAL*8 :: qtmp(-im:im+im)
2074 0 : REAL*8 :: x1(-im:im+im+1)
2075 0 : REAL*8 :: dx1(-im:im+im)
2076 : REAL*8 :: dx
2077 : REAL*8 :: qsum, dlon
2078 : LOGICAL :: found
2079 :
2080 : ! Update
2081 : INTEGER :: n1, n2
2082 : INTEGER :: in
2083 0 : REAL*4, POINTER :: lon2(:)
2084 0 : REAL*8, POINTER :: q2(:,:)
2085 : REAL*4 :: minlon, maxlon
2086 0 : REAL*4 :: lon1s(im+1)
2087 :
2088 : ! Ghost correction
2089 : Logical :: isGlobal
2090 : Real*8 :: xSpan
2091 :
2092 : ! Missing value
2093 : REAL*4 :: miss
2094 :
2095 : ! Initialize
2096 0 : lon2 => NULL()
2097 0 : q2 => NULL()
2098 :
2099 : ! Missing value
2100 0 : miss = miss_r4
2101 0 : if ( present(missval) ) miss = missval
2102 :
2103 : ! XMAP begins here!
2104 0 : do i=1,im+1
2105 0 : x1(i) = lon1(i)
2106 : enddo
2107 :
2108 0 : do i=1,im
2109 0 : dx1(i) = x1(i+1) - x1(i)
2110 : enddo
2111 :
2112 : !===================================================================
2113 : ! define minimum and maximum longitude on output grid
2114 : ! to be used. Remapping will be restricted to this
2115 : ! domain. This procedure allows remapping of nested
2116 : ! domains onto larger (e.g. global) domains.
2117 : ! ckeller, 2/11/15).
2118 : !===================================================================
2119 0 : minlon = minval(lon1)
2120 0 : maxlon = maxval(lon1)
2121 :
2122 : ! check for values > 180.0
2123 0 : if(maxlon > 180.0) then
2124 0 : lon1s = lon1
2125 0 : do while(maxlon > 180.0)
2126 0 : WHERE(lon1s > 180.0) lon1s = lon1s - 360.0
2127 0 : minlon = minval(lon1s)
2128 0 : maxlon = maxval(lon1s)
2129 : enddo
2130 : endif
2131 :
2132 : ! maxlon must represent the easter edge of the grid:
2133 0 : maxlon = maxlon + ( lon1(im+1)-lon1(im) )
2134 :
2135 : ! Reduce input grid
2136 0 : n1 = 1
2137 0 : n2 = iin+1
2138 0 : do i=1,iin+1
2139 0 : if ( ilon2(i) < minlon ) n1 = i
2140 0 : if ( ilon2(iin+2-i) > maxlon ) n2 = iin+2-i
2141 : enddo
2142 0 : in = n2 - n1
2143 0 : lon2 => ilon2(n1:n2)
2144 0 : q2 => iq2(n1:(n2-1),:)
2145 :
2146 : ! Periodic BC only valid if the variable is "global"
2147 0 : xSpan = x1(im+1)-x1(1)
2148 0 : isGlobal = ((xSpan.ge.355.0).and.(xSpan.le.365.0))
2149 :
2150 : !===================================================================
2151 : ! check to see if ghosting is necessary
2152 : ! Western edge:
2153 : !===================================================================
2154 0 : found = .false.
2155 0 : i1 = 1
2156 : do while ( .not. found )
2157 0 : if( lon2(1) .ge. x1(i1) ) then
2158 : found = .true.
2159 : else
2160 0 : i1 = i1 - 1
2161 0 : if (i1 .lt. -im) then
2162 0 : write(6,*) 'Failed in Xmap_R4R8 (regrid_a2a_mod.F90)'
2163 0 : stop
2164 : else
2165 0 : x1(i1) = x1(i1+1) - dx1(im+i1)
2166 0 : dx1(i1) = dx1(im+i1)
2167 : endif
2168 : endif
2169 : enddo
2170 :
2171 : !===================================================================
2172 : ! Eastern edge:
2173 : !===================================================================
2174 : found = .false.
2175 : i2 = im+1
2176 : do while ( .not. found )
2177 0 : if( lon2(in+1) .le. x1(i2) ) then
2178 : found = .true.
2179 : else
2180 0 : i2 = i2 + 1
2181 0 : if (i2 .gt. 2*im) then
2182 0 : write(6,*) 'Failed in Xmap_R4R8 (regrid_a2a_mod.F90)'
2183 0 : stop
2184 : else
2185 0 : dx1(i2-1) = dx1(i2-1-im)
2186 0 : x1(i2) = x1(i2-1) + dx1(i2-1)
2187 : endif
2188 : endif
2189 : enddo
2190 :
2191 : !$OMP PARALLEL DO &
2192 : !$OMP DEFAULT( SHARED ) &
2193 : !$OMP PRIVATE( J, QTMP, I, I0, M, QSUM, DLON, MM, DX )
2194 0 : do 1000 j=1,jm
2195 :
2196 : !=================================================================
2197 : ! Area preserving mapping
2198 : !================================================================
2199 :
2200 0 : qtmp(:) = 0.0d0
2201 0 : do i=1,im
2202 0 : qtmp(i)=q1(i,j)
2203 : enddo
2204 :
2205 : ! SDE 2017-01-07
2206 : ! Only have shadow regions if we are on a global grid. Otherwise, we
2207 : ! should keep the zero boundary conditions.
2208 0 : If (isGlobal) Then
2209 0 : qtmp(0)=q1(im,j)
2210 0 : qtmp(im+1)=q1(1,j)
2211 :
2212 : ! check to see if ghosting is necessary
2213 : ! Western edge
2214 0 : if ( i1 .le. 0 ) then
2215 0 : do i=i1,0
2216 0 : qtmp(i) = qtmp(im+i)
2217 : enddo
2218 : endif
2219 :
2220 : ! Eastern edge:
2221 0 : if ( i2 .gt. im+1 ) then
2222 0 : do i=im+1,i2-1
2223 0 : qtmp(i) = qtmp(i-im)
2224 : enddo
2225 : endif
2226 : End If
2227 :
2228 : i0 = i1
2229 :
2230 0 : do 555 i=1,in
2231 0 : do 100 m=i0,i2-1
2232 :
2233 : !=============================================================
2234 : ! locate the western edge: lon2(i)
2235 : !=============================================================
2236 0 : if(lon2(i) .ge. x1(m) .and. lon2(i) .le. x1(m+1)) then
2237 :
2238 0 : if(lon2(i+1) .le. x1(m+1)) then
2239 :
2240 : ! entire new grid is within the original grid
2241 0 : if( abs(qtmp(m)-miss)>tiny_r4 ) q2(i,j)=qtmp(m)
2242 : i0 = m
2243 : goto 555
2244 : else
2245 :
2246 : ! Left most fractional area
2247 0 : if( abs(qtmp(m)-miss)>tiny_r4 ) then
2248 0 : qsum=(x1(m+1)-lon2(i))*qtmp(m)
2249 : dlon= x1(m+1)-lon2(i)
2250 : else
2251 : qsum=0.0d0
2252 : dlon=0.0d0
2253 : endif
2254 0 : do mm=m+1,i2-1
2255 :
2256 : ! locate the eastern edge: lon2(i+1)
2257 0 : if(lon2(i+1) .gt. x1(mm+1) ) then
2258 :
2259 : ! Whole layer
2260 0 : if( abs(qtmp(mm)-miss)>tiny_r4 ) then
2261 0 : qsum = qsum + dx1(mm)*qtmp(mm)
2262 0 : dlon = dlon + dx1(mm)
2263 : endif
2264 : else
2265 : ! Right most fractional area
2266 0 : if( abs(qtmp(mm)-miss)>tiny_r4 ) then
2267 0 : dx = lon2(i+1)-x1(mm)
2268 0 : qsum=qsum+dx*qtmp(mm)
2269 0 : dlon=dlon+dx
2270 : endif
2271 : i0 = mm
2272 : goto 123
2273 : endif
2274 : enddo
2275 : goto 123
2276 : endif
2277 : endif
2278 0 : 100 continue
2279 0 : 123 if ( ABS( dlon ) > 0.0d0 ) q2(i,j) = qsum / dlon
2280 0 : 555 continue
2281 0 : 1000 continue
2282 : !$OMP END PARALLEL DO
2283 :
2284 : ! Cleanup
2285 0 : lon2 => NULL()
2286 0 : q2 => NULL()
2287 :
2288 0 : END SUBROUTINE xmap_r4r8
2289 : !EOC
2290 : !------------------------------------------------------------------------------
2291 : ! Harmonized Emissions Component (HEMCO) !
2292 : !------------------------------------------------------------------------------
2293 : !BOP
2294 : !
2295 : ! !IROUTINE: Xmap_r8r4
2296 : !
2297 : ! !DESCRIPTION: Routine to perform area preserving mapping in E-W from an
2298 : ! arbitrary resolution to another. The input argument has REAL*8 precision
2299 : ! but the output argument has REAL*4 precision.
2300 : !\\
2301 : !\\
2302 : ! Periodic domain will be assumed, i.e., the eastern wall bounding cell
2303 : ! im is lon1(im+1) = lon1(1); Note the equal sign is true geographysically.
2304 : !\\
2305 : !\\
2306 : ! !INTERFACE:
2307 : !
2308 0 : SUBROUTINE xmap_r8r4(im, jm, lon1, q1, iin, ilon2, iq2, missval)
2309 : !
2310 : ! !INPUT PARAMETERS:
2311 : !
2312 : ! Original E-W dimension
2313 : INTEGER, INTENT(IN) :: im
2314 :
2315 : ! Target E-W dimension
2316 : INTEGER, INTENT(IN) :: iin
2317 :
2318 : ! Original N-S dimension
2319 : INTEGER, INTENT(IN) :: jm
2320 :
2321 : ! Original western edge of the cell
2322 : REAL*4, INTENT(IN) :: lon1(im+1)
2323 :
2324 : ! Original data at center of the cell
2325 : REAL*8, INTENT(IN) :: q1(im,jm)
2326 :
2327 : ! Target cell's western edge
2328 : REAL*4, INTENT(IN), TARGET :: ilon2(iin+1)
2329 : !
2330 : ! !OPTIONAL INPUT PARAMETERS:
2331 : !
2332 : REAL*8, INTENT(IN), OPTIONAL :: missval
2333 : !
2334 : ! !OUTPUT PARAMETERS:
2335 : !
2336 : ! Mapped data at the target resolution
2337 : REAL*4, INTENT(OUT), TARGET :: iq2(iin,jm)
2338 : !
2339 : ! !REMARKS:
2340 : ! lon1(1) < lon1(2) < lon1(3) < ... < lon1(im) < lon1(im+1)
2341 : ! lon2(1) < lon2(2) < lon2(3) < ... < lon2(in) < lon2(in+1)
2342 : !
2343 : ! !AUTHOR:
2344 : ! Developer: Prasad Kasibhatla
2345 : ! March 6, 2012
2346 : !
2347 : ! !REVISION HISTORY
2348 : ! 06 Mar 2012 - P. Kasibhatla - Initial version
2349 : ! See https://github.com/geoschem/hemco for complete history
2350 : !EOP
2351 : !------------------------------------------------------------------------------
2352 : !BOC
2353 : !
2354 : ! !LOCAL VARIABLES:
2355 : !
2356 : INTEGER :: i1, i2, i, i0, m, mm, j
2357 0 : REAL*4 :: qtmp(-im:im+im)
2358 0 : REAL*4 :: x1(-im:im+im+1)
2359 0 : REAL*4 :: dx1(-im:im+im)
2360 : REAL*4 :: dx
2361 : REAL*4 :: qsum, dlon
2362 : LOGICAL :: found
2363 :
2364 : ! Update
2365 : INTEGER :: n1, n2
2366 : INTEGER :: in
2367 0 : REAL*4, POINTER :: lon2(:)
2368 0 : REAL*4, POINTER :: q2(:,:)
2369 : REAL*4 :: minlon, maxlon
2370 0 : REAL*4 :: lon1s(im+1)
2371 :
2372 : ! Ghost correction
2373 : Logical :: isGlobal
2374 : Real*4 :: xSpan
2375 :
2376 : ! Missing value
2377 : REAL*8 :: miss
2378 :
2379 : ! Initialize
2380 0 : lon2 => NULL()
2381 0 : q2 => NULL()
2382 :
2383 : ! Missing value
2384 0 : miss = miss_r8
2385 0 : if ( present(missval) ) miss = missval
2386 :
2387 : ! XMAP begins here!
2388 0 : do i=1,im+1
2389 0 : x1(i) = lon1(i)
2390 : enddo
2391 :
2392 0 : do i=1,im
2393 0 : dx1(i) = x1(i+1) - x1(i)
2394 : enddo
2395 :
2396 : !===================================================================
2397 : ! define minimum and maximum longitude on output grid
2398 : ! to be used. Remapping will be restricted to this
2399 : ! domain. This procedure allows remapping of nested
2400 : ! domains onto larger (e.g. global) domains.
2401 : ! ckeller, 2/11/15).
2402 : !===================================================================
2403 0 : minlon = minval(lon1)
2404 0 : maxlon = maxval(lon1)
2405 :
2406 : ! check for values > 180.0
2407 0 : if(maxlon > 180.0) then
2408 0 : lon1s = lon1
2409 0 : do while(maxlon > 180.0)
2410 0 : WHERE(lon1s > 180.0) lon1s = lon1s - 360.0
2411 0 : minlon = minval(lon1s)
2412 0 : maxlon = maxval(lon1s)
2413 : enddo
2414 : endif
2415 :
2416 : ! maxlon must represent the easter edge of the grid:
2417 0 : maxlon = maxlon + ( lon1(im+1)-lon1(im) )
2418 :
2419 : ! Reduce input grid
2420 0 : n1 = 1
2421 0 : n2 = iin+1
2422 0 : do i=1,iin+1
2423 0 : if ( ilon2(i) < minlon ) n1 = i
2424 0 : if ( ilon2(iin+2-i) > maxlon ) n2 = iin+2-i
2425 : enddo
2426 0 : in = n2 - n1
2427 0 : lon2 => ilon2(n1:n2)
2428 0 : q2 => iq2(n1:(n2-1),:)
2429 :
2430 : ! Periodic BC only valid if the variable is "global"
2431 0 : xSpan = x1(im+1)-x1(1)
2432 0 : isGlobal = ((xSpan.ge.355.0).and.(xSpan.le.365.0))
2433 :
2434 : !===================================================================
2435 : ! check to see if ghosting is necessary
2436 : ! Western edge:
2437 : !===================================================================
2438 0 : found = .false.
2439 0 : i1 = 1
2440 : do while ( .not. found )
2441 0 : if( lon2(1) .ge. x1(i1) ) then
2442 : found = .true.
2443 : else
2444 0 : i1 = i1 - 1
2445 0 : if (i1 .lt. -im) then
2446 0 : write(6,*) 'Failed in Xmap_R4R8 (regrid_a2a_mod.F90)'
2447 0 : stop
2448 : else
2449 0 : x1(i1) = x1(i1+1) - dx1(im+i1)
2450 0 : dx1(i1) = dx1(im+i1)
2451 : endif
2452 : endif
2453 : enddo
2454 :
2455 : !===================================================================
2456 : ! Eastern edge:
2457 : !===================================================================
2458 : found = .false.
2459 : i2 = im+1
2460 : do while ( .not. found )
2461 0 : if( lon2(in+1) .le. x1(i2) ) then
2462 : found = .true.
2463 : else
2464 0 : i2 = i2 + 1
2465 0 : if (i2 .gt. 2*im) then
2466 0 : write(6,*) 'Failed in Xmap_R4R8 (regrid_a2a_mod.F90)'
2467 0 : stop
2468 : else
2469 0 : dx1(i2-1) = dx1(i2-1-im)
2470 0 : x1(i2) = x1(i2-1) + dx1(i2-1)
2471 : endif
2472 : endif
2473 : enddo
2474 :
2475 : !$OMP PARALLEL DO &
2476 : !$OMP DEFAULT( SHARED ) &
2477 : !$OMP PRIVATE( J, QTMP, I, I0, M, QSUM, DLON, MM, DX )
2478 0 : do 1000 j=1,jm
2479 :
2480 : !=================================================================
2481 : ! Area preserving mapping
2482 : !================================================================
2483 :
2484 0 : qtmp(:) = 0.0
2485 0 : do i=1,im
2486 0 : qtmp(i)=q1(i,j)
2487 : enddo
2488 :
2489 : ! SDE 2017-01-07
2490 : ! Only have shadow regions if we are on a global grid. Otherwise, we
2491 : ! should keep the zero boundary conditions.
2492 0 : If (isGlobal) Then
2493 0 : qtmp(0)=q1(im,j)
2494 0 : qtmp(im+1)=q1(1,j)
2495 :
2496 : ! check to see if ghosting is necessary
2497 : ! Western edge
2498 0 : if ( i1 .le. 0 ) then
2499 0 : do i=i1,0
2500 0 : qtmp(i) = qtmp(im+i)
2501 : enddo
2502 : endif
2503 :
2504 : ! Eastern edge:
2505 0 : if ( i2 .gt. im+1 ) then
2506 0 : do i=im+1,i2-1
2507 0 : qtmp(i) = qtmp(i-im)
2508 : enddo
2509 : endif
2510 : End If
2511 :
2512 : i0 = i1
2513 :
2514 0 : do 555 i=1,in
2515 0 : do 100 m=i0,i2-1
2516 :
2517 : !=============================================================
2518 : ! locate the western edge: lon2(i)
2519 : !=============================================================
2520 0 : if(lon2(i) .ge. x1(m) .and. lon2(i) .le. x1(m+1)) then
2521 :
2522 0 : if(lon2(i+1) .le. x1(m+1)) then
2523 :
2524 : ! entire new grid is within the original grid
2525 0 : if( abs(qtmp(m)-miss)>tiny_r8 ) q2(i,j)=qtmp(m)
2526 : i0 = m
2527 : goto 555
2528 : else
2529 :
2530 : ! Left most fractional area
2531 0 : if( abs(qtmp(m)-miss)>tiny_r8 ) then
2532 0 : qsum=(x1(m+1)-lon2(i))*qtmp(m)
2533 : dlon= x1(m+1)-lon2(i)
2534 : else
2535 : qsum=0.0
2536 : dlon=0.0
2537 : endif
2538 0 : do mm=m+1,i2-1
2539 :
2540 : ! locate the eastern edge: lon2(i+1)
2541 0 : if(lon2(i+1) .gt. x1(mm+1) ) then
2542 :
2543 : ! Whole layer
2544 0 : if( abs(qtmp(mm)-miss)>tiny_r8 ) then
2545 0 : qsum = qsum + dx1(mm)*qtmp(mm)
2546 0 : dlon = dlon + dx1(mm)
2547 : endif
2548 : else
2549 : ! Right most fractional area
2550 0 : if( abs(qtmp(m)-miss)>tiny_r8 ) then
2551 0 : dx = lon2(i+1)-x1(mm)
2552 0 : qsum=qsum+dx*qtmp(mm)
2553 0 : dlon=dlon+dx
2554 : endif
2555 : i0 = mm
2556 : goto 123
2557 : endif
2558 : enddo
2559 : goto 123
2560 : endif
2561 : endif
2562 0 : 100 continue
2563 0 : 123 if( ABS( dlon ) > 0.0e0 ) q2(i,j) = qsum / dlon
2564 0 : 555 continue
2565 0 : 1000 continue
2566 : !$OMP END PARALLEL DO
2567 :
2568 : ! Cleanup
2569 0 : lon2 => NULL()
2570 0 : q2 => NULL()
2571 :
2572 0 : END SUBROUTINE xmap_r8r4
2573 : !EOC
2574 : !------------------------------------------------------------------------------
2575 : ! Harmonized Emissions Component (HEMCO) !
2576 : !------------------------------------------------------------------------------
2577 : !BOP
2578 : !
2579 : ! !IROUTINE: Init_Map_A2A
2580 : !
2581 : ! !DESCRIPTION: Subroutine Init\_Map\_A2A initializes all module variables.
2582 : ! This allows us to keep "shadow" copies of variables that are defined
2583 : ! elsewhere in GEOS-Chem. This also helps us from having dependencies to
2584 : ! GEOS-Chem modules in the HEMCO core modules.
2585 : !\\
2586 : !\\
2587 : ! !INTERFACE:
2588 : !
2589 0 : SUBROUTINE Init_Map_A2A( NX, NY, LONS, SINES, AREAS, DIR )
2590 : !
2591 : ! !INPUT PARAMETERS:
2592 : !
2593 : INTEGER, INTENT(IN) :: NX ! # of longitudes
2594 : INTEGER, INTENT(IN) :: NY ! # of latitudes
2595 : REAL(fp), INTENT(IN) :: LONS (NX+1 ) ! Longitudes
2596 : REAL(fp), INTENT(IN) :: SINES(NY+1 ) ! Sines of latitudes
2597 : REAL(fp), INTENT(IN) :: AREAS(NX,NY) ! Surface areas [m2]
2598 : CHARACTER(LEN=*), INTENT(IN) :: DIR ! Dir for netCDF files w/
2599 : ! grid definitions
2600 : !
2601 : ! !REVISION HISTORY:
2602 : ! 14 Jul 2014 - R. Yantosca - Initial version
2603 : ! See https://github.com/geoschem/hemco for complete history
2604 : !EOP
2605 : !------------------------------------------------------------------------------
2606 : !BOC
2607 : !
2608 : ! !LOCAL VARIABLES:
2609 : !
2610 : INTEGER :: AS
2611 :
2612 : !------------------------------------------
2613 : ! Allocate module variables
2614 : !------------------------------------------
2615 0 : IF ( .not. ALLOCATED( OUTLON ) ) THEN
2616 0 : ALLOCATE( OUTLON( NX+1 ), STAT=AS )
2617 0 : IF ( AS /= 0 ) THEN
2618 0 : PRINT*, '### Could not allocate OUTLON (regrid_a2a_mod.F90)'
2619 0 : STOP
2620 : ENDIF
2621 : ENDIF
2622 :
2623 0 : IF ( .not. ALLOCATED( OUTSIN ) ) THEN
2624 0 : ALLOCATE( OUTSIN( NY+1 ), STAT=AS )
2625 0 : IF ( AS /= 0 ) THEN
2626 0 : PRINT*, '### Could not allocate OUTSIN (regrid_a2a_mod.F90)'
2627 0 : STOP
2628 : ENDIF
2629 : ENDIF
2630 :
2631 0 : IF ( .not. ALLOCATED( OUTAREA ) ) THEN
2632 0 : ALLOCATE( OUTAREA( NX, NY ), STAT=AS )
2633 0 : IF ( AS /= 0 ) THEN
2634 0 : PRINT*, '### Could not allocate OUTAREA (regrid_a2a_mod.F90)'
2635 0 : STOP
2636 : ENDIF
2637 : ENDIF
2638 :
2639 : !------------------------------------------
2640 : ! Store values in local shadow variables
2641 : !------------------------------------------
2642 0 : OUTNX = NX
2643 0 : OUTNY = NY
2644 0 : OUTLON = LONS
2645 0 : OUTSIN = SINES
2646 0 : OUTAREA = AREAS
2647 0 : NC_DIR = DIR
2648 :
2649 0 : END SUBROUTINE Init_Map_A2A
2650 : !EOC
2651 : !------------------------------------------------------------------------------
2652 : ! Harmonized Emissions Component (HEMCO) !
2653 : !------------------------------------------------------------------------------
2654 : !BOP
2655 : !
2656 : ! !IROUTINE: Cleanup_Map_A2A
2657 : !
2658 : ! !DESCRIPTION: Subroutine Cleanup\_Map\_A2A deallocates all module variables.
2659 : !\\
2660 : !\\
2661 : ! !INTERFACE:
2662 : !
2663 0 : SUBROUTINE Cleanup_Map_A2A()
2664 : !
2665 : ! !REVISION HISTORY:
2666 : ! 14 Jul 2014 - R. Yantosca - Initial version
2667 : ! See https://github.com/geoschem/hemco for complete history
2668 : !EOP
2669 : !------------------------------------------------------------------------------
2670 : !BOC
2671 : ! Cleanup module variables
2672 0 : IF ( ALLOCATED( OUTLON ) ) DEALLOCATE( OUTLON )
2673 0 : IF ( ALLOCATED( OUTSIN ) ) DEALLOCATE( OUTSIN )
2674 0 : IF ( ALLOCATED( OUTAREA ) ) DEALLOCATE( OUTAREA )
2675 :
2676 0 : END SUBROUTINE Cleanup_Map_A2A
2677 : !EOC
2678 : END MODULE HCO_Regrid_A2A_Mod
|
