Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hco_interp_mod.F90
7 : !
8 : ! !DESCRIPTION: Module HCO\_INTERP\_MOD contains routines to interpolate
9 : ! input data onto the HEMCO grid. This module contains routine for
10 : ! horizontal regridding between regular grids (MAP\_A2A), as well as
11 : ! vertical interpolation amongst GEOS model levels (full <--> reduced).
12 : !\\
13 : !\\
14 : ! Regridding is supported for concentration quantities (default) and
15 : ! index-based values. For the latter, the values in the regridded grid
16 : ! boxes correspond to the value of the original grid that contrbutes most
17 : ! to the given box.
18 : !\\
19 : !\\
20 : ! !INTERFACE:
21 : !
22 : MODULE HCO_Interp_Mod
23 : !
24 : ! !USES:
25 : !
26 : USE HCO_Types_Mod
27 : USE HCO_Error_Mod
28 : USE HCO_State_Mod, ONLY : Hco_State
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 : !
33 : ! !PUBLIC MEMBER FUNCTIONS:
34 : !
35 : PUBLIC :: ModelLev_Check
36 : PUBLIC :: ModelLev_Interpolate
37 : PUBLIC :: REGRID_MAPA2A
38 : !
39 : ! !PRIVATE MEMBER FUNCTIONS:
40 : !
41 : PRIVATE :: COLLAPSE
42 : !
43 : ! !REVISION HISTORY:
44 : ! 30 Dec 2014 - C. Keller - Initialization
45 : ! See https://github.com/geoschem/hemco for complete history
46 : !EOP
47 : !------------------------------------------------------------------------------
48 : !BOC
49 : !
50 : ! !PRIVATE VARIABLES:
51 : !
52 : ! AP parameter of native GEOS-5 grid. Needed to remap GEOS-5 data from native
53 : ! onto the reduced vertical grid.
54 : REAL(hp), TARGET :: G5_EDGE_NATIVE(73) = (/ &
55 : 0.000000e+00_hp, 4.804826e-02_hp, 6.593752e+00_hp, 1.313480e+01_hp, &
56 : 1.961311e+01_hp, 2.609201e+01_hp, 3.257081e+01_hp, 3.898201e+01_hp, &
57 : 4.533901e+01_hp, 5.169611e+01_hp, 5.805321e+01_hp, 6.436264e+01_hp, &
58 : 7.062198e+01_hp, 7.883422e+01_hp, 8.909992e+01_hp, 9.936521e+01_hp, &
59 : 1.091817e+02_hp, 1.189586e+02_hp, 1.286959e+02_hp, 1.429100e+02_hp, &
60 : 1.562600e+02_hp, 1.696090e+02_hp, 1.816190e+02_hp, 1.930970e+02_hp, &
61 : 2.032590e+02_hp, 2.121500e+02_hp, 2.187760e+02_hp, 2.238980e+02_hp, &
62 : 2.243630e+02_hp, 2.168650e+02_hp, 2.011920e+02_hp, 1.769300e+02_hp, &
63 : 1.503930e+02_hp, 1.278370e+02_hp, 1.086630e+02_hp, 9.236572e+01_hp, &
64 : 7.851231e+01_hp, 6.660341e+01_hp, 5.638791e+01_hp, 4.764391e+01_hp, &
65 : 4.017541e+01_hp, 3.381001e+01_hp, 2.836781e+01_hp, 2.373041e+01_hp, &
66 : 1.979160e+01_hp, 1.645710e+01_hp, 1.364340e+01_hp, 1.127690e+01_hp, &
67 : 9.292942e+00_hp, 7.619842e+00_hp, 6.216801e+00_hp, 5.046801e+00_hp, &
68 : 4.076571e+00_hp, 3.276431e+00_hp, 2.620211e+00_hp, 2.084970e+00_hp, &
69 : 1.650790e+00_hp, 1.300510e+00_hp, 1.019440e+00_hp, 7.951341e-01_hp, &
70 : 6.167791e-01_hp, 4.758061e-01_hp, 3.650411e-01_hp, 2.785261e-01_hp, &
71 : 2.113490e-01_hp, 1.594950e-01_hp, 1.197030e-01_hp, 8.934502e-02_hp, &
72 : 6.600001e-02_hp, 4.758501e-02_hp, 3.270000e-02_hp, 2.000000e-02_hp, &
73 : 1.000000e-02_hp /)
74 :
75 : ! AP parameter of native 102-layer GISS grid
76 : REAL(hp), TARGET :: E102_EDGE_NATIVE(103) = (/ &
77 : 0.0000000, 2.7871507, 5.5743014, 8.3614521, 11.1486028, 13.9357536, &
78 : 16.7229043, 19.5100550, 22.2972057, 25.0843564, 27.8715071, 30.6586578, &
79 : 33.4458085, 36.2329593, 39.0201100, 41.8087123, 44.6089278, 47.4534183, &
80 : 50.4082336, 53.5662786, 57.0095710, 60.7533531, 64.7323011, 68.8549615, &
81 : 73.0567364, 77.2969797, 81.5364973, 85.7346430, 89.8565776, 93.8754457, &
82 : 97.7709243, 101.5277712, 105.1350991, 108.5878272, 111.8859556, 115.0302100, &
83 : 118.0249453, 120.8854039, 123.6326345, 126.2811535, 128.8360417, 131.2987506, &
84 : 133.6736353, 135.9708571, 138.2013035, 140.3700552, 142.4814670, 144.5457005, &
85 : 146.5692881, 148.5464231, 150.4712991, 152.3497225, 154.1875000, 144.5468750, &
86 : 135.1875000, 126.0781250, 117.1914062, 108.5859375, 100.3671875, 92.5898438, &
87 : 85.2265625, 78.2226562, 71.5546875, 65.2226562, 59.2226562, 53.5546875, &
88 : 48.2226562, 43.2226562, 38.5546875, 34.2226562, 30.2226562, 26.5507812, &
89 : 23.1875000, 20.0781250, 17.1896562, 14.5684375, 12.2865742, 10.3573086, &
90 : 8.7353750, 7.3664922, 6.2100156, 5.2343633, 4.4119297, 3.7186797, &
91 : 3.1341479, 2.6404328, 2.2207877, 1.8587369, 1.5477125, 1.2782115, &
92 : 1.0427319, 0.8367716, 0.6514691, 0.4772511, 0.3168814, 0.1785988, &
93 : 0.1000000, 0.0560000, 0.0320000, 0.0180000, 0.0100000, 0.0050000, &
94 : 0.0020000 /)
95 :
96 : CONTAINS
97 : !EOC
98 : !------------------------------------------------------------------------------
99 : ! Harmonized Emissions Component (HEMCO) !
100 : !------------------------------------------------------------------------------
101 : !BOP
102 : !
103 : ! !IROUTINE: Regrid_MAPA2A
104 : !
105 : ! !DESCRIPTION: Subroutine Regrid\_MAPA2A regrids input array NcArr onto
106 : ! the simulation grid and stores the data in list container Lct. Horizontal
107 : ! regridding is performed using MAP\_A2A algorithm. Vertical interpolation
108 : ! between GEOS levels (full vs. reduced, GEOS-5 vs. GEOS-4), is also
109 : ! supported.
110 : !\\
111 : !\\
112 : ! This routine can remap concentrations and index-based quantities.
113 : !\\
114 : !\\
115 : ! !INTERFACE:
116 : !
117 0 : SUBROUTINE REGRID_MAPA2A( HcoState, NcArr, LonE, LatE, Lct, RC )
118 : !
119 : ! !USES:
120 : !
121 : USE HCO_REGRID_A2A_Mod, ONLY : MAP_A2A
122 : USE HCO_FileData_Mod, ONLY : FileData_ArrCheck
123 : USE HCO_UNIT_MOD, ONLY : HCO_IsIndexData
124 : !
125 : ! !INPUT PARAMETERS:
126 : !
127 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
128 : REAL(sp), POINTER :: NcArr(:,:,:,:) ! 4D input data
129 : REAL(hp), POINTER :: LonE(:) ! Input grid longitude edges
130 : REAL(hp), POINTER :: LatE(:) ! Input grid latitude edges
131 : !
132 : ! !INPUT/OUTPUT PARAMETERS:
133 : !
134 : TYPE(ListCont), POINTER :: Lct ! HEMCO list container
135 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
136 : !
137 : ! !REVISION HISTORY:
138 : ! 03 Feb 2015 - C. Keller - Initial version
139 : ! See https://github.com/geoschem/hemco for complete history
140 : !EOP
141 : !------------------------------------------------------------------------------
142 : !BOC
143 : !
144 : ! !LOCAL VARIABLES:
145 : !
146 : INTEGER :: nLonEdge, nLatEdge
147 : INTEGER :: NX, NY, NZ, NLEV, NTIME, NCELLS
148 : INTEGER :: I, J, L, T, AS, I2
149 : INTEGER :: nIndex
150 0 : REAL(sp), ALLOCATABLE :: LonEdgeI(:)
151 0 : REAL(sp), ALLOCATABLE :: LatEdgeI(:)
152 0 : REAL(sp) :: LonEdgeO(HcoState%NX+1)
153 0 : REAL(sp) :: LatEdgeO(HcoState%NY+1)
154 :
155 0 : REAL(sp), POINTER :: ORIG_2D(:,:)
156 0 : REAL(sp), POINTER :: REGR_2D(:,:)
157 0 : REAL(sp), POINTER :: REGR_4D(:,:,:,:)
158 :
159 0 : REAL(sp), ALLOCATABLE, TARGET :: FRACS(:,:,:,:)
160 0 : REAL(hp), ALLOCATABLE :: REGFRACS(:,:,:,:)
161 0 : REAL(hp), ALLOCATABLE :: MAXFRACS(:,:,:,:)
162 0 : REAL(hp), ALLOCATABLE :: INDECES(:,:,:,:)
163 0 : REAL(hp), ALLOCATABLE :: UNIQVALS(:)
164 : REAL(hp) :: IVAL
165 : LOGICAL :: IsIndex
166 :
167 : LOGICAL :: VERB
168 : CHARACTER(LEN=255) :: MSG
169 : CHARACTER(LEN=255) :: LOC = 'ModelLev_Interpolate (hco_interp_mod.F90)'
170 :
171 : !=================================================================
172 : ! REGRID_MAPA2A begins here
173 : !=================================================================
174 :
175 : ! Init
176 0 : ORIG_2D => NULL()
177 0 : REGR_2D => NULL()
178 0 : REGR_4D => NULL()
179 :
180 : ! Check for verbose mode
181 0 : verb = HCO_IsVerb( HcoState%Config%Err )
182 :
183 : ! get longitude / latitude sizes
184 0 : nLonEdge = SIZE(LonE,1)
185 0 : nLatEdge = SIZE(LatE,1)
186 :
187 : ! Write input grid edges to shadow variables so that map_a2a accepts them
188 : ! as argument.
189 : ! Also, for map_a2a, latitudes have to be sines...
190 0 : ALLOCATE(LonEdgeI(nlonEdge), LatEdgeI(nlatEdge), STAT=AS )
191 0 : IF ( AS /= 0 ) THEN
192 0 : CALL HCO_ERROR( 'alloc error LonEdgeI/LatEdgeI', RC, THISLOC=LOC )
193 0 : RETURN
194 : ENDIF
195 0 : LonEdgeI(:) = LonE
196 0 : LatEdgeI(:) = SIN( LatE * HcoState%Phys%PI_180 )
197 :
198 : ! Get output grid edges from HEMCO state
199 0 : LonEdgeO(:) = HcoState%Grid%XEDGE%Val(:,1)
200 0 : LatEdgeO(:) = HcoState%Grid%YSIN%Val(1,:)
201 :
202 : ! Get input array sizes
203 0 : NX = size(ncArr,1)
204 0 : NY = size(ncArr,2)
205 0 : NLEV = size(ncArr,3)
206 0 : NTIME = size(ncArr,4)
207 0 : NCELLS = NX * NY * NLEV * NTIME
208 :
209 : ! Are these index-based data? If so, need to remap the fraction (1 or 0)
210 : ! of every value independently. For every grid box, the value with the
211 : ! highest overlap (closest to 1) is taken.
212 0 : IsIndex = HCO_IsIndexData(Lct%Dct%Dta%OrigUnit)
213 :
214 0 : IF ( IsIndex ) THEN
215 :
216 : ! Allocate working arrays:
217 : ! - FRACS contains the fractions on the original grid. These are
218 : ! binary (1 or 0).
219 : ! - MAXFRACS stores the highest used fraction for each output grid
220 : ! box. Will be updated continously.
221 : ! - INDECES is the output array holding the index-based remapped
222 : ! values. Will be updated continuously.
223 : ! - UNIQVALS is a vector holding all unique values of the input
224 : ! array (NINDEX is the number of unique values).
225 : !
226 : ! ckeller, 9/24/15: Extend vertical axis of MAXFRACS, REGFRACS, and
227 : ! INDECES to HcoState%NZ+1 for fields that are on edges instead of
228 : ! mid-points.
229 0 : ALLOCATE( FRACS(NX,NY,NLEV,NTIME), STAT=AS )
230 0 : IF ( AS /= 0 ) THEN
231 0 : CALL HCO_ERROR( 'alloc error FRACS', RC, THISLOC=LOC )
232 0 : RETURN
233 : ENDIF
234 0 : ALLOCATE( MAXFRACS(HcoState%NX,HcoState%NY,HcoState%NZ+1,NTIME), STAT=AS )
235 0 : IF ( AS /= 0 ) THEN
236 0 : CALL HCO_ERROR( 'alloc error MAXFRACS', RC, THISLOC=LOC )
237 0 : RETURN
238 : ENDIF
239 0 : ALLOCATE( REGFRACS(HcoState%NX,HcoState%NY,HcoState%NZ+1,NTIME), STAT=AS )
240 0 : IF ( AS /= 0 ) THEN
241 0 : CALL HCO_ERROR( 'alloc error INDECES', RC, THISLOC=LOC )
242 0 : RETURN
243 : ENDIF
244 0 : ALLOCATE( INDECES(HcoState%NX,HcoState%NY,HcoState%NZ+1,NTIME), STAT=AS )
245 0 : IF ( AS /= 0 ) THEN
246 0 : CALL HCO_ERROR( 'alloc error INDECES', RC, THISLOC=LOC )
247 0 : RETURN
248 : ENDIF
249 0 : ALLOCATE( UNIQVALS(NCELLS), STAT=AS )
250 0 : IF ( AS /= 0 ) THEN
251 0 : CALL HCO_ERROR( 'alloc error INDECES', RC, THISLOC=LOC )
252 0 : RETURN
253 : ENDIF
254 0 : FRACS = 0.0_sp
255 0 : REGFRACS = 0.0_hp
256 0 : MAXFRACS = 0.0_hp
257 0 : INDECES = 0.0_hp
258 0 : UNIQVALS = 0.0_hp
259 :
260 : ! Get unique values. Loop over all input data values and add
261 : ! them to UNIQVALS vector if UNIQVALS doesn't hold that same value
262 : ! yet.
263 0 : NINDEX = 0
264 0 : DO T = 1, NTIME
265 0 : DO L = 1, NLEV
266 0 : DO J = 1, NY
267 0 : DO I = 1, NX
268 :
269 : ! Current value
270 0 : IVAL = NcArr(I,J,L,T)
271 :
272 : ! Check if value already exists in UNIQVALS
273 0 : IF ( NINDEX > 0 ) THEN
274 0 : IF ( ANY(UNIQVALS(1:NINDEX) == IVAL) ) CYCLE
275 : ENDIF
276 :
277 : ! Add to UNIQVALS
278 0 : NINDEX = NINDEX + 1
279 0 : UNIQVALS(NINDEX) = IVAL
280 : ENDDO
281 : ENDDO
282 : ENDDO
283 : ENDDO
284 :
285 : ! Verbose mode
286 0 : IF ( verb ) THEN
287 0 : MSG = 'Do index based regridding for field ' // TRIM(Lct%Dct%cName)
288 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
289 0 : WRITE(MSG,*) ' - Number of indeces: ', NINDEX
290 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
291 : ENDIF
292 :
293 : ELSE
294 0 : NINDEX = 1
295 : ENDIF
296 :
297 : ! Define array to put horizontally regridded data onto. If this
298 : ! is 3D data, we first regrid all vertical levels horizontally
299 : ! and then pass these data to the list container. In this second
300 : ! step, levels may be deflated/collapsed.
301 :
302 : ! 2D data is directly passed to the data container
303 0 : IF ( Lct%Dct%Dta%SpaceDim <= 2 ) THEN
304 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, &
305 0 : HcoState%NX, HcoState%NY, NTIME, RC )
306 0 : IF ( RC /= 0 ) RETURN
307 : ENDIF
308 :
309 : ! 3D data and index data is first written into a temporary array,
310 : ! REGR_4D.
311 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 .OR. IsIndex ) THEN
312 0 : ALLOCATE( REGR_4D(HcoState%NX,HcoState%NY,NLEV,NTIME), STAT=AS )
313 : IF ( AS /= 0 ) THEN
314 0 : CALL HCO_ERROR( 'alloc error REGR_4D', RC, THISLOC=LOC )
315 0 : RETURN
316 : ENDIF
317 0 : REGR_4D = 0.0_hp
318 : ENDIF
319 :
320 : ! Do regridding for every index value. If it's not index data, this loop
321 : ! is executed only once (NINDEX=1).
322 0 : DO I = 1, NINDEX
323 :
324 : ! For index based data, create fractions array for the given index.
325 0 : IF ( IsIndex ) THEN
326 0 : IVAL = UNIQVALS(I)
327 0 : WHERE( ncArr == IVAL )
328 : FRACS = 1.0_sp
329 : ELSEWHERE
330 : FRACS = 0.0_sp
331 : END WHERE
332 : ENDIF
333 :
334 : ! Regrid horizontally
335 0 : DO T = 1, NTIME
336 0 : DO L = 1, NLEV
337 :
338 : ! Point to 2D slices to be regridded:
339 : ! - Original 2D array
340 0 : IF ( IsIndex ) THEN
341 0 : ORIG_2D => FRACS(:,:,L,T)
342 : ELSE
343 0 : ORIG_2D => ncArr(:,:,L,T)
344 : ENDIF
345 :
346 : ! - Regridded 2D array
347 0 : IF ( Lct%Dct%Dta%SpaceDim <= 2 .AND. .NOT. IsIndex ) THEN
348 0 : REGR_2D => Lct%Dct%Dta%V2(T)%Val(:,:)
349 : ELSE
350 0 : REGR_2D => REGR_4D(:,:,L,T)
351 : ENDIF
352 :
353 : ! Do the regridding
354 : CALL MAP_A2A( NX, NY, LonEdgeI, LatEdgeI, ORIG_2D, &
355 : HcoState%NX, HcoState%NY, LonEdgeO, LatEdgeO, &
356 0 : REGR_2D, 0, 0, HCO_MISSVAL )
357 0 : ORIG_2D => NULL()
358 0 : REGR_2D => NULL()
359 :
360 : ENDDO !L
361 : ENDDO !T
362 :
363 : ! Eventually collapse levels onto simulation levels.
364 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
365 0 : CALL ModelLev_Interpolate( HcoState, REGR_4D, Lct, RC )
366 0 : IF ( RC /= HCO_SUCCESS ) THEN
367 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
368 0 : RETURN
369 : ENDIF
370 : ENDIF
371 :
372 : ! For index based data, map fractions back to corresponding value.
373 : ! Array INDECES holds the index-based remapped values. Set INDECES
374 : ! to current index value in every grid box where the regridded
375 : ! fraction of this index is higher than any previous fraction
376 : ! (array MAXFRACS stores the highest used fraction in each grid box).
377 0 : IF ( IsIndex ) THEN
378 :
379 : ! Reset
380 0 : REGFRACS = 0.0_hp
381 :
382 : ! 3D data written to Lct needs to be mapped back onto REGR_4D.
383 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
384 0 : DO T = 1, NTIME
385 0 : NZ = SIZE(Lct%Dct%Dta%V3(T)%Val,3)
386 0 : REGFRACS(:,:,1:NZ,T) = Lct%Dct%Dta%V3(T)%Val(:,:,:)
387 : ENDDO
388 : ELSE
389 0 : REGFRACS(:,:,1:NLEV,:) = REGR_4D(:,:,:,:)
390 : ENDIF
391 :
392 : ! REGR_4D are the remapped fractions.
393 0 : DO T = 1, NTIME
394 0 : DO L = 1, HcoState%NZ
395 0 : DO J = 1, HcoState%NY
396 0 : DO I2 = 1, HcoState%NX
397 0 : IF ( REGFRACS(I2,J,L,T) > MAXFRACS(I2,J,L,T) ) THEN
398 0 : MAXFRACS(I2,J,L,T) = REGR_4D(I2,J,L,T)
399 0 : INDECES (I2,J,L,T) = IVAL
400 : ENDIF
401 : ENDDO
402 : ENDDO
403 : ENDDO
404 : ENDDO
405 :
406 : !------------------------------------------------------------------------------
407 : ! Prior to 9/29/16:
408 : ! ! This code is preblematic in Gfortran. Replace it with the
409 : ! ! explicit DO loops above. Leave this here for reference.
410 : ! ! (sde, bmy, 9/21/16)
411 : ! WHERE ( REGFRACS > MAXFRACS )
412 : ! MAXFRACS = REGR_4D
413 : ! INDECES = IVAL
414 : ! END WHERE
415 : !------------------------------------------------------------------------------
416 : ENDIF
417 :
418 : ENDDO !I
419 :
420 : ! For index values, pass index data to data container.
421 0 : IF ( IsIndex ) THEN
422 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
423 0 : DO T = 1, NTIME
424 0 : NZ = SIZE(Lct%Dct%Dta%V3(T)%Val,3)
425 0 : Lct%Dct%Dta%V3(T)%Val(:,:,:) = INDECES(:,:,1:NZ,T)
426 : ENDDO
427 : ELSE
428 0 : DO T = 1, NTIME
429 0 : Lct%Dct%Dta%V2(T)%Val(:,:) = INDECES(:,:,1,T)
430 : ENDDO
431 : ENDIF
432 : ENDIF
433 :
434 : ! Cleanup
435 0 : DEALLOCATE(LonEdgeI, LatEdgeI)
436 0 : IF ( ASSOCIATED( REGR_4D ) ) DEALLOCATE( REGR_4D )
437 0 : IF ( ALLOCATED ( FRACS ) ) DEALLOCATE( FRACS )
438 0 : IF ( ALLOCATED ( REGFRACS ) ) DEALLOCATE( REGFRACS )
439 0 : IF ( ALLOCATED ( MAXFRACS ) ) DEALLOCATE( MAXFRACS )
440 0 : IF ( ALLOCATED ( INDECES ) ) DEALLOCATE( INDECES )
441 0 : IF ( ALLOCATED ( UNIQVALS ) ) DEALLOCATE( UNIQVALS )
442 :
443 : ! Return w/ success
444 0 : RC = HCO_SUCCESS
445 :
446 0 : END SUBROUTINE REGRID_MAPA2A
447 : !EOC
448 : !------------------------------------------------------------------------------
449 : ! Harmonized Emissions Component (HEMCO) !
450 : !------------------------------------------------------------------------------
451 : !BOP
452 : !
453 : ! !IROUTINE: ModelLev_Check
454 : !
455 : ! !DESCRIPTION: Subroutine ModelLev\_Check checks if the passed number of
456 : ! vertical levels indicates that these are model levels or not.
457 : !\\
458 : !\\
459 : ! !INTERFACE:
460 : !
461 0 : SUBROUTINE ModelLev_Check( HcoState, nLev, IsModelLev, RC )
462 : !
463 : ! !USES:
464 : !
465 : USE HCO_FileData_Mod, ONLY : FileData_ArrCheck
466 : !
467 : ! !INPUT PARAMETERS:
468 : !
469 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
470 : INTEGER, INTENT(IN ) :: nlev ! number of levels
471 : !
472 : ! !INPUT/OUTPUT PARAMETERS:
473 : !
474 : LOGICAL, INTENT(INOUT) :: IsModelLev ! Are these model levels?
475 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
476 : !
477 : ! !REVISION HISTORY:
478 : ! 29 Sep 2015 - C. Keller - Initial version
479 : ! See https://github.com/geoschem/hemco for complete history
480 : !EOP
481 : !------------------------------------------------------------------------------
482 : !BOC
483 : !
484 : ! !LOCAL VARIABLES:
485 : !
486 : INTEGER :: nz
487 :
488 : !=================================================================
489 : ! ModelLev_Check begins here
490 : !=================================================================
491 :
492 : ! Assume success until otherwise
493 0 : RC = HCO_SUCCESS
494 :
495 : ! Shadow number of vertical levels on grid
496 0 : nz = HcoState%NZ
497 :
498 : ! Assume model levels if input data levels correspond to # of grid
499 : ! levels or levels + 1 (edges)
500 0 : IF ( nlev == nz .OR. nlev == nz + 1 ) THEN
501 0 : IsModelLev = .TRUE.
502 :
503 : ! If input is 72 layer (or 3/11/36 layer) and output is 47 layer
504 0 : ELSEIF ( nz == 47 ) THEN
505 0 : IsModelLev = ( nlev == 72 .OR. nlev == 73 .OR. nlev == 36 .OR. nlev == 3 .OR. nlev == 11)
506 :
507 : ! If input is 102 layer and output is 74 layer
508 0 : ELSEIF ( nz == 74 ) THEN
509 0 : IsModelLev = ( nlev == 102 .OR. nlev == 103 )
510 :
511 : ! If input is 47 layer (or 3/11/36 layer) and output is 72 layer
512 0 : ELSEIF ( nz == 72 ) THEN
513 0 : IsModelLev = ( nlev == 47 .OR. nlev == 48 .OR. nlev == 36 .OR. nlev == 3 .OR. nlev == 11)
514 :
515 : ELSE
516 0 : IsModelLev = .FALSE.
517 : ENDIF
518 :
519 0 : END SUBROUTINE ModelLev_Check
520 : !EOC
521 : !------------------------------------------------------------------------------
522 : ! Harmonized Emissions Component (HEMCO) !
523 : !------------------------------------------------------------------------------
524 : !BOP
525 : !
526 : ! !IROUTINE: ModelLev_Interpolate
527 : !
528 : ! !DESCRIPTION: Subroutine ModelLev\_Interpolate puts 3D data from an
529 : ! arbitrary number of model levels onto the vertical levels of the simulation
530 : ! grid. Since the input data is already on model levels, this is only to
531 : ! collapse fields between native/reduced vertical levels, e.g. from
532 : ! 72 native GEOS-5 levels onto the reduced 47 levels.
533 : !
534 : !
535 : ! The input data (REGR\_4D) is expected to be already regridded horizontally.
536 : ! The 4th dimension of REGR\_4D denotes time.
537 : !
538 : !
539 : ! The 3rd dimension of REGR\_3D holds the vertical levels. It is assumed that
540 : ! these are model levels, starting at the surface (level 1). If the input
541 : ! data holds 72/73 input levels, this is interpreted as native data and will
542 : ! be collapsed onto the reduced GEOS-5 grid. If the input holds 102/103 input
543 : ! levels, this is interpreted as native data and will be collapsed onto the
544 : ! reduced GISS grid. If the input holds 47/48 input levels, this is interpreted
545 : ! as reduced GEOS-5 data and it will be inflated to the native GEOS-5 grid
546 : ! (with a warning, as this is not recommended). If the input holds N input levels,
547 : ! (where N = 3, 11, or 36 to account for NEI and AEIC emissions), this is assumed
548 : ! to be the first N levels of the GEOS-5 grid, meaning they will be written as
549 : ! levels 1-N of a 47 or 72 level output grid (with the remaining values left to
550 : ! be zero) (nbalasus, 8/29/2023).
551 : !
552 : !
553 : ! Currently, this routine can remap the following combinations:
554 : !
555 : ! * Native GEOS-5 onto reduced GEOS-5 (72 --> 47 levels, 73 --> 48 edges)
556 : ! * Native GISS onto reduced GISS (102 --> 74 levels, 103 --> 75 edges)
557 : ! * Reduced GEOS-5 onto native GEOS-5 (47 --> 72 levels, 48 --> 73 edges)
558 : ! * N (N = 3/11/36) levels onto native/reduced GEOS-5 (N --> levels 1-N levels of 47/72 level grid, rest are 0)
559 : !
560 : ! !INTERFACE:
561 : !
562 0 : SUBROUTINE ModelLev_Interpolate( HcoState, REGR_4D, Lct, RC )
563 : !
564 : ! !USES:
565 : !
566 : USE HCO_FileData_Mod, ONLY : FileData_ArrCheck
567 : !
568 : ! !INPUT PARAMETERS:
569 : !
570 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
571 : REAL(sp), POINTER :: REGR_4D(:,:,:,:) ! 4D input data
572 : !
573 : ! !INPUT/OUTPUT PARAMETERS:
574 : !
575 : TYPE(ListCont), POINTER :: Lct ! HEMCO list container
576 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
577 : !
578 : ! !REVISION HISTORY:
579 : ! 30 Dec 2014 - C. Keller - Initial version
580 : ! See https://github.com/geoschem/hemco for complete history
581 : !EOP
582 : !------------------------------------------------------------------------------
583 : !BOC
584 : !
585 : ! !LOCAL VARIABLES:
586 : !
587 : INTEGER :: nx, ny, nz, nt
588 : INTEGER :: fineIDX, coarseIDX
589 : INTEGER :: minlev, nlev, nout
590 : INTEGER :: L, T, NL, I
591 : INTEGER :: OS
592 : LOGICAL :: verb, infl, clps
593 : LOGICAL :: DONE
594 : CHARACTER(LEN=255) :: MSG, LOC
595 :
596 : !=================================================================
597 : ! ModelLev_Interpolate begins here
598 : !=================================================================
599 0 : LOC = 'ModelLev_Interpolate (HCO_INTERP_MOD.F90)'
600 :
601 : ! Enter
602 0 : CALL HCO_ENTER (HcoState%Config%Err, LOC, RC )
603 0 : IF ( RC /= HCO_SUCCESS ) THEN
604 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
605 0 : RETURN
606 : ENDIF
607 :
608 : ! Check for verbose mode
609 0 : verb = HCO_IsVerb( HcoState%Config%Err )
610 0 : IF ( verb ) THEN
611 0 : MSG = 'Vertically interpolate model levels: '//TRIM(Lct%Dct%cName)
612 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
613 : ENDIF
614 :
615 : ! Get HEMCO grid dimensions
616 0 : nx = HcoState%NX
617 0 : ny = HcoState%NY
618 0 : nz = HcoState%NZ
619 :
620 : ! Input data must be on horizontal HEMCO grid
621 0 : IF ( SIZE(REGR_4D,1) /= nx ) THEN
622 0 : WRITE(MSG,*) 'x dimension mismatch ', TRIM(Lct%Dct%cName), &
623 0 : ': ', nx, SIZE(REGR_4D,1)
624 0 : CALL HCO_ERROR( MSG, RC )
625 0 : RETURN
626 : ENDIF
627 0 : IF ( SIZE(REGR_4D,2) /= ny ) THEN
628 0 : WRITE(MSG,*) 'y dimension mismatch ', TRIM(Lct%Dct%cName), &
629 0 : ': ', ny, SIZE(REGR_4D,2)
630 0 : CALL HCO_ERROR( MSG, RC )
631 0 : RETURN
632 : ENDIF
633 :
634 : ! Get vertical and time dimension of input data
635 0 : nlev = SIZE(REGR_4D,3)
636 0 : nt = SIZE(REGR_4D,4)
637 :
638 : ! Check to make sure ModelLev_Interpolate should have been called
639 : IF ( ( ( nlev == nz ) .OR. ( nlev == nz+1 ) ) .OR. & ! write data without doing anything
640 : ( ( nz == 47 ) .AND. ( ( nlev == 72 ) .OR. ( nlev == 73 ) ) ) .OR. & ! collapse native to reduced GEOS-5
641 : ( ( nz == 74 ) .AND. ( ( nlev == 102 ) .OR. ( nlev == 103 ) ) ) .OR. & ! collapse native to reduced GISS
642 0 : ( ( nz == 72 ) .AND. ( ( nlev == 47 ) .OR. ( nlev == 48 ) ) ) .OR. & ! inflate reduced to native GEOS-5
643 : ( ( ( nz == 72 ) .OR. ( nz == 47 ) ) .AND. ( ( nlev == 3 ) .OR. ( nlev == 11 ) .OR. ( nlev == 36 ) ) ) ) THEN ! write 3/11/36 levels to reduced/native GEOS-5
644 : ! do nothing
645 : ELSE
646 0 : WRITE(MSG,*) 'ModelLev_Interpolate was called but MESSy should have been used: ',TRIM(Lct%Dct%cName)
647 0 : CALL HCO_ERROR( MSG, RC )
648 0 : RETURN
649 : ENDIF
650 :
651 : ! Vertical interpolation done?
652 0 : DONE = .FALSE.
653 :
654 : !===================================================================
655 : ! If no vertical interpolation is needed, then (1) save the 4D
656 : ! input data array to to the HEMCO list container object and
657 : ! (2) exit this subroutine.
658 : !===================================================================
659 0 : IF ( ( nlev == nz ) .OR. ( nlev == nz+1 ) ) THEN
660 :
661 0 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nlev, nt, RC )
662 0 : IF ( RC /= HCO_SUCCESS ) THEN
663 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
664 0 : RETURN
665 : ENDIF
666 :
667 0 : DO T = 1, nt
668 0 : Lct%Dct%Dta%V3(T)%Val(:,:,:) = REGR_4D(:,:,:,T)
669 : ENDDO
670 :
671 : ! Verbose
672 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
673 0 : MSG = '# of input levels = # of output levels - passed as is.'
674 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
675 : ENDIF
676 :
677 : ! Done!
678 : DONE = .TRUE.
679 : ENDIF
680 :
681 : !===================================================================
682 : ! Do vertical regridding:
683 : !===================================================================
684 : IF ( .NOT. DONE ) THEN
685 :
686 : !----------------------------------------------------------------
687 : ! Native to reduced GEOS-5 levels
688 : !----------------------------------------------------------------
689 0 : IF ( nz == 47 ) THEN
690 :
691 : ! Determine if the variable is on model levels or edges
692 : IF ( nlev == 72 ) THEN
693 0 : NL = 36
694 0 : nout = 47
695 : ELSEIF ( nlev == 73 ) THEN
696 0 : NL = 37
697 0 : nout = 48
698 : ELSEIF ( nlev == 3 ) THEN
699 0 : NL = 3
700 0 : nout = 47
701 : ELSEIF ( nlev == 11 ) THEN
702 0 : NL = 11
703 0 : nout = 47
704 : ELSEIF ( nlev == 36 ) THEN
705 0 : NL = 36
706 0 : nout = 47
707 : ELSE
708 : MSG = 'Can only remap from native onto reduced GEOS-5 if '// &
709 0 : 'input data has exactly 72, 73, 3, 11, or 36 levels: '//TRIM(Lct%Dct%cName)
710 0 : CALL HCO_ERROR( MSG, RC )
711 0 : RETURN
712 : ENDIF ! nlev == (72,73,3,11,36ELSE)
713 :
714 : ! Make sure output array is allocated
715 0 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )
716 :
717 : ! Do for every time slice
718 0 : DO T = 1, nt
719 :
720 : ! Levels that are passed level-by-level.
721 0 : DO L = 1, NL
722 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
723 : ENDDO !L
724 :
725 : ! If remapping model grid layers, collapse layers
726 0 : IF ( nlev == 72 ) THEN
727 : ! Collapse two levels (e.g. levels 37-38 into level 37):
728 0 : CALL COLLAPSE( Lct, REGR_4D, 37, 37, 2, T, 5, RC )
729 0 : CALL COLLAPSE( Lct, REGR_4D, 38, 39, 2, T, 5, RC )
730 0 : CALL COLLAPSE( Lct, REGR_4D, 39, 41, 2, T, 5, RC )
731 0 : CALL COLLAPSE( Lct, REGR_4D, 40, 43, 2, T, 5, RC )
732 : ! Collapse four levels:
733 0 : CALL COLLAPSE( Lct, REGR_4D, 41, 45, 4, T, 5, RC )
734 0 : CALL COLLAPSE( Lct, REGR_4D, 42, 49, 4, T, 5, RC )
735 0 : CALL COLLAPSE( Lct, REGR_4D, 43, 53, 4, T, 5, RC )
736 0 : CALL COLLAPSE( Lct, REGR_4D, 44, 57, 4, T, 5, RC )
737 0 : CALL COLLAPSE( Lct, REGR_4D, 45, 61, 4, T, 5, RC )
738 0 : CALL COLLAPSE( Lct, REGR_4D, 46, 65, 4, T, 5, RC )
739 0 : CALL COLLAPSE( Lct, REGR_4D, 47, 69, 4, T, 5, RC )
740 : ! If remapping model grid edges, sample at edges
741 0 : ELSEIF ( nlev == 73 ) THEN
742 0 : Lct%Dct%Dta%V3(T)%Val(:,:,38) = REGR_4D(:,:,39,T)
743 0 : Lct%Dct%Dta%V3(T)%Val(:,:,39) = REGR_4D(:,:,41,T)
744 0 : Lct%Dct%Dta%V3(T)%Val(:,:,40) = REGR_4D(:,:,43,T)
745 0 : Lct%Dct%Dta%V3(T)%Val(:,:,41) = REGR_4D(:,:,45,T)
746 0 : Lct%Dct%Dta%V3(T)%Val(:,:,42) = REGR_4D(:,:,49,T)
747 0 : Lct%Dct%Dta%V3(T)%Val(:,:,43) = REGR_4D(:,:,53,T)
748 0 : Lct%Dct%Dta%V3(T)%Val(:,:,44) = REGR_4D(:,:,57,T)
749 0 : Lct%Dct%Dta%V3(T)%Val(:,:,45) = REGR_4D(:,:,61,T)
750 0 : Lct%Dct%Dta%V3(T)%Val(:,:,46) = REGR_4D(:,:,65,T)
751 0 : Lct%Dct%Dta%V3(T)%Val(:,:,47) = REGR_4D(:,:,69,T)
752 0 : Lct%Dct%Dta%V3(T)%Val(:,:,48) = REGR_4D(:,:,73,T)
753 : ! If the input is N (N = 3/11/36) levels, levels N+1 to 47 are set to 0
754 0 : ELSEIF ( nlev == 3 ) THEN
755 0 : DO L = 4,47
756 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
757 : ENDDO !L
758 0 : ELSEIF ( nlev == 11 ) THEN
759 0 : DO L = 12,47
760 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
761 : ENDDO !L
762 0 : ELSEIF ( nlev == 36 ) THEN
763 0 : DO L = 37,47
764 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
765 : ENDDO !L
766 : ENDIF ! nlev == (72,73,36)
767 :
768 : ENDDO ! T
769 :
770 : ! Verbose
771 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
772 0 : WRITE(MSG,*) 'Mapped ', nlev, ' levels onto reduced GEOS-5 levels.'
773 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
774 0 : IF ( nlev == 3 .OR. nlev == 11 .OR. nlev == 36 ) THEN
775 0 : WRITE(MSG,*) 'The input variable has ', nlev, 'L, which were written to be L 1-', nlev, ' on the output 47 L grid (remaining values set to 0).'
776 : ELSE
777 0 : WRITE(MSG,*) 'Pressure-weighted vertical regridding was done - consider if this is appropriate for the variable units.'
778 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
779 : ENDIF
780 : ENDIF
781 :
782 : ! Done!
783 : DONE = .TRUE.
784 :
785 : !----------------------------------------------------------------
786 : ! Native to reduced GISS levels
787 : !----------------------------------------------------------------
788 0 : ELSEIF ( nz == 74 ) THEN
789 :
790 : ! Determine if the variable is on model levels or edges
791 0 : IF ( nlev == 102 ) THEN
792 0 : NL = 60
793 0 : nout = 74
794 0 : ELSEIF ( nlev == 103 ) THEN
795 0 : NL = 61
796 0 : nout = 75
797 : ELSE
798 : MSG = 'Can only remap from native onto reduced GISS if '// &
799 0 : 'input data has exactly 102 or 103 levels: '//TRIM(Lct%Dct%cName)
800 0 : CALL HCO_ERROR( MSG, RC )
801 0 : RETURN
802 : ENDIF ! nlev == (102,103,ELSE)
803 :
804 : ! Make sure output array is allocated
805 0 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )
806 :
807 : ! Do for every time slice
808 0 : DO T = 1, nt
809 :
810 : ! Levels that are passed level-by-level.
811 0 : DO L = 1, NL
812 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
813 : ENDDO !L
814 :
815 : ! If remapping model grid layers, collapse layers
816 0 : IF ( nlev == 102 ) THEN
817 : ! Collapse two levels (e.g. levels 61-62 into level 61):
818 0 : CALL COLLAPSE( Lct, REGR_4D, 61, 61, 2, T, 22, RC )
819 0 : CALL COLLAPSE( Lct, REGR_4D, 62, 63, 2, T, 22, RC )
820 0 : CALL COLLAPSE( Lct, REGR_4D, 63, 65, 2, T, 22, RC )
821 0 : CALL COLLAPSE( Lct, REGR_4D, 64, 67, 2, T, 22, RC )
822 0 : CALL COLLAPSE( Lct, REGR_4D, 65, 69, 2, T, 22, RC )
823 0 : CALL COLLAPSE( Lct, REGR_4D, 66, 71, 2, T, 22, RC )
824 0 : CALL COLLAPSE( Lct, REGR_4D, 67, 73, 2, T, 22, RC )
825 : ! Collapse four levels:
826 0 : CALL COLLAPSE( Lct, REGR_4D, 68, 75, 4, T, 22, RC )
827 0 : CALL COLLAPSE( Lct, REGR_4D, 69, 79, 4, T, 22, RC )
828 0 : CALL COLLAPSE( Lct, REGR_4D, 70, 83, 4, T, 22, RC )
829 0 : CALL COLLAPSE( Lct, REGR_4D, 71, 87, 4, T, 22, RC )
830 0 : CALL COLLAPSE( Lct, REGR_4D, 72, 91, 4, T, 22, RC )
831 0 : CALL COLLAPSE( Lct, REGR_4D, 73, 95, 4, T, 22, RC )
832 0 : CALL COLLAPSE( Lct, REGR_4D, 74, 99, 4, T, 22, RC )
833 : ! If remapping model grid edges, sample at the edges
834 : ELSE
835 0 : Lct%Dct%Dta%V3(T)%Val(:,:,62) = REGR_4D(:,:,63,T)
836 0 : Lct%Dct%Dta%V3(T)%Val(:,:,63) = REGR_4D(:,:,65,T)
837 0 : Lct%Dct%Dta%V3(T)%Val(:,:,64) = REGR_4D(:,:,67,T)
838 0 : Lct%Dct%Dta%V3(T)%Val(:,:,65) = REGR_4D(:,:,69,T)
839 0 : Lct%Dct%Dta%V3(T)%Val(:,:,66) = REGR_4D(:,:,71,T)
840 0 : Lct%Dct%Dta%V3(T)%Val(:,:,67) = REGR_4D(:,:,73,T)
841 0 : Lct%Dct%Dta%V3(T)%Val(:,:,68) = REGR_4D(:,:,75,T)
842 0 : Lct%Dct%Dta%V3(T)%Val(:,:,69) = REGR_4D(:,:,79,T)
843 0 : Lct%Dct%Dta%V3(T)%Val(:,:,70) = REGR_4D(:,:,83,T)
844 0 : Lct%Dct%Dta%V3(T)%Val(:,:,71) = REGR_4D(:,:,87,T)
845 0 : Lct%Dct%Dta%V3(T)%Val(:,:,72) = REGR_4D(:,:,91,T)
846 0 : Lct%Dct%Dta%V3(T)%Val(:,:,73) = REGR_4D(:,:,95,T)
847 0 : Lct%Dct%Dta%V3(T)%Val(:,:,74) = REGR_4D(:,:,99,T)
848 0 : Lct%Dct%Dta%V3(T)%Val(:,:,75) = REGR_4D(:,:,103,T)
849 : ENDIF ! nlev == (102,ELSE)
850 : ENDDO ! T
851 :
852 : ! Verbose
853 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
854 0 : WRITE(MSG,*) 'Mapped ', nlev, ' levels onto reduced GISS levels.'
855 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
856 0 : WRITE(MSG,*) 'Pressure-weighted vertical regridding was done - consider if this is appropriate for the variable units.'
857 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
858 : ENDIF
859 :
860 : ! Done!
861 : DONE = .TRUE.
862 :
863 : !----------------------------------------------------------------
864 : ! Reduced to native GEOS-5 levels
865 : !----------------------------------------------------------------
866 0 : ELSEIF ( nz == 72 ) THEN
867 :
868 : ! Determine if the variable is on model levels or edges
869 : IF ( nlev == 47 ) THEN
870 0 : NL = 36
871 0 : nout = 72
872 : ELSEIF ( nlev == 48 ) THEN
873 0 : NL = 37
874 0 : nout = 73
875 : ELSEIF ( nlev == 3 ) THEN
876 0 : NL = 3
877 0 : nout = 72
878 : ELSEIF ( nlev == 11 ) THEN
879 0 : NL = 11
880 0 : nout = 72
881 : ELSEIF ( nlev == 36 ) THEN
882 0 : NL = 36
883 0 : nout = 72
884 : ELSE
885 : MSG = 'Can only remap from reduced onto native GEOS-5 if '// &
886 0 : 'input data has exactly 47, 48, 3, 11, or 36 levels: '//TRIM(Lct%Dct%cName)
887 0 : CALL HCO_ERROR( MSG, RC )
888 0 : RETURN
889 : ENDIF ! nlev == (48,48,3,11,36,ELSE)
890 :
891 : ! Make sure output array is allocated
892 0 : CALL FileData_ArrCheck( HcoState%Config, Lct%Dct%Dta, nx, ny, nout, nt, RC )
893 :
894 : ! Do for every time slice
895 0 : DO T = 1, nt
896 :
897 : ! Levels that are passed level-by-level.
898 0 : DO L = 1, NL
899 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = REGR_4D(:,:,L,T)
900 : ENDDO !L
901 :
902 : ! If remapping model grid layers, inflate layers
903 0 : IF ( nlev == 47 ) THEN
904 :
905 : ! Inflate two levels (e.g. levels 37-38 on the fine grid are copies of level 37 on the coarse grid):
906 : coarseIDX = 36
907 0 : DO I = 1,8
908 0 : fineIDX = 36 + I
909 0 : IF ( MOD(I,2) /= 0 ) THEN
910 0 : coarseIDX = coarseIDX + 1
911 : ENDIF
912 0 : Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
913 : ENDDO ! I
914 :
915 : ! Inflate four levels (e.g. levels 45-48 on the fine grid are copies of level 41 on the coarse grid)
916 : coarseIDX = 40
917 0 : DO I = 1,28
918 0 : fineIDX = 44 + i
919 0 : IF ( MOD(I-1,4) == 0 ) THEN
920 0 : coarseIDX = coarseIDX + 1
921 : ENDIF
922 0 : Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
923 : ENDDO ! I
924 :
925 : ! If remapping model grid edges, inflate edges
926 0 : ELSEIF ( nlev == 48 ) THEN
927 :
928 : ! Sample every two edges (e.g. edges 38-39 on the fine grid are copies of edge 38 on the coarse grid)
929 : coarseIDX = 37
930 0 : DO I = 1,8
931 0 : fineIDX = 37 + I
932 0 : IF ( MOD(I,2) /= 0 ) THEN
933 0 : coarseIDX = coarseIDX + 1
934 : ENDIF
935 0 : Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
936 : ENDDO ! I
937 :
938 : ! Sample every four edges (e.g. edges 46-49 on the fine grid are copies of edge 42 on the coarse grid)
939 : coarseIDX = 40
940 0 : DO I = 1,28
941 0 : fineIDX = 44 + i
942 0 : IF ( MOD(I-1,4) == 0 ) THEN
943 0 : coarseIDX = coarseIDX + 1
944 : ENDIF
945 0 : Lct%Dct%Dta%V3(T)%Val(:,:,fineIDX) = REGR_4D(:,:,coarseIDX,T)
946 : ENDDO ! I
947 :
948 0 : ELSEIF ( nlev == 3 ) THEN
949 0 : DO L = 4,72
950 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
951 : ENDDO !L
952 :
953 0 : ELSEIF ( nlev == 11 ) THEN
954 0 : DO L = 12,72
955 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
956 : ENDDO !L
957 :
958 0 : ELSEIF ( nlev == 36 ) THEN
959 0 : DO L = 37,72
960 0 : Lct%Dct%Dta%V3(T)%Val(:,:,L) = 0.0_hp
961 : ENDDO !L
962 :
963 : ENDIF ! nlev == (47,48,3,11,36)
964 : ENDDO ! T
965 :
966 : ! Verbose
967 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
968 0 : WRITE(MSG,*) 'Mapped ', nlev, ' levels onto native GEOS-5 levels.'
969 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
970 0 : IF ( nlev == 3 .OR. nlev == 11 .OR. nlev == 36 ) THEN
971 0 : WRITE(MSG,*) 'The input variable has ', nlev, 'L, which were written to be L 1-', nlev, ' on the output 72 L grid (remaining values set to 0).'
972 : ELSE
973 0 : WRITE(MSG,*) 'Inflating from 47/48 to 72/73 levels is not recommended and is likely not mass-conserving.'
974 : ENDIF
975 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
976 : ENDIF
977 :
978 : ! Done!
979 : DONE = .TRUE.
980 :
981 : ENDIF ! nz == (47,74,72)
982 :
983 : ENDIF ! Vertical regridding required
984 :
985 : !===================================================================
986 : ! Error check / verbose mode
987 : !===================================================================
988 : IF ( DONE ) THEN
989 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) ) THEN
990 0 : WRITE(MSG,*) 'Did vertical regridding for ',TRIM(Lct%Dct%cName),':'
991 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
992 0 : WRITE(MSG,*) 'Number of original levels: ', nlev
993 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
994 0 : WRITE(MSG,*) 'Number of output levels: ', SIZE(Lct%Dct%Dta%V3(1)%Val,3)
995 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
996 : ENDIF
997 : ELSE
998 0 : WRITE(MSG,*) 'Vertical regridding failed: ',TRIM(Lct%Dct%cName)
999 0 : CALL HCO_ERROR( MSG, RC )
1000 0 : RETURN
1001 : ENDIF
1002 :
1003 : ! Return w/ success
1004 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
1005 :
1006 : END SUBROUTINE ModelLev_Interpolate
1007 : !EOC
1008 : !------------------------------------------------------------------------------
1009 : ! Harmonized Emissions Component (HEMCO) !
1010 : !------------------------------------------------------------------------------
1011 : !BOP
1012 : !
1013 : ! !IROUTINE: COLLAPSE
1014 : !
1015 : ! !DESCRIPTION: Helper routine to collapse input levels onto the output grid.
1016 : ! The input data is weighted by the grid box thicknesses defined on top of
1017 : ! this module. The input parameter T determines the time slice to be considered,
1018 : ! and MET denotes the met field type of the input data (22 = GISS, 5= GEOS-5).
1019 : !\\
1020 : !\\
1021 : ! !INTERFACE:
1022 : !
1023 0 : SUBROUTINE COLLAPSE ( Lct, REGR_4D, OutLev, InLev1, NLEV, T, MET, RC )
1024 : !
1025 : ! !INPUT PARAMETERS:
1026 : !
1027 : REAL(sp), POINTER :: REGR_4D(:,:,:,:) ! 4D input data
1028 : INTEGER, INTENT(IN) :: OutLev
1029 : INTEGER, INTENT(IN) :: InLev1
1030 : INTEGER, INTENT(IN) :: NLEV
1031 : INTEGER, INTENT(IN) :: T
1032 : INTEGER, INTENT(IN) :: MET ! 22=GISS E2.2, else GEOS-5
1033 : !
1034 : ! !INPUT/OUTPUT PARAMETERS:
1035 : !
1036 : TYPE(ListCont), POINTER :: Lct ! HEMCO list container
1037 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
1038 : !
1039 : ! !REVISION HISTORY:
1040 : ! 30 Dec 2014 - C. Keller - Initial version
1041 : ! See https://github.com/geoschem/hemco for complete history
1042 : !EOP
1043 : !------------------------------------------------------------------------------
1044 : !BOC
1045 : INTEGER :: I, NZ, ILEV, TOPLEV
1046 : REAL(hp) :: THICK
1047 0 : REAL(hp), POINTER :: EDG(:)
1048 0 : REAL(hp), ALLOCATABLE :: WGT(:)
1049 : CHARACTER(LEN=255) :: MSG
1050 : CHARACTER(LEN=255) :: LOC = 'ModelLev_Interpolate (hco_interp_mod.F90)'
1051 :
1052 : !=================================================================
1053 : ! COLLAPSE begins here
1054 : !=================================================================
1055 :
1056 : ! Init
1057 0 : EDG => NULL()
1058 :
1059 : ! Reset
1060 0 : Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = 0.0_hp
1061 :
1062 : ! Don't do anything if there are not enough levels in REGR_4D
1063 0 : NZ = SIZE(REGR_4D,3)
1064 0 : IF ( NZ < InLev1 ) RETURN
1065 :
1066 : ! Get maximum level to be used for pressure thickness calculations.
1067 0 : TOPLEV = InLev1 + NLEV
1068 :
1069 : ! Get pointer to grid edges on the native input grid
1070 0 : IF ( Met == 22 ) THEN
1071 0 : EDG => E102_EDGE_NATIVE(InLev1:TOPLEV)
1072 0 : ELSEIF ( Met == 5 ) THEN
1073 0 : EDG => G5_EDGE_NATIVE(InLev1:TOPLEV)
1074 : ELSE
1075 0 : WRITE(MSG,*) 'The Met value given was not valid: ', Met
1076 0 : CALL HCO_ERROR( MSG, RC, THISLOC=LOC )
1077 0 : RETURN
1078 : ENDIF
1079 :
1080 : ! Thickness of output level
1081 0 : THICK = EDG(1) - EDG(NLEV+1)
1082 :
1083 : ! Get level weights
1084 0 : ALLOCATE(WGT(NLEV))
1085 0 : WGT = 0.0
1086 0 : DO I = 1, NLEV
1087 0 : WGT(I) = ( EDG(I) - EDG(I+1) ) / THICK
1088 : ENDDO
1089 :
1090 : ! Pass levels to output data, one after each other
1091 0 : Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = REGR_4D(:,:,InLev1,T) * WGT(1)
1092 0 : DO I = 1, NLEV-1
1093 0 : ILEV = InLev1 + I
1094 0 : IF ( NZ < ILEV ) EXIT
1095 0 : Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) = Lct%Dct%Dta%V3(T)%Val(:,:,OutLev) &
1096 0 : + ( REGR_4D(:,:,ILEV,T) * WGT(I+1) )
1097 : ENDDO
1098 :
1099 : ! Cleanup
1100 0 : DEALLOCATE(WGT)
1101 0 : EDG => NULL()
1102 :
1103 0 : END SUBROUTINE COLLAPSE
1104 : !EOC
1105 : END MODULE HCO_Interp_Mod
|
