Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hcoio_messy_mod.F90
7 : !
8 : ! !DESCRIPTION: Module HCOIO\_MESSY\_MOD interfaces HEMCO with the regridding
9 : ! tool NCREGRID of the Modular Earth Submodel System (MESSy). This regridding
10 : ! scheme is used for vertical regridding and/or for index data, i.e. data with
11 : ! discrete values (e.g. land type integers). This code currently only works for
12 : ! rectilinear (regular lon-lat) grids but can be extended to support curvilinear
13 : ! grids.
14 : !\\
15 : ! REFERENCES:
16 : ! \begin{itemize}
17 : ! \item Joeckel, P. Technical note: Recursive rediscretisation of geo-
18 : ! scientific data in the Modular Earth Submodel System (MESSy), ACP, 6,
19 : ! 3557--3562, 2006.
20 : ! \end{itemize}
21 : ! !INTERFACE:
22 : !
23 : MODULE HCOIO_MESSY_MOD
24 : !
25 : ! !USES:
26 : !
27 : USE HCO_ERROR_MOD
28 : USE HCO_TYPES_MOD, ONLY : ListCont
29 : USE HCO_STATE_MOD, ONLY : Hco_State
30 : USE MESSY_NCREGRID_BASE, ONLY : NARRAY, AXIS
31 :
32 : IMPLICIT NONE
33 : PRIVATE
34 : !
35 : ! !PUBLIC MEMBER FUNCTIONS:
36 : !
37 : PUBLIC :: HCO_MESSY_REGRID
38 : !
39 : ! !PRIVATE MEMBER FUNCTIONS:
40 : !
41 : ! Set this value to TRUE if you want to reduce the output array
42 : ! to the minimum required number of vertical levels.
43 : LOGICAL, PARAMETER :: ReduceVert = .FALSE.
44 : !
45 : ! !MODULE INTERFACES:
46 : !
47 : ! !REVISION HISTORY:
48 : ! 24 Jun 2014 - C. Keller - Initial version
49 : ! See https://github.com/geoschem/hemco for complete history
50 : !EOP
51 : !------------------------------------------------------------------------------
52 : !BOC
53 : !
54 : ! !MODULE VARIABLES:
55 :
56 : CONTAINS
57 : !EOC
58 : !------------------------------------------------------------------------------
59 : ! Harmonized Emissions Component (HEMCO) !
60 : !------------------------------------------------------------------------------
61 : !BOP
62 : !
63 : ! !IROUTINE: Hco_Messy_Regrid
64 : !
65 : ! !DESCRIPTION: This is the wrapper routine to regrid a 4D input array
66 : ! NcArr (x,y,z,t) onto the HEMCO emissions grid (defined in HcoState)
67 : ! using the regridding tool NCREGRID. LonEdge, LatEdge and LevEdge are
68 : ! the grid point edges of the input grid. The data is written into list
69 : ! container Lct.
70 : !\\
71 : !\\
72 : ! If the input grid is 2D (horizontal only), LevEdge must not be specified
73 : ! (null pointer) and the data is regridded in the horizontal only. If the
74 : ! input grid has only one vertical level, it is assumed that this is the
75 : ! surface level and the output data is 3D but with only one vertical level.
76 : !\\
77 : !\\
78 : ! For input data with more than one vertical level, the data is mapped
79 : ! onto the entire 3D grid. The module parameter ReduceVert can be used to
80 : ! cap the output data at the lowest possible level. For example, if the
81 : ! input grid only covers three surface levels with a minimum sigma value
82 : ! of 0.75, vertical regridding is performed within this sigma range (1-0.75)
83 : ! and the output grid is reduced accordingly. This option is not used in the
84 : ! standard HEMCO setup because problems can arise if the data array of a
85 : ! given container suddently changes its size (i.e. when updated data covers
86 : ! more/less vertical levels than the data beforehand).
87 : !\\
88 : !\\
89 : ! The input argument IsModelLev denotes whether or not the vertical
90 : ! coordinates of the input data are on model levels. If set to yes and LevEdge
91 : ! is not provided (i.e. a nullified pointer), the MESSy regridding routines are
92 : ! only used for the horizontal remapping and subroutine ModelLev\_Interpolate
93 : ! (module hco\_interp\_mod.F90) is used for the vertical remapping.
94 : !\\
95 : !\\
96 : ! !INTERFACE:
97 : !
98 0 : SUBROUTINE HCO_MESSY_REGRID ( HcoState, NcArr, &
99 : LonEdge, LatEdge, LevEdge, &
100 : Lct, IsModelLev, RC )
101 : !
102 : ! !USES:
103 : !
104 : USE HCO_FILEDATA_MOD, ONLY : FileData_ArrCheck
105 : USE HCO_UNIT_MOD, ONLY : HCO_IsIndexData
106 : USE HCO_INTERP_MOD, ONLY : ModelLev_Interpolate
107 : USE MESSY_NCREGRID_BASE, ONLY : RG_INT, RG_IDX
108 : USE MESSY_NCREGRID_BASE, ONLY : NREGRID
109 : USE MESSY_NCREGRID_BASE, ONLY : INIT_NARRAY
110 : !
111 : ! !INPUT/OUTPUT PARAMETERS:
112 : !
113 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO obj.
114 : REAL(sp), POINTER :: ncArr(:,:,:,:) ! Input array(x,y,z,t)
115 : REAL(hp), POINTER :: LonEdge(:) ! lon edges
116 : REAL(hp), POINTER :: LatEdge(:) ! lat edges
117 : REAL(hp), POINTER :: LevEdge(:,:,:) ! sigma level edges
118 : TYPE(ListCont), POINTER :: Lct ! Target list container
119 : LOGICAL, INTENT(IN ) :: IsModelLev ! Are these model levels?
120 : INTEGER, INTENT(INOUT) :: RC ! Return code
121 : !
122 : ! !REVISION HISTORY:
123 : ! 27 Jun 2014 - C. Keller - Initial version
124 : ! See https://github.com/geoschem/hemco for complete history
125 : !EOP
126 : !------------------------------------------------------------------------------
127 : !BOC
128 : !
129 : ! !ROUTINE ARGUMENTS:
130 : !
131 0 : TYPE(narray), POINTER :: narr_src(:)
132 0 : TYPE(narray), POINTER :: narr_dst(:)
133 0 : TYPE(axis), POINTER :: axis_src(:)
134 0 : TYPE(axis), POINTER :: axis_dst(:)
135 0 : INTEGER, POINTER :: rg_type(:)
136 0 : INTEGER, POINTER :: rcnt (:)
137 0 : REAL(dp), POINTER :: sovl (:)
138 0 : REAL(dp), POINTER :: dovl (:)
139 0 : REAL(hp), POINTER :: lon (:)
140 0 : REAL(hp), POINTER :: lat (:)
141 0 : REAL(hp), POINTER :: sigma (:,:,:)
142 0 : REAL(sp), POINTER :: ArrIn (:,:,:,:)
143 0 : REAL(sp), POINTER :: ArrOut (:,:,:,:)
144 0 : REAL(hp), ALLOCATABLE, TARGET :: sigout (:,:,:)
145 : REAL(hp) :: sigMin
146 : INTEGER :: NZIN, NZOUT, NTIME
147 : INTEGER :: NXIN, NYIN
148 : INTEGER :: I, L, AS
149 : INTEGER :: NCALLS
150 : CHARACTER(LEN=255) :: MSG, LOC
151 : LOGICAL :: SameGrid, verb
152 :
153 : !=================================================================
154 : ! HCO_MESSY_REGRID begins here
155 : !=================================================================
156 :
157 : ! For error handling
158 0 : LOC = 'HCO_MESSY_REGRID (HCOIO_MESSY_MOD.F90)'
159 0 : CALL HCO_ENTER ( HcoState%Config%Err, LOC, RC )
160 0 : IF ( RC /= HCO_SUCCESS ) THEN
161 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
162 0 : RETURN
163 : ENDIF
164 :
165 : ! Init
166 0 : narr_src => NULL()
167 0 : narr_dst => NULL()
168 0 : axis_src => NULL()
169 0 : axis_dst => NULL()
170 0 : rg_type => NULL()
171 0 : rcnt => NULL()
172 0 : sovl => NULL()
173 0 : dovl => NULL()
174 0 : lon => NULL()
175 0 : lat => NULL()
176 0 : sigma => NULL()
177 0 : ArrIn => NULL()
178 0 : ArrOut => NULL()
179 :
180 : ! verbose?
181 0 : verb = HCO_IsVerb( HcoState%Config%Err )
182 :
183 : ! Horizontal dimension of input data
184 0 : NXIN = SIZE(NcArr,1)
185 0 : NYIN = SIZE(NcArr,2)
186 :
187 : ! Number of vertical levels of input data
188 0 : NZIN = SIZE(NcArr,3)
189 :
190 : ! Number of time slices. All time slices will be regridded
191 : ! simultaneously.
192 0 : NTIME = SIZE(NcArr,4)
193 :
194 : ! Error check: data must be 2D or 3D.
195 0 : IF ( Lct%Dct%Dta%SpaceDim /= 2 .AND. &
196 : Lct%Dct%Dta%SpaceDim /= 3 ) THEN
197 0 : MSG = 'Can only regrid 2D or 3D data: ' // TRIM(Lct%Dct%cName)
198 0 : CALL HCO_ERROR ( MSG, RC )
199 0 : RETURN
200 : ENDIF
201 :
202 : !-----------------------------------------------------------------
203 : ! Shortcut if input field is already on output grid: directly
204 : ! pass data to Lct.
205 : ! NOTE: if the number of input levels matches the number of output
206 : ! levels (and the horizontal dimensions agree as well), it is
207 : ! assumed that they are on the same grid!
208 : !-----------------------------------------------------------------
209 :
210 : ! Input grid = output grid?
211 0 : SameGrid = .FALSE.
212 :
213 : ! Horizontal dimensions have to match
214 0 : IF ( (NXIN == HcoState%NX) .AND. (NYIN == HcoState%NY) ) THEN
215 :
216 : ! Vertical dimensions have to match or be 1
217 0 : IF ( NZIN == 1 .OR. NZIN == HcoState%NZ ) THEN
218 :
219 : ! Assume same grid
220 0 : SameGrid = .TRUE.
221 :
222 : ! Check for same boundaries. Otherwise falsify SameGrid
223 : ! Lon ...
224 0 : IF ( MINVAL(LonEdge) /= MINVAL(HcoState%Grid%XEDGE%Val) .OR. &
225 : MAXVAL(LonEdge) /= MAXVAL(HcoState%Grid%XEDGE%Val) ) THEN
226 0 : SameGrid = .FALSE.
227 : ENDIF
228 : ! ... Lat ...
229 0 : IF ( MINVAL(LatEdge) /= MINVAL(HcoState%Grid%YEDGE%Val) .OR. &
230 : MAXVAL(LatEdge) /= MAXVAL(HcoState%Grid%YEDGE%Val) ) THEN
231 0 : SameGrid = .FALSE.
232 : ENDIF
233 : ! ... Lev
234 : IF ( NZIN > 1 ) THEN
235 : ! TODO: Eventually need to add level check here.
236 : ! For now, assume that input levels are equal to
237 : ! output level if they have the same dimension!
238 : ENDIF
239 : ENDIF
240 : ENDIF
241 :
242 : !-----------------------------------------------------------------
243 : ! Define number of vertical levels on output grid
244 : !-----------------------------------------------------------------
245 :
246 : ! Error check. If we are not on the same grid, LevEdge must be
247 : ! provided or IsModelLev must be set to true for 3D data.
248 0 : IF ( NZIN > 1 .AND. .NOT. SameGrid ) THEN
249 0 : IF ( .NOT. ASSOCIATED(LevEdge) .AND. .NOT. IsModelLev ) THEN
250 : MSG = 'Cannot regrid '//TRIM(Lct%Dct%cName)//'. Either level '//&
251 0 : 'edges must be provided or data must be on model levels.'
252 0 : CALL HCO_ERROR ( MSG, RC )
253 0 : RETURN
254 : ENDIF
255 : ENDIF
256 :
257 0 : NZOUT = NZIN
258 0 : IF ( .NOT. SameGrid .AND. ASSOCIATED(LevEdge) ) THEN
259 :
260 : ! Calculate sigma level for each grid point on the output grid.
261 : ! This is the pressure at location i,j,l normalized by surface
262 : ! pressure @ i,j: sigma(i,j,l) = p(i,j,l) / ps(i,j)
263 0 : ALLOCATE(sigout(HcoState%NX,HcoState%NY,HcoState%NZ+1),STAT=AS)
264 0 : IF ( AS/= 0 ) THEN
265 0 : CALL HCO_ERROR( 'Cannot allocate sigout', RC )
266 0 : RETURN
267 : ENDIF
268 0 : DO l = 1, HcoState%NZ+1
269 0 : sigout(:,:,l) = HcoState%Grid%PEDGE%Val(:,:,l) &
270 0 : / HcoState%Grid%PEDGE%Val(:,:,1)
271 : ENDDO
272 :
273 : ! Now find first level on output grid where all sigma levels are
274 : ! lower than the lowest sigma value on the input grid. This is
275 : ! the highest level that needs to be considered for regridding.
276 : IF ( ReduceVert ) THEN
277 : sigMin = MINVAL(LevEdge)
278 : DO l = 1, HcoState%NZ+1
279 : IF ( MINVAL(sigout(:,:,l)) < sigMin ) EXIT
280 : ENDDO
281 :
282 : ! The output grid is at grid center, so use l-1. Must be at least one.
283 : NZOUT = max(1,l-1)
284 :
285 : ! Use full vertical grid if vertical levels shall not be restricted
286 : ! to range of input data (default).
287 : ELSE
288 0 : NZOUT = HcoState%NZ
289 : ENDIF
290 : ENDIF
291 :
292 : ! verbose mode
293 0 : IF ( verb ) THEN
294 0 : MSG = 'Do MESSy regridding: ' // TRIM(Lct%Dct%cName)
295 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
296 0 : WRITE(MSG,*) ' - SameGrid ? ', SameGrid
297 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
298 0 : WRITE(MSG,*) ' - Model levels ? ', IsModelLev
299 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
300 : ENDIF
301 :
302 : !-----------------------------------------------------------------
303 : ! Make sure output array is defined & allocated
304 : !-----------------------------------------------------------------
305 0 : IF ( Lct%Dct%Dta%SpaceDim == 2 ) THEN
306 : CALL FileData_ArrCheck( HcoState%Config, &
307 : Lct%Dct%Dta, HcoState%NX, HcoState%NY, &
308 0 : NTIME, RC )
309 0 : IF ( RC /= HCO_SUCCESS ) THEN
310 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
311 0 : RETURN
312 : ENDIF
313 0 : ELSEIF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
314 : CALL FileData_ArrCheck( HcoState%Config, &
315 : Lct%Dct%Dta, HcoState%NX, HcoState%NY, &
316 0 : NZOUT, NTIME, RC )
317 0 : IF ( RC /= HCO_SUCCESS ) THEN
318 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
319 0 : RETURN
320 : ENDIF
321 : ENDIF
322 :
323 : !-----------------------------------------------------------------
324 : ! Do straight-forward mapping if input grid = output grid
325 : !-----------------------------------------------------------------
326 0 : IF ( SameGrid ) THEN
327 : MSG = 'Input grid seems to match output grid. ' // &
328 0 : 'No regridding is performed: ' // TRIM(Lct%Dct%cName)
329 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC )
330 :
331 : ! For every time slice...
332 0 : DO I = 1, NTIME
333 0 : DO L = 1, NZOUT
334 :
335 : ! 3D data
336 0 : IF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
337 0 : Lct%Dct%Dta%V3(I)%Val(:,:,L) = NcArr(:,:,L,I)
338 :
339 : ! 2D data
340 : ELSE
341 0 : Lct%Dct%Dta%V2(I)%Val(:,:) = NcArr(:,:,L,I)
342 : ENDIF
343 : ENDDO
344 : ENDDO
345 :
346 : ! All done!
347 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
348 0 : RETURN
349 : ENDIF
350 :
351 : !=================================================================
352 : ! MESSy regridding follows below
353 : !=================================================================
354 :
355 : !-----------------------------------------------------------------
356 : ! Source grid description.
357 : ! This creates a MESSy axis object for the source grid.
358 : !-----------------------------------------------------------------
359 0 : lon => LonEdge
360 0 : lat => LatEdge
361 0 : sigma => LevEdge
362 :
363 0 : CALL AXIS_CREATE( HcoState, lon, lat, sigma, axis_src, RC )
364 0 : IF ( RC /= HCO_SUCCESS ) THEN
365 0 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
366 0 : RETURN
367 : ENDIF
368 :
369 : ! Free pointer
370 0 : lon => NULL()
371 0 : lat => NULL()
372 0 : sigma => NULL()
373 :
374 : !-----------------------------------------------------------------
375 : ! Destination grid description.
376 : ! This creates a MESSy axis object for the target (=HEMCO) grid.
377 : !-----------------------------------------------------------------
378 :
379 : ! Get horizontal grid directly from HEMCO state
380 0 : lon => HcoState%Grid%XEDGE%Val(:,1)
381 0 : lat => HcoState%Grid%YEDGE%Val(1,:)
382 0 : IF( ASSOCIATED(LevEdge) ) sigma => sigout(:,:,1:NZOUT+1)
383 :
384 0 : CALL AXIS_CREATE( HcoState, lon, lat, sigma, axis_dst, RC )
385 0 : IF ( RC /= HCO_SUCCESS ) THEN
386 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
387 0 : RETURN
388 : ENDIF
389 :
390 : ! Free pointer
391 0 : lon => NULL()
392 0 : lat => NULL()
393 0 : sigma => NULL()
394 0 : IF ( ALLOCATED(sigout) ) DEALLOCATE(sigout)
395 :
396 : !-----------------------------------------------------------------
397 : ! Set all other regridding parameter
398 : !-----------------------------------------------------------------
399 :
400 : ! rg_type denotes the regridding type for each array (i.e. time
401 : ! slice). Set to 'intensive quantity' for all concentrations (incl.
402 : ! unitless) data. Set to 'index distribution' for data marked as
403 : ! index data in the configuration file. This will remap discrete
404 : ! values without interpolation, i.e. each grid box on the new
405 : ! grid holds the value with most overlap in the original grid.
406 0 : ALLOCATE(rg_type(NTIME))
407 0 : IF ( HCO_IsIndexData(Lct%Dct%Dta%OrigUnit) ) THEN
408 0 : rg_type(:) = RG_IDX
409 0 : IF ( verb ) THEN
410 0 : MSG = ' - Remap as index data.'
411 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
412 : ENDIF
413 : ELSE
414 0 : rg_type(:) = RG_INT
415 0 : IF ( verb ) THEN
416 0 : MSG = ' - Remap as concentration data.'
417 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
418 : ENDIF
419 : ENDIF
420 :
421 : !-----------------------------------------------------------------
422 : ! Number of times the regridding need to be performed. If input
423 : ! data is on model levels, the data is only horizontally regridded,
424 : ! e.g. the regridding routine is called for every horizontal level
425 : ! separately. The vertical interpolation is done afterwards using
426 : ! routine ModelLev_Interpolate.
427 : !-----------------------------------------------------------------
428 0 : NCALLS = 1
429 0 : IF ( IsModelLev .AND. .NOT. ASSOCIATED(LevEdge) .AND. NZIN > 1 ) THEN
430 0 : NCALLS = NZIN
431 0 : ALLOCATE(ArrOut(HcoState%NX,HcoState%NY,NZIN,NTIME),STAT=AS)
432 : IF ( AS /= 0 ) THEN
433 0 : CALL HCO_ERROR( 'Cannot allocate ArrOut', RC )
434 0 : RETURN
435 : ENDIF
436 0 : ArrOut = 0.0_sp
437 : ENDIF
438 :
439 : ! Do for all level batches ...
440 0 : DO I = 1, NCALLS
441 :
442 : ! ArrIn is the input array to be used
443 0 : IF ( NCALLS /= 1 ) THEN
444 0 : ArrIn => NcArr(:,:,I:I,:)
445 : ELSE
446 0 : ArrIn => NcArr
447 : ENDIF
448 :
449 : !-----------------------------------------------------------------
450 : ! Map input array onto MESSy array. Different time slices are
451 : ! stored as individual vector elements of narr_src.
452 : !-----------------------------------------------------------------
453 0 : CALL HCO2MESSY( HcoState, ArrIn, narr_src, axis_src, RC )
454 0 : IF ( RC /= HCO_SUCCESS ) THEN
455 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
456 0 : RETURN
457 : ENDIF
458 :
459 : !-----------------------------------------------------------------
460 : ! Do the regridding
461 : !-----------------------------------------------------------------
462 : CALL NREGRID(s=narr_src, sax=axis_src, dax=axis_dst, d=narr_dst, &
463 0 : rg_type=rg_type, sovl=sovl, dovl=dovl, rcnt=rcnt )
464 :
465 : !-----------------------------------------------------------------
466 : ! Map the destination array narr_dst onto the data vector in the
467 : ! HEMCO list container or onto the temporary array ArrOut.
468 : !-----------------------------------------------------------------
469 0 : CALL MESSY2HCO( HcoState, narr_dst, Lct, I, ArrOut, RC )
470 0 : IF ( RC /= HCO_SUCCESS ) THEN
471 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
472 0 : RETURN
473 : ENDIF
474 :
475 : ! Cleanup
476 0 : ArrIn => NULL()
477 :
478 : ENDDO !NCALLS
479 :
480 : !-----------------------------------------------------------------
481 : ! If these are model levels, do vertical interpolation now
482 : !-----------------------------------------------------------------
483 0 : IF ( ASSOCIATED(ArrOut) ) THEN
484 0 : CALL ModelLev_Interpolate( HcoState, ArrOut, Lct, RC )
485 0 : IF ( RC /= HCO_SUCCESS ) THEN
486 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
487 0 : RETURN
488 : ENDIF
489 0 : DEALLOCATE(ArrOut)
490 : ENDIF
491 :
492 : !-----------------------------------------------------------------
493 : ! Cleanup
494 : !-----------------------------------------------------------------
495 0 : DEALLOCATE( sovl, dovl, rcnt, STAT=AS)
496 0 : IF(AS/=0) THEN
497 0 : CALL HCO_ERROR('DEALLOCATION ERROR 1', RC )
498 0 : RETURN
499 : ENDIF
500 0 : NULLIFY(sovl, dovl, rcnt)
501 :
502 0 : DO I=1, SIZE(narr_dst)
503 0 : CALL INIT_NARRAY(narr_dst(I))
504 : ENDDO
505 0 : DEALLOCATE(narr_dst, STAT=AS)
506 : IF(AS/=0) THEN
507 0 : CALL HCO_ERROR('DEALLOCATION ERROR 3', RC )
508 0 : RETURN
509 : ENDIF
510 : NULLIFY(narr_dst)
511 :
512 0 : DO I=1, SIZE(narr_src)
513 0 : CALL INIT_NARRAY(narr_src(I))
514 : ENDDO
515 0 : DEALLOCATE(narr_src, STAT=AS)
516 : IF(AS/=0) THEN
517 0 : CALL HCO_ERROR('DEALLOCATION ERROR 2', RC )
518 0 : RETURN
519 : ENDIF
520 : NULLIFY(narr_src)
521 :
522 0 : DEALLOCATE( rg_type )
523 : rg_type => NULL()
524 :
525 0 : CALL AXIS_DELETE( axis_src, axis_dst, RC )
526 0 : IF ( RC /= HCO_SUCCESS ) THEN
527 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
528 0 : RETURN
529 : ENDIF
530 :
531 : ! Return w/ success
532 0 : CALL HCO_LEAVE ( HcoState%Config%Err, RC )
533 :
534 0 : END SUBROUTINE HCO_MESSY_REGRID
535 : !EOC
536 : !------------------------------------------------------------------------------
537 : ! Harmonized Emissions Component (HEMCO) !
538 : !------------------------------------------------------------------------------
539 : !BOP
540 : !
541 : ! !IROUTINE: Axis_Create
542 : !
543 : ! !DESCRIPTION: Subroutine AXIS\_CREATE creates a MESSy axis type
544 : ! from the grid defined by mid points Lon, Lat, Lev. Lev must be in
545 : ! (unitless) sigma coordinates: sigma(i,j,l) = p(i,j,l) / p\_surface(i,j)
546 : !\\
547 : !\\
548 : ! !INTERFACE:
549 : !
550 0 : SUBROUTINE AXIS_CREATE( HcoState, lon, lat, lev, ax, RC )
551 : !
552 : ! !USES:
553 : !
554 : USE MESSY_NCREGRID_BASE, ONLY : INIT_AXIS
555 : !
556 : ! !INPUT PARAMETERS:
557 : !
558 : TYPE(HCO_State), POINTER :: HcoState
559 : TYPE(ListCont), POINTER :: Lct
560 : REAL(hp), POINTER :: Lon(:)
561 : REAL(hp), POINTER :: Lat(:)
562 : REAL(hp), POINTER :: Lev(:,:,:)
563 : !
564 : ! !INPUT/OUTPUT PARAMETERS:
565 : !
566 : TYPE(axis), POINTER :: ax(:)
567 : INTEGER, INTENT(INOUT) :: RC
568 : !
569 : ! !REVISION HISTORY:
570 : ! 22 Jun 2014 - C. Keller - Initial version (from messy_ncregrid_geohyb.f90)
571 : ! See https://github.com/geoschem/hemco for complete history
572 : !EOP
573 : !------------------------------------------------------------------------------
574 : !BOC
575 : !
576 : ! !ROUTINE ARGUMENTS:
577 : !
578 : INTEGER :: fID, I, J, N, status
579 : INTEGER :: LOW, UPP
580 : INTEGER :: XLON, YLAT, XLEV, YLEV, ZLEV
581 : INTEGER :: ndp, nlev, cnt
582 : INTEGER :: vtype
583 : INTEGER :: ndep_lon, ndep_lat
584 : CHARACTER(LEN=255) :: MSG, LOC
585 :
586 : !=================================================================
587 : ! AXIS_CREATE begins here
588 : !=================================================================
589 :
590 : ! For error handling
591 0 : LOC = 'AXIS_CREATE (HCOIO_MESSY_MOD.F90)'
592 0 : CALL HCO_ENTER ( HcoState%Config%Err, LOC, RC )
593 0 : IF ( RC /= HCO_SUCCESS ) THEN
594 0 : CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
595 0 : RETURN
596 : ENDIF
597 :
598 : ! ----------------------------------------------------------------
599 : ! Pass horizontal grid dimensions to local variables. For now,
600 : ! assume all grids to be regular, i.e. no curvilinear grids.
601 : ! In this case, lon and lat dimensions are 1D-vectors.
602 : ! ----------------------------------------------------------------
603 :
604 : ! ndep_lon and ndep_lat are the axis number of lon and lat. Needed
605 : ! if we have to set dependencies for the vertical axis.
606 0 : ndep_lon = 0
607 0 : ndep_lat = 0
608 :
609 : ! Count axes and get grid dimensions
610 : ! Last axis is always 'free' dimension. This is required for the
611 : ! regridding to work properly.
612 0 : N = 1
613 0 : IF ( ASSOCIATED(lon) ) N = N + 1
614 0 : IF ( ASSOCIATED(lat) ) N = N + 1
615 0 : IF ( ASSOCIATED(lev) ) N = N + 1
616 :
617 : ! ALLOCATE AXIS
618 0 : ALLOCATE(ax(N), STAT=status)
619 0 : IF ( status /= 0 ) THEN
620 0 : CALL HCO_ERROR ( 'Cannot allocate axis', RC )
621 0 : RETURN
622 : ENDIF
623 0 : DO I=1, N
624 0 : CALL INIT_AXIS(ax(I))
625 : ENDDO
626 :
627 : ! N is the current axis
628 0 : N = 0
629 :
630 : ! ----------------------------------------------------------------
631 : ! Assign longitude: this is always the first dimension
632 : ! ----------------------------------------------------------------
633 0 : IF ( ASSOCIATED(lon) ) THEN
634 0 : N = N + 1
635 0 : ax(N)%lm = .true. ! LONGITUDE IS MODULO AXIS
636 :
637 : ! Axis dimension
638 0 : XLON = SIZE(lon,1)
639 :
640 : ! FOR NOW, ASSUME NO DEPENDENCIES. NEED TO EDIT HERE
641 : ! IF WE WANT TO USE CURVILINEAR GRIDS
642 0 : ax(N)%ndp = 1 ! LONGITUDE IS ...
643 0 : ALLOCATE(ax(N)%dep(1), STAT=status)
644 : IF ( status/= 0 ) THEN
645 0 : CALL HCO_ERROR ( 'Cannot allocate lon dependencies', RC )
646 0 : RETURN
647 : ENDIF
648 0 : ax(N)%dep(1) = N ! ... INDEPENDENT
649 0 : ndep_lon = N
650 :
651 0 : ax(N)%dat%n = 1 ! 1 dimension
652 0 : ALLOCATE(ax(N)%dat%dim(ax(N)%dat%n), STAT=status)
653 : IF ( status/= 0 ) THEN
654 0 : CALL HCO_ERROR( 'Cannot allocate lon dimensions', RC )
655 0 : RETURN
656 : ENDIF
657 0 : ax(N)%dat%dim(:) = XLON
658 :
659 0 : ALLOCATE(ax(N)%dat%vd(XLON),STAT=status)
660 : IF ( status/= 0 ) THEN
661 0 : CALL HCO_ERROR( 'Cannot allocate lon axis', RC )
662 0 : RETURN
663 : ENDIF
664 0 : ax(N)%dat%vd(:) = lon
665 :
666 : ENDIF !lon
667 :
668 : ! ----------------------------------------------------------------
669 : ! Assign latitude: this is always the second dimension
670 : ! ----------------------------------------------------------------
671 0 : IF ( ASSOCIATED(lat) ) THEN
672 0 : N = N + 1
673 0 : ax(N)%lm = .false. ! LATITUDE IS NON-MODULO AXIS
674 :
675 : ! Axis dimension
676 0 : YLAT = SIZE(lat,1)
677 :
678 : ! FOR NOW, ASSUME NO DEPENDENCIES. NEED TO EDIT HERE
679 : ! IF WE WANT TO USE CURVILINEAR GRIDS
680 0 : ax(N)%ndp = 1 ! LATITUDE IS ...
681 0 : ALLOCATE(ax(N)%dep(1), STAT=status)
682 : IF ( status/= 0 ) THEN
683 0 : CALL HCO_ERROR( 'Cannot allocate lat dependencies', RC )
684 0 : RETURN
685 : ENDIF
686 0 : ax(N)%dep(1) = N ! ... INDEPENDENT
687 0 : ndep_lat = N
688 :
689 0 : ax(N)%dat%n = 1 ! 1 dimension
690 0 : ALLOCATE(ax(N)%dat%dim(ax(N)%dat%n), STAT=status)
691 : IF ( status/= 0 ) THEN
692 0 : CALL HCO_ERROR( 'Cannot allocate lat dimensions', RC )
693 0 : RETURN
694 : ENDIF
695 0 : ax(N)%dat%dim(:) = YLAT
696 :
697 0 : ALLOCATE(ax(N)%dat%vd(YLAT),STAT=status)
698 : IF ( status/= 0 ) THEN
699 0 : CALL HCO_ERROR( 'Cannot allocate lat axis', RC )
700 0 : RETURN
701 : ENDIF
702 0 : ax(N)%dat%vd(:) = lat
703 :
704 : ! TAKE INTO ACCOUNT SPHERICAL GEOMETRY ...
705 0 : ax(N)%dat%vd = COS( ( (ax(N)%dat%vd - 90.0_dp) / 180.0_dp ) * &
706 0 : HcoState%Phys%PI )
707 :
708 : ENDIF !lat
709 :
710 : ! ----------------------------------------------------------------
711 : ! Assign vertical levels (if defined): this is the 3rd dimension.
712 : ! The vertical axis is assumed to be in sigma coordinates.
713 : ! ----------------------------------------------------------------
714 0 : IF ( ASSOCIATED(lev) ) THEN
715 :
716 : ! -------------------------------------------------------------
717 : ! Initialize vertical axis. Set dependencies on other axis.
718 : ! Define dimensions in the following order: lev, lat, lon.
719 : ! (first dimension has to be the 'own' dimension!).
720 : ! -------------------------------------------------------------
721 0 : N = N + 1
722 0 : ax(N)%lm = .false. ! VERTICAL AXIS IS NON-MODULO AXIS
723 :
724 : ! Axis dimension
725 0 : XLEV = SIZE(lev,1)
726 0 : YLEV = SIZE(lev,2)
727 0 : ZLEV = SIZE(lev,3)
728 :
729 : ! Sanity check: if XLEV and YLEV are > 1, they must correspond
730 : ! to the lon/lat axis defined above
731 0 : IF ( XLEV > 1 .AND. XLEV /= (XLON-1) ) THEN
732 0 : CALL HCO_ERROR ( 'level lon has wrong dimension', RC )
733 0 : RETURN
734 : ENDIF
735 0 : IF ( YLEV > 1 .AND. YLEV /= (YLAT-1) ) THEN
736 0 : CALL HCO_ERROR ( 'level lat has wrong dimension', RC )
737 0 : RETURN
738 : ENDIF
739 :
740 : ! Set dependencies. First dimension must be vertical axis!
741 0 : ndp = 1
742 0 : IF ( YLEV > 1 ) ndp = ndp + 1
743 0 : IF ( XLEV > 1 ) ndp = ndp + 1
744 :
745 0 : ax(N)%ndp = ndp
746 0 : ALLOCATE(ax(N)%dep(ax(N)%ndp), STAT=status)
747 : IF ( status /= 0 ) THEN
748 0 : CALL HCO_ERROR( 'Cannot allocate lev dependencies', RC )
749 0 : RETURN
750 : ENDIF
751 :
752 : ! The variable dat%n holds the number of axis that level depends
753 : ! upon, and dat%dim are the corresponding axis dimensions (lengths).
754 0 : ax(N)%dat%n = ndp
755 0 : ALLOCATE(ax(N)%dat%dim(ax(N)%dat%n), STAT=status)
756 : IF ( status/= 0 ) THEN
757 0 : CALL HCO_ERROR( 'Cannot allocate lat dimensions', RC )
758 0 : RETURN
759 : ENDIF
760 :
761 : ! Set axis indices and lengths:
762 : ! - First dimension is vertical axis
763 0 : cnt = 1
764 0 : ax(N)%dep(cnt) = N
765 0 : ax(N)%dat%dim(cnt) = ZLEV
766 0 : nlev = ZLEV ! number of grid cells
767 : ! - Next dimension is latitude
768 0 : IF ( YLEV > 1 ) THEN
769 0 : cnt = cnt + 1
770 0 : ax(N)%dep(cnt) = ndep_lat
771 0 : ax(N)%dat%dim(cnt) = YLEV
772 0 : nlev = nlev * YLEV
773 : ENDIF
774 : ! - Next dimension is longitude
775 0 : IF ( XLEV > 1 ) THEN
776 0 : cnt = cnt + 1
777 0 : ax(N)%dep(cnt) = ndep_lon
778 0 : ax(N)%dat%dim(cnt) = XLEV
779 0 : nlev = nlev * XLEV
780 : ENDIF
781 :
782 0 : ALLOCATE(ax(N)%dat%vd(nlev),STAT=status)
783 : IF ( status/= 0 ) THEN
784 0 : CALL HCO_ERROR( 'Cannot allocate lat axis', RC )
785 0 : RETURN
786 : ENDIF
787 :
788 : ! -------------------------------------------------------------
789 : ! Pass vertical sigma coordinates to vector
790 : ! -------------------------------------------------------------
791 : LOW = 1
792 0 : UPP = 1
793 0 : DO J = 1, YLEV !lat
794 0 : DO I = 1, XLEV !lon
795 0 : UPP = LOW + ZLEV - 1 ! Upper index
796 0 : ax(N)%dat%vd(LOW:UPP) = lev(I,J,:) ! Pass to vector
797 0 : LOW = UPP + 1 ! Next lower index
798 : ENDDO
799 : ENDDO
800 :
801 : END IF !lev
802 :
803 0 : CALL HCO_LEAVE ( HcoState%config%Err, RC )
804 :
805 : END SUBROUTINE AXIS_CREATE
806 : !EOC
807 : !------------------------------------------------------------------------------
808 : ! Harmonized Emissions Component (HEMCO) !
809 : !------------------------------------------------------------------------------
810 : !BOP
811 : !
812 : ! !IROUTINE: Axis_Delete
813 : !
814 : ! !DESCRIPTION: Subroutine AXIS\_DELETE deletes the specified MESSy
815 : ! axis.
816 : !\\
817 : !\\
818 : ! !INTERFACE:
819 : !
820 0 : SUBROUTINE AXIS_DELETE ( ax1, ax2, RC )
821 : !
822 : ! !USES:
823 : !
824 : USE MESSY_NCREGRID_BASE, ONLY : INIT_AXIS
825 : !
826 : ! !INPUT/OUTPUT PARAMETERS:
827 : !
828 : TYPE(axis), POINTER :: ax1(:)
829 : TYPE(axis), POINTER :: ax2(:)
830 : INTEGER, INTENT(INOUT) :: RC
831 : !
832 : ! !REVISION HISTORY:
833 : ! 28 Aug 2013 - C. Keller - Initial version
834 : ! See https://github.com/geoschem/hemco for complete history
835 : !EOP
836 : !------------------------------------------------------------------------------
837 : !BOC
838 : !
839 : ! !ROUTINE ARGUMENTS:
840 : !
841 : INTEGER :: I, status
842 :
843 : !=================================================================
844 : ! AXIS_DELETE begins here
845 : !=================================================================
846 :
847 0 : DO I=1, SIZE(ax1)
848 0 : CALL INIT_AXIS(ax1(I))
849 0 : CALL INIT_AXIS(ax2(I))
850 : ENDDO
851 0 : DEALLOCATE(ax1, ax2, STAT=status)
852 0 : IF(status/=0) THEN
853 0 : CALL HCO_ERROR('AXIS DEALLOCATION ERROR', RC )
854 0 : RETURN
855 : ENDIF
856 0 : NULLIFY(ax1)
857 0 : NULLIFY(ax2)
858 :
859 : END SUBROUTINE AXIS_DELETE
860 : !EOC
861 : !------------------------------------------------------------------------------
862 : ! Harmonized Emissions Component (HEMCO) !
863 : !------------------------------------------------------------------------------
864 : !BOP
865 : !
866 : ! !IROUTINE: Hco2Messy
867 : !
868 : ! !DESCRIPTION: Subroutine HCO2MESSY converts a HEMCO data array into a
869 : ! messy array-structure.
870 : !\\
871 : !\\
872 : ! !INTERFACE:
873 : !
874 0 : SUBROUTINE HCO2MESSY( HcoState, InArr, narr, ax, RC )
875 : !
876 : ! !USES:
877 : !
878 : !
879 : ! !INPUT/OUTPUT PARAMETERS:
880 : !
881 : TYPE(HCO_State), POINTER :: HcoState
882 : REAL(sp), POINTER :: InArr(:,:,:,:)
883 : TYPE(narray), POINTER :: narr(:)
884 : TYPE(axis), POINTER :: ax(:)
885 : INTEGER, INTENT(INOUT) :: RC
886 : !
887 : ! !REVISION HISTORY:
888 : ! 27 Jun 2014 - C. Keller - Initial version
889 : ! See https://github.com/geoschem/hemco for complete history
890 : !EOP
891 : !------------------------------------------------------------------------------
892 : !BOC
893 : !
894 : ! !ROUTINE ARGUMENTS:
895 : !
896 : INTEGER :: NCELLS, NX, NY, NZ, NT
897 : INTEGER :: TMP
898 : INTEGER :: status
899 : INTEGER :: J, L, T, LOW, UPP
900 : CHARACTER(LEN=255) :: MSG, LOC
901 :
902 : !=================================================================
903 : ! HCO2MESSY begins here
904 : !=================================================================
905 :
906 : ! For error handling
907 0 : LOC = 'HCO2MESSY (HCOIO_MESSY_MOD.F90)'
908 :
909 : ! ----------------------------------------------------------------
910 : ! Number of grid cells
911 : ! ----------------------------------------------------------------
912 0 : NX = SIZE(InArr,1)
913 0 : NY = SIZE(InArr,2)
914 0 : NZ = SIZE(InArr,3)
915 0 : NT = SIZE(InArr,4)
916 0 : NCELLS = NX * NY * NZ
917 :
918 : ! create
919 0 : ALLOCATE(narr(NT),STAT=status)
920 0 : IF(status/=0) THEN
921 0 : CALL HCO_ERROR( 'narr allocation error', RC, THISLOC=LOC )
922 0 : RETURN
923 : ENDIF
924 :
925 : ! Set MESSy vector dimensions. Copy from source axis object.
926 : ! Array dimensions are always one less than axis interface dimensions!
927 : ! For last ('free') dimension, set dimension length to 1.
928 0 : DO T = 1, NT
929 0 : narr(T)%n = size(ax)
930 0 : ALLOCATE(narr(T)%dim(narr(T)%n),STAT=status)
931 : IF(status/=0) THEN
932 0 : CALL HCO_ERROR( 'Cannot allocate array dims', RC, THISLOC=LOC )
933 0 : RETURN
934 : ENDIF
935 :
936 0 : DO J = 1, narr(T)%n
937 0 : IF ( ax(J)%dat%n > 0 ) THEN
938 0 : tmp = ax(J)%dat%dim(1)-1
939 : ELSE
940 : tmp = 1
941 : ENDIF
942 0 : narr(T)%dim(J) = tmp
943 : ENDDO
944 :
945 0 : ALLOCATE(narr(T)%vd(NCELLS),STAT=status)
946 0 : IF(status/=0) THEN
947 0 : CALL HCO_ERROR( 'Cannot allocate array', RC, THISLOC=LOC )
948 0 : RETURN
949 : ENDIF
950 : ENDDO !T
951 :
952 : ! ----------------------------------------------------------------
953 : ! Pass HEMCO data array in slices to MESSy data vector.
954 : ! The MESSy data vector is an 'unfolded' HEMCO data array, with
955 : ! the longitude dimension changing the fastest, then followed by
956 : ! latitude, and level. Hence, we pass the HEMCO array in slices
957 : ! along the longitude axis to the MESSy vector, starting with the
958 : ! slice at lat/lev position 1/1, followed by 2/1, 3/1, etc.
959 : ! For now, the MESSy vector is always of type double. In future,
960 : ! we may want to set it to the same type as the HEMCO array and
961 : ! use pointers to the HEMCO data slices!
962 : ! ----------------------------------------------------------------
963 :
964 : ! Vector index counter
965 : LOW = 1
966 0 : UPP = 1
967 :
968 : ! Loop over all higher dimensions. Don't loop over lon because we pass all
969 : ! lon-values in slices.
970 0 : DO L = 1, NZ ! NZ is 1 for 2D arrays
971 0 : DO J = 1, NY
972 0 : UPP = LOW + NX - 1 ! Upper index of slice to be filled
973 0 : DO T = 1, NT
974 0 : narr(T)%vd(LOW:UPP) = InArr(:,J,L,T) ! Pass this slice to vector
975 : ENDDO
976 0 : LOW = UPP + 1 ! Next slice begins right after this one
977 : ENDDO !J
978 : ENDDO !L
979 :
980 : ! Return w/ success
981 0 : RC = HCO_SUCCESS
982 :
983 : END SUBROUTINE HCO2MESSY
984 : !EOC
985 : !------------------------------------------------------------------------------
986 : ! Harmonized Emissions Component (HEMCO) !
987 : !------------------------------------------------------------------------------
988 : !BOP
989 : !
990 : ! !IROUTINE: Messy2Hco
991 : !
992 : ! !DESCRIPTION: Subroutine MESSY2HCO converts a MESSy array structure
993 : ! into a HEMCO data array. This is basically the reverse function of
994 : ! HCO2MESSY.
995 : !\\
996 : !\\
997 : ! !INTERFACE:
998 : !
999 0 : SUBROUTINE MESSY2HCO( HcoState, narr, Lct, LEV, Ptr4D, RC )
1000 : !
1001 : ! !USES:
1002 : !
1003 : !
1004 : ! !INPUT/OUTPUT PARAMETERS:
1005 : !
1006 : TYPE(HCO_State), POINTER :: HcoState
1007 : TYPE(narray), POINTER :: narr(:)
1008 : TYPE(ListCont), POINTER :: Lct
1009 : INTEGER, INTENT(IN ) :: LEV
1010 : REAL(sp), POINTER :: Ptr4D(:,:,:,:)
1011 : INTEGER, INTENT(INOUT) :: RC
1012 : !
1013 : ! !REVISION HISTORY:
1014 : ! 27 Jun 2014 - C. Keller - Initial version
1015 : ! See https://github.com/geoschem/hemco for complete history
1016 : !EOP
1017 : !------------------------------------------------------------------------------
1018 : !BOC
1019 : !
1020 : ! !ROUTINE ARGUMENTS:
1021 : !
1022 : INTEGER :: J, L, T
1023 : INTEGER :: LOW, UPP
1024 : INTEGER :: NX, NY, NZ, NT, Z1, Z2
1025 : CHARACTER(LEN=255) :: MSG
1026 :
1027 : !=================================================================
1028 : ! MESSY2HCO begins here
1029 : !=================================================================
1030 :
1031 : ! Grid dimensions
1032 0 : NX = HcoState%NX
1033 0 : NY = HcoState%NY
1034 0 : NT = SIZE(narr)
1035 :
1036 : ! Vertical levels to be filled on output array. Only level LEV
1037 : ! is filled if NC is different than one!
1038 0 : IF ( Lct%Dct%Dta%SpaceDim == 2 ) THEN
1039 : Z1 = 1
1040 : Z2 = 1
1041 0 : ELSEIF ( Lct%Dct%Dta%SpaceDim == 3 ) THEN
1042 0 : IF ( ASSOCIATED(Ptr4D) ) THEN
1043 0 : Z1 = LEV
1044 0 : Z2 = LEV
1045 : ELSE
1046 0 : Z1 = 1
1047 0 : Z2 = SIZE(Lct%Dct%Dta%V3(1)%Val,3)
1048 : ENDIF
1049 : ENDIF
1050 :
1051 : ! Vector index counter
1052 0 : LOW = 1
1053 0 : UPP = 1
1054 :
1055 : ! If temporary array Ptr4D is provided, make sure that its
1056 : ! dimensions are correct
1057 0 : IF ( ASSOCIATED(Ptr4D) ) THEN
1058 0 : NZ = Z2 - Z1 + 1
1059 : IF ( ( SIZE(Ptr4D,1) /= NX ) .OR. &
1060 : ( SIZE(Ptr4D,2) /= NY ) .OR. &
1061 0 : ( SIZE(Ptr4D,3) /= NZ ) .OR. &
1062 : ( SIZE(Ptr4D,4) /= NT ) ) THEN
1063 0 : WRITE(MSG,*) 'Temporary pointer has wrong dimensions: ', &
1064 0 : TRIM(Lct%Dct%cName), NX, NY, NZ, NT, SIZE(Ptr4D,1), SIZE(Ptr4D,2), SIZE(Ptr4D,3), SIZE(Ptr4D,4)
1065 : CALL HCO_ERROR ( MSG, RC, &
1066 0 : THISLOC='MESSY2HCO (hcoio_messy_mod.F90)' )
1067 0 : RETURN
1068 : ENDIF
1069 : ENDIF
1070 :
1071 : ! Loop over all higher dimensions
1072 0 : DO L = Z1, Z2
1073 0 : DO J = 1, NY
1074 :
1075 : ! Upper index of slice to be filled
1076 0 : UPP = LOW + NX - 1
1077 :
1078 : ! Do for every time slice
1079 0 : DO T = 1, NT
1080 :
1081 : ! If provided, write into temporary array
1082 0 : IF ( ASSOCIATED(Ptr4D) ) THEN
1083 0 : Ptr4D(:,J,L,T) = narr(T)%vd(LOW:UPP)
1084 :
1085 : ! If temporary array Ptr4D is not provided, write directly
1086 : ! into list container array
1087 : ELSE
1088 :
1089 : ! 2D
1090 0 : IF ( Lct%Dct%Dta%SpaceDim == 2 ) THEN
1091 0 : Lct%Dct%Dta%V2(T)%Val(:,J) = narr(T)%vd(LOW:UPP)
1092 : ! 3D
1093 : ELSE
1094 0 : Lct%Dct%Dta%V3(T)%Val(:,J,L) = narr(T)%vd(LOW:UPP)
1095 : ENDIF
1096 :
1097 : ENDIF
1098 : ENDDO !NT
1099 :
1100 : ! Next slice begins right after this one
1101 0 : LOW = UPP + 1
1102 : ENDDO !J
1103 : ENDDO !L
1104 :
1105 : ! Return w/ success
1106 0 : RC = HCO_SUCCESS
1107 :
1108 : END SUBROUTINE MESSY2HCO
1109 : !EOC
1110 : END MODULE HCOIO_MESSY_MOD
|