Line data Source code
1 : !BOC
2 : #if defined ( MODEL_GCCLASSIC ) || defined( MODEL_WRF ) || defined( MODEL_CESM ) || defined( HEMCO_STANDALONE )
3 : ! The 'standard' HEMCO I/O module is used for:
4 : ! - HEMCO Standalone (HEMCO_STANDALONE)
5 : ! - GEOS-Chem 'Classic' (MODEL_GCCLASSIC)
6 : ! - WRF-GC (MODEL_WRF)
7 : ! - CESM-GC and CAM-Chem / HEMCO-CESM (MODEL_CESM)
8 : !EOC
9 : !------------------------------------------------------------------------------
10 : ! Harmonized Emissions Component (HEMCO) !
11 : !------------------------------------------------------------------------------
12 : !BOP
13 : !
14 : ! !MODULE: hcoio_write_mod.F90
15 : !
16 : ! !DESCRIPTION: Module HCOIO\_write\_mod.F90 is the HEMCO data output
17 : ! interface for the 'standard' model environment. It contains routines to
18 : ! write out diagnostics into a netCDF file.
19 : !\\
20 : !\\
21 : ! !INTERFACE:
22 : !
23 : MODULE HCOIO_Write_Mod
24 : !
25 : ! !USES:
26 : !
27 : USE HCO_ERROR_MOD
28 : USE HCO_DIAGN_MOD
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 : !
33 : ! !PUBLIC MEMBER FUNCTIONS:
34 : !
35 : PUBLIC :: HCOIO_Write
36 : !
37 : ! !PRIVATE MEMBER FUNCTIONS:
38 : !
39 : PRIVATE :: ConstructTimeStamp
40 : !
41 : ! !REMARKS:
42 : ! HEMCO diagnostics are still in testing mode. We will fully activate them
43 : ! at a later time. They will be turned on when debugging & unit testing.
44 : !
45 : ! !REVISION HISTORY:
46 : ! 04 May 2014 - C. Keller - Initial version
47 : ! See https://github.com/geoschem/hemco for complete history
48 : !EOP
49 : !------------------------------------------------------------------------------
50 : !BOC
51 : !
52 : ! !DEFINED PARAMETERS:
53 : !
54 : ! Fill value used in HEMCO diagnostics netCDF files.
55 : ! REAL(hp), PARAMETER :: FillValue = 1.e-31_hp
56 : REAL(sp), PARAMETER :: FillValue = HCO_MISSVAL
57 :
58 : CONTAINS
59 : !EOC
60 : !------------------------------------------------------------------------------
61 : ! Harmonized Emissions Component (HEMCO) !
62 : !------------------------------------------------------------------------------
63 : !BOP
64 : !
65 : ! !IROUTINE: HCOIO_write_std
66 : !
67 : ! !DESCRIPTION: Subroutine HCOIO\_write\_std writes diagnostics to
68 : ! netCDF file. If the ForceWrite flag is set to TRUE, all diagnostics are
69 : ! written out except they have already been written out during this time
70 : ! step. This option is usually only used at the end of a simulation run.
71 : ! If ForceWrite is False, only the diagnostics that are at the end of their
72 : ! time averaging interval are written. For example, if the current month
73 : ! is different from the previous (emissions) month, all diagnostics with
74 : ! hourly, daily and monthly time averaging intervals are written out.
75 : ! If the optional argument OnlyIfFirst is set to TRUE, diagnostics will
76 : ! only be written out if its nnGetCalls is 1. This can be used to avoid
77 : ! that diagnostics will be written out twice. The nnGetCalls is reset to
78 : ! zero the first time a diagnostics is updated. For diagnostics that
79 : ! point to data stored somewhere else (i.e. that simply contain a data
80 : ! pointer, nnGetCalls is never reset and keeps counting.
81 : !\\
82 : !\\
83 : ! !INTERFACE:
84 : !
85 0 : SUBROUTINE HCOIO_Write( HcoState, ForceWrite, &
86 : RC, PREFIX, UsePrevTime, &
87 : OnlyIfFirst, COL )
88 : !
89 : ! !USES:
90 : !
91 : USE HCO_m_netCDF_io_define
92 : USE HCO_m_netcdf_io_read
93 : USE HCO_m_netcdf_io_open
94 : USE HCO_Ncdf_Mod, ONLY : NC_Open
95 : USE HCO_Ncdf_Mod, ONLY : NC_Read_Time
96 : USE HCO_Ncdf_Mod, ONLY : NC_Read_Arr
97 : USE HCO_Ncdf_Mod, ONLY : NC_Create
98 : USE HCO_Ncdf_Mod, ONLY : NC_Close
99 : USE HCO_Ncdf_Mod, ONLY : NC_Var_Def
100 : USE HCO_Ncdf_Mod, ONLY : NC_Var_Write
101 : USE HCO_Ncdf_Mod, ONLY : NC_Get_RefDateTime
102 : USE HCO_CHARPAK_Mod, ONLY : TRANLC
103 : USE HCO_Chartools_Mod, ONLY : HCO_CharParse
104 : USE HCO_State_Mod, ONLY : HCO_State
105 : USE HCO_JulDay_Mod, ONLY : JulDay
106 : USE HCO_EXTLIST_MOD, ONLY : GetExtOpt, CoreNr
107 : USE HCO_Types_Mod, ONLY : DiagnCont
108 : USE HCO_Clock_Mod
109 :
110 : ! Parameters for netCDF routines
111 : include "netcdf.inc"
112 : !
113 : ! !INPUT PARAMETERS:
114 : !
115 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state object
116 : LOGICAL, INTENT(IN ) :: ForceWrite ! Write all diagnostics?
117 : CHARACTER(LEN=*), OPTIONAL, INTENT(IN ) :: PREFIX ! File prefix
118 : LOGICAL, OPTIONAL, INTENT(IN ) :: UsePrevTime ! Use previous time
119 : LOGICAL, OPTIONAL, INTENT(IN ) :: OnlyIfFirst ! Only write if nnDiagn is 1
120 : INTEGER, OPTIONAL, INTENT(IN ) :: COL ! Collection Nr.
121 : !
122 : ! !INPUT/OUTPUT PARAMETERS:
123 : !
124 :
125 : INTEGER, INTENT(INOUT) :: RC ! Failure or success
126 : !
127 : ! !REVISION HISTORY:
128 : ! 12 Sep 2013 - C. Keller - Initial version
129 : ! See https://github.com/geoschem/hemco for complete history
130 : !EOP
131 : !------------------------------------------------------------------------------
132 : !BOC
133 : !
134 : ! !LOCAL VARIABLES:
135 : !
136 : INTEGER :: I, PS, CNT, levIdTmp, indexL, indexR
137 : REAL(dp) :: GMT, JD1, JD1985, JD_DELTA, THISDAY, P0
138 : REAL(sp) :: TMP, JD_DELTA_RND
139 : INTEGER :: YYYY, MM, DD, h, m, s
140 0 : REAL(sp), POINTER :: nctime(:)
141 0 : REAL(dp), POINTER :: Arr1D(:)
142 0 : INTEGER, POINTER :: Int1D(:)
143 0 : REAL(sp), POINTER :: Arr3D(:,:,:)
144 0 : REAL(sp), POINTER :: Arr4D(:,:,:,:)
145 0 : REAL(sp), POINTER :: Arr4DOld(:,:,:,:)
146 0 : REAL*8, POINTER :: timeVec(:)
147 0 : REAL(hp), POINTER :: hyam(:)
148 0 : REAL(hp), POINTER :: hybm(:)
149 : TYPE(DiagnCont), POINTER :: ThisDiagn
150 : INTEGER :: FLAG
151 : CHARACTER(LEN=255) :: ncFile
152 : CHARACTER(LEN=255) :: Pfx, title, Reference, Contact
153 : CHARACTER(LEN=255) :: myLName, mySName, myFterm
154 : CHARACTER(LEN=255) :: MSG
155 : CHARACTER(LEN=255) :: RefTime
156 : CHARACTER(LEN=4 ) :: Yrs
157 : CHARACTER(LEN=2 ) :: Mts, Dys, hrs, mns
158 : CHARACTER(LEN=31) :: myName, myUnit, OutOper
159 : CHARACTER(LEN=63) :: timeunit
160 : INTEGER :: fId, lonId, latId, levId, TimeId
161 : INTEGER :: VarCt
162 : INTEGER :: nLon, nLat, nLev, nLevTmp, nTime
163 : INTEGER :: Prc, L
164 : INTEGER :: lymd, lhms
165 : INTEGER :: refYYYY, refMM, refDD, refh, refm, refs
166 : LOGICAL :: EOI, DoWrite, PrevTime, FOUND
167 : LOGICAL :: NoLevDim, DefMode
168 : LOGICAL :: IsOldFile
169 :
170 : CHARACTER(LEN=255), PARAMETER :: LOC = 'HCOIO_WRITE_STD (hcoio_write_std_mod.F90)'
171 :
172 : !=================================================================
173 : ! HCOIO_WRITE_STD begins here!
174 : !=================================================================
175 :
176 : ! Init
177 0 : RC = HCO_SUCCESS
178 0 : CNT = 0
179 0 : Arr1D => NULL()
180 0 : Int1D => NULL()
181 0 : Arr3D => NULL()
182 0 : Arr4D => NULL()
183 0 : Arr4DOld => NULL()
184 0 : timeVec => NULL()
185 0 : nctime => NULL()
186 0 : ThisDiagn => NULL()
187 :
188 : ! Collection number
189 0 : PS = HcoState%Diagn%HcoDiagnIDDefault
190 0 : IF ( PRESENT(COL) ) PS = COL
191 :
192 : ! Check if it's time to write out this collection. Also set the
193 : ! end-of-interval EOI flag accordingly. This will be used lateron
194 : ! when calling Diagn_Get. Since all diagnostic containers in a
195 : ! given collection have the same output frequency, this is somewhat
196 : ! redundant (because we already check here if it is time to write
197 : ! out this particular collection). Keep it here for backwards
198 : ! consistency (ckeller, 8/6/2015).
199 0 : IF ( ForceWrite ) THEN
200 0 : DoWrite = .TRUE.
201 0 : EOI = .FALSE.
202 : ELSE
203 0 : DoWrite = DiagnCollection_IsTimeToWrite( HcoState, PS )
204 0 : EOI = .TRUE.
205 : ENDIF
206 :
207 : ! Create current time stamps (to be used to archive time stamps)
208 : CALL HcoClock_Get( HcoState%Clock,sYYYY=YYYY,sMM=MM,&
209 0 : sDD=DD,sH=h,sM=m,sS=s,RC=RC)
210 0 : IF ( RC /= HCO_SUCCESS ) THEN
211 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
212 0 : RETURN
213 : ENDIF
214 0 : lymd = YYYY*10000 + MM*100 + DD
215 0 : lhms = h *10000 + m *100 + s
216 :
217 : ! Leave here if it's not time to write diagnostics. On the first
218 : ! time step, set lastYMD and LastHMS to current dates.
219 0 : IF ( .NOT. DoWrite ) THEN
220 0 : IF ( .NOT. DiagnCollection_LastTimesSet(HcoState%Diagn,PS) ) THEN
221 : CALL DiagnCollection_Set ( HcoState%Diagn, COL=PS, &
222 0 : LastYMD=lymd, LastHMS=lhms, RC=RC )
223 : ENDIF
224 0 : RETURN
225 : ENDIF
226 :
227 : ! Inherit precision from HEMCO
228 0 : Prc = HP
229 :
230 : ! Get PrevTime flag from input argument or set to default (=> TRUE)
231 0 : IF ( PRESENT(UsePrevTime) ) THEN
232 0 : PrevTime = UsePrevTime
233 : ELSE
234 0 : PrevTime = .TRUE.
235 : ENDIF
236 :
237 : !-----------------------------------------------------------------
238 : ! Don't define level dimension if there are no 3D fields to write
239 : ! This is an optional feature. By default, all diagnostics have
240 : ! the full dimension definitions (lon,lat,lev,time) even if all
241 : ! output fields are only 2D. If the flag DiagnNoLevDim is
242 : ! enabled, the lev dimension is not defined if there are no 3D
243 : ! fields on the file.
244 : !-----------------------------------------------------------------
245 0 : NoLevDim = .FALSE.
246 : CALL GetExtOpt ( HcoState%Config, CoreNr, 'DiagnNoLevDim', &
247 0 : OptValBool=NoLevDim, Found=Found, RC=RC )
248 0 : IF ( RC /= HCO_SUCCESS ) THEN
249 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
250 0 : RETURN
251 : ENDIF
252 0 : IF ( Found ) THEN
253 0 : IF ( NoLevDim ) THEN
254 :
255 : ! Loop over all diagnostics to see if any is 3D
256 0 : ThisDiagn => NULL()
257 : DO WHILE ( .TRUE. )
258 :
259 : ! Get next diagnostics in list. This will return the next
260 : ! diagnostics container that contains content.
261 : CALL Diagn_Get ( HcoState, EOI, &
262 0 : ThisDiagn, FLAG, RC, COL=PS )
263 0 : IF ( RC /= HCO_SUCCESS ) THEN
264 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
265 0 : RETURN
266 : ENDIF
267 0 : IF ( FLAG /= HCO_SUCCESS ) EXIT
268 :
269 : ! If this is a 3D diagnostics, we must write the level
270 : ! coordinate
271 0 : IF ( ThisDiagn%SpaceDim == 3 ) THEN
272 0 : NoLevDim = .FALSE.
273 0 : EXIT
274 : ENDIF
275 : ENDDO
276 : ENDIF
277 : ENDIF
278 :
279 : !-----------------------------------------------------------------
280 : ! Create output file
281 : !-----------------------------------------------------------------
282 :
283 : ! Define grid dimensions
284 0 : nLon = HcoState%NX
285 0 : nLat = HcoState%NY
286 0 : nLev = HcoState%NZ
287 0 : nTime = 1
288 :
289 : ! Initialize mirror variables
290 0 : allocate(Arr4D(nlon,nlat,nlev,ntime))
291 0 : allocate(Arr3D(nlon,nlat,ntime))
292 0 : Arr3D = 0.0_sp
293 0 : Arr4D = 0.0_sp
294 :
295 : ! Construct filename: diagnostics will be written into file
296 : ! PREFIX.YYYYMMDDhm.nc, where PREFIX is the input argument or
297 : ! (if not present) obtained from the HEMCO configuration file.
298 : CALL ConstructTimeStamp ( HcoState, PS, PrevTime, &
299 0 : YYYY, MM, DD, h, m, RC )
300 0 : IF ( RC /= HCO_SUCCESS ) THEN
301 0 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
302 0 : RETURN
303 : ENDIF
304 :
305 : ! Write datetime
306 0 : WRITE( Yrs, '(i4.4)' ) YYYY
307 0 : WRITE( Mts, '(i2.2)' ) MM
308 0 : WRITE( Dys, '(i2.2)' ) DD
309 0 : WRITE( hrs, '(i2.2)' ) h
310 0 : WRITE( mns, '(i2.2)' ) m
311 :
312 : ! Get prefix
313 0 : IF ( PRESENT(PREFIX) ) THEN
314 0 : Pfx = PREFIX
315 : ELSE
316 0 : CALL DiagnCollection_Get( HcoState%Diagn, PS, PREFIX=Pfx, RC=RC )
317 0 : IF ( RC /= HCO_SUCCESS ) THEN
318 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
319 0 : RETURN
320 : ENDIF
321 : ENDIF
322 0 : ncFile = TRIM(Pfx)//'.'//Yrs//Mts//Dys//hrs//mns//'.nc'
323 :
324 : ! Place HEMCO restart files in the Restarts folder of the run directory
325 0 : IF ( PS == HcoState%Diagn%HcoDiagnIDRestart ) THEN
326 0 : ncFile = 'Restarts/' // TRIM( ncFile )
327 : ENDIF
328 :
329 : ! Multiple time slice update. Comment out for now since it causes
330 : ! timestamping the filename twice (ewl, 10/19/18)
331 : ! Add default time stamp if no time tokens are in the file template.
332 : ! This also ensures backward compatibility.
333 : !IF ( INDEX(TRIM(ncFile),'$') <= 0 ) THEN
334 : ! ncFile = TRIM(ncFile)//'.$YYYY$MM$DD$HH$MN.nc'
335 : !ENDIF
336 : !CALL HCO_CharParse ( HcoState%Config, ncFile, YYYY, MM, DD, h, m, RC )
337 : !IF ( RC /= HCO_SUCCESS ) RETURN
338 :
339 : ! Use filename prefix for title, replacing '_' with spaces
340 : ! NOTE: Prefix can only contain up to two underscores
341 0 : indexL = SCAN( Pfx, '_', .FALSE. ) ! Return left-most position
342 0 : indexR = SCAN( Pfx, '_', .TRUE. ) ! Return right-most position
343 0 : IF ( indexL > 0 .AND. indexR > 0 ) THEN
344 : title = Pfx(1:indexL-1) // ' ' // &
345 : Pfx(indexL+1:indexR-1) // ' ' // &
346 0 : Pfx(indexR+1:)
347 0 : ELSE IF ( indexL > 0 .AND. indexR == 0 ) THEN
348 0 : title = Pfx(1:indexL-1) // ' ' // Pfx(indexL+1:)
349 : ELSE
350 0 : title = Pfx
351 : ENDIF
352 :
353 : ! verbose
354 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) .AND. PS==1 ) THEN
355 0 : MSG = 'Write diagnostics into file '//TRIM(ncFile)
356 0 : CALL HCO_MSG( HcoState%Config%Err, MSG )
357 : ENDIF
358 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) .AND. PS==1 ) THEN
359 0 : WRITE(MSG,*) '--> write level dimension: ', .NOT.NoLevDim
360 0 : CALL HCO_MSG( HcoState%Config%Err, MSG )
361 : ENDIF
362 :
363 : ! Check if file already exists. If so, add new diagnostics to this file
364 : ! (instead of creating a new one)
365 0 : INQUIRE( FILE=ncFile, EXIST=IsOldFile )
366 :
367 : ! Disable multiple time slice update since causes an issue writing
368 : ! restart files. Re-enable when restart files are written via HISTORY
369 : ! rather than HEMCO by deleting the forcing of IsOldFile below.
370 : ! (ewl, 10/19/18)
371 0 : IsOldFile = .FALSE.
372 :
373 : ! If file exists, open file and get time dimension
374 : IF ( IsOldFile ) THEN
375 : CALL Ncop_Wr( fID, ncFile )
376 : CALL NC_READ_TIME( fID, ntime, timeunit, timeVec, RC=RC )
377 :
378 : ! new file will have one more time dimension
379 : ntime = ntime + 1
380 :
381 : ! Create output file
382 : ELSE
383 :
384 : ! Define a variable for the number of levels, which will either be -1
385 : ! (if all 2D data) or the number of levels in the grid (for 3D data).
386 0 : IF ( NoLevDim ) THEN
387 0 : nLevTmp = -1
388 : ELSE
389 0 : nLevTmp = nLev
390 : ENDIF
391 :
392 : ! Define extra metadata for global attributes
393 0 : Reference = 'http://wiki.geos-chem.org/The_HEMCO_Users_Guide'
394 0 : Contact = 'GEOS-Chem Support Team (geos-chem-support@as.harvard.edu)'
395 :
396 : ! Create output file
397 : ! Pass CREATE_NC4 to make file format netCDF-4 (mps, 3/3/16)
398 : ! Now create netCDF file with time dimension as UNLIMITED (bmy, 3/8/17)
399 : CALL NC_Create( NcFile = NcFile, &
400 : Title = Title, &
401 : Reference = Reference, &
402 : Contact = Contact, &
403 : nLon = nLon, &
404 : nLat = nLat, &
405 : nLev = nLevTmp, &
406 : nTime = NF_UNLIMITED, &
407 : fId = fId, &
408 : lonId = lonId, &
409 : latId = latId, &
410 : levId = levId, &
411 : timeId = timeId, &
412 : VarCt = VarCt, &
413 0 : CREATE_NC4 =.TRUE. )
414 :
415 : ENDIF
416 :
417 : !-----------------------------------------------------------------
418 : ! Write grid dimensions (incl. time)
419 : !-----------------------------------------------------------------
420 0 : IF ( .NOT. IsOldFile ) THEN
421 :
422 : ! Write longitude axis variable ("lon") to file
423 : CALL NC_Var_Def( fId = fId, &
424 : lonId = lonId, &
425 : latId = -1, &
426 : levId = -1, &
427 : timeId = -1, &
428 : VarName = 'lon', &
429 : VarLongName = 'Longitude', &
430 : VarUnit = 'degrees_east', &
431 : Axis = 'X', &
432 : DataType = dp, &
433 : VarCt = VarCt, &
434 0 : Compress = .TRUE. )
435 0 : ALLOCATE( Arr1D( nLon ) )
436 0 : Arr1D = HcoState%Grid%XMID%Val(:,1)
437 0 : CALL NC_Var_Write( fId, 'lon', Arr1D=Arr1D )
438 0 : DEALLOCATE( Arr1D )
439 :
440 : ! Write latitude axis variable ("lat") to file
441 : CALL NC_Var_Def( fId = fId, &
442 : lonId = -1, &
443 : latId = latId, &
444 : levId = -1, &
445 : timeId = -1, &
446 : VarName = 'lat', &
447 : VarLongName = 'Latitude', &
448 : VarUnit = 'degrees_north', &
449 : Axis = 'Y', &
450 : DataType = dp, &
451 : VarCt = VarCt, &
452 0 : Compress = .TRUE. )
453 0 : ALLOCATE( Arr1D( nLat ) )
454 0 : Arr1D = HcoState%Grid%YMID%Val(1,:)
455 0 : CALL NC_Var_Write( fId, 'lat', Arr1D=Arr1D )
456 0 : DEALLOCATE( Arr1D )
457 :
458 : ! Write vertical grid parameters to file (if necessary)
459 0 : IF ( .NOT. NoLevDim ) THEN
460 :
461 : ! Reference pressure [Pa]
462 0 : P0 = 1.0e+05_dp
463 :
464 : ! Allocate vertical coordinate arrays
465 0 : ALLOCATE( Arr1D( nLev ) )
466 0 : ALLOCATE( hyam ( nLev ) )
467 0 : ALLOCATE( hybm ( nLev ) )
468 :
469 : ! Construct vertical level coordinates
470 0 : DO L = 1, nLev
471 :
472 : ! A parameter at grid midpoints
473 0 : hyam(L) = ( HcoState%Grid%zGrid%Ap(L) &
474 0 : + HcoState%Grid%zGrid%Ap(L+1) ) * 0.5_dp
475 :
476 : ! B parameter at grid midpoints
477 0 : hybm(L) = ( HcoState%Grid%zGrid%Bp(L) &
478 0 : + HcoState%Grid%zGrid%Bp(L+1) ) * 0.5_dp
479 :
480 : ! Vertical level coordinate
481 0 : Arr1d(L) = ( hyam(L) / P0 ) + hybm(L)
482 :
483 : ENDDO
484 :
485 : ! Write level axis variable ("lev") to file
486 : ! Define extra metadata for calls to NC_Var_Def
487 0 : myLName = 'hybrid level at midpoints ((A/P0)+B)'
488 0 : mySName = 'atmosphere_hybrid_sigma_pressure_coordinate'
489 0 : myFTerm = 'a: hyai b: hybi p0: P0 ps: PS'
490 : CALL NC_Var_Def( fId = fId, &
491 : lonId = -1, &
492 : latId = -1, &
493 : levId = levId, &
494 : timeId = -1, &
495 : VarName = 'lev', &
496 : VarLongName = MyLName, &
497 : StandardName = MySName, &
498 : FormulaTerms = myFTerm, &
499 : VarUnit = 'level', &
500 : Axis = 'Z', &
501 : Positive = 'up', &
502 : DataType = dp, &
503 : VarCt = VarCt, &
504 0 : Compress = .TRUE. )
505 0 : CALL NC_Var_Write( fId, 'lev', Arr1D=Arr1D )
506 :
507 : ! Write hybrid A coordinate ("hyam") to file
508 : ! Define extra metadata for calls to NC_Var_Def
509 0 : myLName = 'hybrid A coefficient at layer midpoints'
510 : CALL NC_Var_Def( fId = fId, &
511 : lonId = -1, &
512 : latId = -1, &
513 : levId = levId, &
514 : timeId = -1, &
515 : VarName = 'hyam', &
516 : VarLongName = MyLName, &
517 : VarUnit = 'Pa', &
518 : DataType = dp, &
519 : VarCt = VarCt, &
520 0 : Compress = .TRUE. )
521 0 : CALL NC_Var_Write ( fId, 'hyam', Arr1D=hyam )
522 :
523 : ! Write hybrid B coordinate ("hybm") to file
524 : ! Define extra metadata for calls to NC_Var_Def
525 0 : myLName = 'hybrid B coefficient at layer midpoints'
526 : CALL NC_Var_Def( fId = fId, &
527 : lonId = -1, &
528 : latId = -1, &
529 : levId = levId, &
530 : timeId = -1, &
531 : VarName = 'hybm', &
532 : VarLongName = MyLName, &
533 : VarUnit = '1', &
534 : DataType = dp, &
535 : VarCt = VarCt, &
536 0 : Compress = .TRUE. )
537 0 : CALL NC_Var_Write( fId, 'hybm', Arr1D=hybm )
538 :
539 : ! Write out reference pressure (P0) to file
540 : CALL NC_Var_Def( fId = fId, &
541 : lonId = -1, &
542 : latId = -1, &
543 : levId = -1, &
544 : timeId = -1, &
545 : VarName = 'P0', &
546 : VarLongName = 'Reference pressure', &
547 : VarUnit = 'Pa', &
548 : DataType = dp, &
549 : VarCt = VarCt, &
550 0 : Compress = .FALSE. )
551 0 : CALL NC_Var_Write( fId, 'P0', P0 )
552 :
553 : ! Deallocate arrays
554 0 : DEALLOCATE( Arr1d )
555 0 : DEALLOCATE( hyam )
556 0 : DEALLOCATE( hybm )
557 :
558 : ENDIF
559 : ENDIF
560 :
561 : !------------------------------------------------------------------------
562 : ! Write time axis variable ("time") to file
563 : !------------------------------------------------------------------------
564 :
565 : ! JD1 is the julian day of the data slice
566 0 : GMT = REAL(h,dp) + (REAL(m,dp)/60.0_dp) + (REAL(s,dp)/3600.0_dp)
567 0 : THISDAY = DD + ( GMT / 24.0_dp )
568 0 : JD1 = JULDAY ( YYYY, MM, THISDAY )
569 :
570 : ! Check if reference time is given in HEMCO configuration file
571 : CALL GetExtOpt ( HcoState%Config, CoreNr, 'DiagnRefTime', &
572 0 : OptValChar=RefTime, Found=Found, RC=RC )
573 0 : IF ( RC /= HCO_SUCCESS ) THEN
574 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
575 0 : RETURN
576 : ENDIF
577 :
578 : ! Use specified reference time (if available)
579 0 : IF ( Found ) THEN
580 0 : timeunit = ADJUSTL(TRIM(RefTime))
581 0 : CALL TRANLC( timeunit )
582 : CALL NC_GET_REFDATETIME( timeunit, refYYYY, refMM, refDD, refh, &
583 0 : refm, refs, RC )
584 0 : refs = 0
585 0 : IF ( RC /= HCO_SUCCESS ) THEN
586 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
587 0 : RETURN
588 : ENDIF
589 : GMT = REAL(MAX(refh,0),dp) + (REAL(MAX(refm,0),dp)/60.0_dp) + &
590 0 : (REAL(MAX(refs,0),dp)/3600.0_dp)
591 0 : THISDAY = refDD + ( GMT / 24.0_dp )
592 0 : JD1985 = JULDAY ( refYYYY, refMM, THISDAY )
593 :
594 : ! Use current time if not found
595 : ELSE
596 0 : WRITE(timeunit,100) YYYY,MM,DD,h,m,s
597 0 : JD1985 = JD1
598 : ENDIF
599 : 100 FORMAT ( 'hours since ',i4.4,'-',i2.2,'-',i2.2,' ',i2.2,':',i2.2,':',i2.2,' GMT' )
600 :
601 : ! Calculate time value
602 0 : JD_DELTA = (JD1 - JD1985 )
603 :
604 : ! Default is 'days since'. Adjust for 'hours since', 'minutes since',
605 : ! 'seconds since'.
606 0 : IF ( timeunit(1:4) == 'days' ) THEN
607 : ! all ok
608 0 : ELSEIF ( timeunit(1:5) == 'hours' ) THEN
609 0 : JD_DELTA = JD_DELTA * 24.0_dp
610 0 : ELSEIF ( timeunit(1:7) == 'minutes' ) THEN
611 0 : JD_DELTA = JD_DELTA * 24.0_dp * 60.0_dp
612 0 : ELSEIF ( timeunit(1:7) == 'seconds' ) THEN
613 0 : JD_DELTA = JD_DELTA * 24.0_dp * 3600.0_dp
614 : ELSE
615 : MSG = 'Unrecognized output reference time, will ' // &
616 0 : 'assume `days since`: '//TRIM(timeunit)
617 0 : CALL HCO_WARNING( HcoState%Config%Err, MSG, RC, THISLOC=LOC )
618 : ENDIF
619 :
620 : ! Special case where we have an old file but it has the same time stamp: in
621 : ! that case simply overwrite the current values
622 : ! Comment out code for single precision rounded time (ewl, 10/18/18)
623 : !IF ( IsOldFile .AND. ntime == 2 .AND. timeVec(1) == JD_DELTA_RND ) THEN
624 0 : IF ( IsOldFile .AND. ntime == 2 ) THEN
625 0 : IF ( timeVec(1) == JD_DELTA ) THEN
626 0 : ntime = 1
627 : ENDIF
628 : ENDIF
629 0 : ALLOCATE( nctime(ntime) )
630 0 : IF ( IsOldFile .AND. ntime > 1 ) THEN
631 0 : nctime(1:ntime-1) = timeVec(:)
632 : ENDIF
633 0 : nctime(ntime) = JD_DELTA
634 :
635 0 : IF ( .NOT. IsOldFile ) THEN
636 : CALL NC_Var_Def( fId = fId, &
637 : lonId = -1, &
638 : latId = -1, &
639 : levId = -1, &
640 : timeId = timeId, &
641 : VarName = 'time', &
642 : VarLongName = 'Time', &
643 : VarUnit = TimeUnit, &
644 : Axis = 'T', &
645 : Calendar = 'gregorian', &
646 : DataType = 8, &
647 : VarCt = VarCt, &
648 0 : Compress = .TRUE. )
649 : ENDIF
650 0 : CALL NC_VAR_WRITE( fId, 'time', Arr1D=nctime )
651 0 : DEALLOCATE( nctime )
652 0 : IF ( ASSOCIATED(timeVec) ) DEALLOCATE( timeVec )
653 :
654 : !-----------------------------------------------------------------
655 : ! Write out grid box areas
656 : !-----------------------------------------------------------------
657 :
658 0 : IF ( .NOT. IsOldFile ) THEN
659 : CALL NC_Var_Def( fId = fId, &
660 : lonId = lonId, &
661 : latId = latId, &
662 : levId = -1, &
663 : timeId = -1, &
664 : VarName = 'AREA', &
665 : VarLongName = 'Grid box area', &
666 : VarUnit = 'm2', &
667 : DataType = Prc, &
668 : VarCt = VarCt, &
669 0 : Compress = .TRUE. )
670 0 : CALL NC_Var_Write ( fId, 'AREA', Arr2D=HcoState%Grid%Area_M2%Val )
671 : ENDIF
672 :
673 : !-----------------------------------------------------------------
674 : ! Write diagnostics
675 : !-----------------------------------------------------------------
676 :
677 : ! Run this section twice, first in define mode for metadata, then in
678 : ! data mode to write variables
679 0 : DO I=1,2
680 :
681 : ! Skip definition mode for existing file
682 0 : IF ( I==1 .AND. IsOldFile ) CYCLE
683 :
684 0 : IF (I==1) THEN
685 : ! Open netCDF define mode
686 0 : CALL NcBegin_Def( fID )
687 0 : DefMode=.TRUE.
688 : ELSE
689 : ! IF ( .NOT. IsOldFile ) THEN
690 : ! Close netCDF define mode
691 0 : CALL NcEnd_Def( fID )
692 : ! ENDIF
693 0 : DefMode=.False.
694 : ENDIF
695 :
696 : ! Loop over all diagnostics in diagnostics list
697 0 : ThisDiagn => NULL()
698 0 : DO WHILE ( .TRUE. )
699 :
700 : ! Get next diagnostics in list. This will return the next
701 : ! diagnostics container that contains content.
702 0 : CALL Diagn_Get ( HcoState, EOI, ThisDiagn, FLAG, RC, COL=PS )
703 0 : IF ( RC /= HCO_SUCCESS ) THEN
704 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
705 0 : RETURN
706 : ENDIF
707 0 : IF ( FLAG /= HCO_SUCCESS ) EXIT
708 :
709 : ! Only write diagnostics if this is the first Diagn_Get call for
710 : ! this container and time step.
711 0 : IF ( PRESENT( OnlyIfFirst ) ) THEN
712 0 : IF ( OnlyIfFirst .AND. ThisDiagn%nnGetCalls > 1 ) CYCLE
713 : ENDIF
714 :
715 : ! Define variable
716 0 : myName = ThisDiagn%cName
717 0 : myUnit = ThisDiagn%OutUnit
718 0 : IF ( ThisDiagn%SpaceDim == 3 ) THEN
719 0 : levIdTmp = levId
720 : ELSE
721 0 : levIdTmp = -1
722 : ENDIF
723 :
724 : ! Error check: this should never happen!
725 0 : IF ( levIdTmp > 0 .AND. NoLevDim ) THEN
726 : MSG = 'Level dimension undefined but 3D container found: ' &
727 0 : // TRIM(myName)
728 0 : CALL HCO_ERROR(MSG,RC,THISLOC=LOC)
729 0 : RETURN
730 : ENDIF
731 :
732 0 : IF (DefMode) THEN
733 :
734 : !------------------------------------
735 : ! Define variables in define mode
736 : !------------------------------------
737 :
738 : ! Define variable as single precision
739 : CALL NC_Var_Def( fId = fId, &
740 : lonId = lonId, &
741 : latId = latId, &
742 : levId = levIdTmp, &
743 : timeId = timeId, &
744 : VarName = TRIM(myName), &
745 : VarLongName = ThisDiagn%long_name, &
746 : VarUnit = TRIM(myUnit), &
747 : AvgMethod = ThisDiagn%AvgName, &
748 : MissingValue = FillValue, &
749 : DataType = sp, &
750 : VarCt = VarCt, &
751 : DefMode = DefMode, &
752 0 : Compress = .True. )
753 :
754 : ELSE
755 :
756 : !------------------------------------
757 : ! Write variables in data mode
758 : !------------------------------------
759 :
760 0 : IF ( IsOldFile .AND. ntime > 1 ) THEN
761 0 : IF ( ThisDiagn%SpaceDim == 3 ) THEN
762 : CALL NC_READ_ARR( fID, TRIM(myName), 1, nlon, 1, nlat, &
763 0 : 1, nlev, 1, ntime-1, ncArr=Arr4DOld, RC=RC )
764 0 : Arr4D(:,:,:,1:ntime-1) = Arr4DOld(:,:,:,:)
765 : ELSE
766 : CALL NC_READ_ARR( fID, TRIM(myName), 1, nlon, 1, nlat, &
767 0 : -1, -1, 1, ntime-1, ncArr=Arr4DOld, RC=RC )
768 0 : Arr3D(:,:,1:ntime-1) = Arr4DOld(:,:,1,:)
769 : ENDIF
770 0 : IF ( ASSOCIATED(Arr4DOld) ) DEALLOCATE(Arr4DOld)
771 : ENDIF
772 :
773 : ! Mirror data and write to file. The mirroring is required in
774 : ! order to add the time dimension. Otherwise, the data would
775 : ! have no time information!
776 0 : IF ( ThisDiagn%SpaceDim == 3 ) THEN
777 0 : IF ( ASSOCIATED(ThisDiagn%Arr3D) ) THEN
778 0 : Arr4D(:,:,:,ntime) = ThisDiagn%Arr3D%Val
779 0 : Arr4D(:,:,:,1) = ThisDiagn%Arr3D%Val
780 : ENDIF
781 0 : CALL NC_VAR_WRITE ( fId, TRIM(myName), Arr4D=Arr4D )
782 : ELSE
783 0 : IF ( ASSOCIATED(ThisDiagn%Arr2D) ) THEN
784 0 : Arr3D(:,:,ntime) = ThisDiagn%Arr2D%Val
785 0 : Arr3D(:,:,1) = ThisDiagn%Arr2D%Val
786 : ENDIF
787 0 : CALL NC_VAR_WRITE ( fId, TRIM(myName), Arr3D=Arr3D )
788 : ENDIF
789 :
790 : ! verbose
791 0 : IF ( HCO_IsVerb(HcoState%Config%Err ) .AND. PS==1 ) THEN
792 0 : MSG = '--- Added diagnostics: '//TRIM(myName)
793 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
794 : ENDIF
795 : ENDIF
796 : ENDDO
797 : ENDDO
798 :
799 : !-----------------------------------------------------------------
800 : ! Cleanup
801 : !-----------------------------------------------------------------
802 :
803 : ! Close file
804 0 : CALL NC_CLOSE ( fId )
805 :
806 : ! Cleanup local variables
807 0 : Deallocate(Arr3D,Arr4D)
808 0 : ThisDiagn => NULL()
809 :
810 : ! Archive time stamp
811 : CALL DiagnCollection_Set ( HcoState%Diagn, COL=PS, &
812 0 : LastYMD=lymd, LastHMS=lhms, RC=RC )
813 :
814 : ! Return
815 0 : RC = HCO_SUCCESS
816 :
817 0 : END SUBROUTINE HCOIO_Write
818 : !EOC
819 : !------------------------------------------------------------------------------
820 : ! Harmonized Emissions Component (HEMCO) !
821 : !------------------------------------------------------------------------------
822 : !BOP
823 : !
824 : ! !IROUTINE: ConstructTimeStamp
825 : !
826 : ! !DESCRIPTION: Subroutine ConstructTimeStamp is a helper routine to construct
827 : ! the time stamp of a given diagnostics collection.
828 : !\\
829 : !\\
830 : ! !INTERFACE:
831 : !
832 0 : SUBROUTINE ConstructTimeStamp ( HcoState, PS, PrevTime, Yr, Mt, Dy, hr, mn, RC )
833 : !
834 : ! !USES:
835 : !
836 : USE HCO_State_Mod, ONLY : HCO_State
837 : USE HCO_Clock_Mod
838 : USE HCO_JULDAY_MOD
839 : !
840 : ! !INPUT/OUTPUT PARAMETERS:
841 : !
842 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state obj
843 : INTEGER, INTENT(IN ) :: PS ! collecion ID
844 : LOGICAL, INTENT(IN ) :: PrevTime ! Use previous time?
845 : !
846 : ! !INPUT/OUTPUT PARAMETERS:
847 : !
848 : INTEGER, INTENT(INOUT) :: RC ! Return code
849 : !
850 : ! !OUTPUT PARAMETERS:
851 : !
852 : INTEGER, INTENT( OUT) :: Yr
853 : INTEGER, INTENT( OUT) :: Mt
854 : INTEGER, INTENT( OUT) :: Dy
855 : INTEGER, INTENT( OUT) :: hr
856 : INTEGER, INTENT( OUT) :: mn
857 : !
858 : ! !REVISION HISTORY:
859 : ! 06 Nov 2015 - C. Keller - Initial version
860 : ! See https://github.com/geoschem/hemco for complete history
861 : !EOP
862 : !------------------------------------------------------------------------------
863 : !BOC
864 : !
865 : ! !LOCAL VARIABLES:
866 : !
867 : INTEGER :: Y2, M2, D2, h2, n2, s2
868 : INTEGER :: Y1, M1, D1, h1, n1, s1
869 : INTEGER :: LastYMD, LastHMS
870 : INTEGER :: YYYYMMDD, HHMMSS
871 : INTEGER :: OutTimeStamp
872 : REAL(dp) :: DAY, UTC, JD1, JD2, JDMID
873 : CHARACTER(LEN=255) :: MSG
874 : CHARACTER(LEN=255) :: LOC = 'ConstuctTimeStamp (hcoi_diagn_mod.F90)'
875 :
876 : !=================================================================
877 : ! ConstructTimeStamp begins here!
878 : !=================================================================
879 :
880 : ! Use HEMCO clock to create timestamp used in filename. Use previous
881 : ! time step if this option is selected.
882 0 : IF ( .NOT. PrevTime ) THEN
883 : CALL HcoClock_Get(HcoState%Clock,sYYYY=Y2,sMM=M2,&
884 0 : sDD=D2,sH=h2,sM=n2,sS=s2,RC=RC)
885 0 : IF ( RC /= HCO_SUCCESS ) THEN
886 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
887 0 : RETURN
888 : ENDIF
889 : ELSE
890 : CALL HcoClock_Get(HcoState%Clock,pYYYY=Y2,pMM=M2,&
891 0 : pDD=D2,pH=h2,pM=n2,pS=s2,RC=RC)
892 0 : IF ( RC /= HCO_SUCCESS ) THEN
893 0 : CALL HCO_ERROR( 'ERROR 9', RC, THISLOC=LOC )
894 0 : RETURN
895 : ENDIF
896 : ENDIF
897 :
898 : ! Get timestamp location for this collection
899 : CALL DiagnCollection_Get( HcoState%Diagn, PS, OutTimeStamp=OutTimeStamp, &
900 0 : LastYMD=LastYMD, LastHMS=LastHMS, RC=RC )
901 0 : IF ( RC /= HCO_SUCCESS ) THEN
902 0 : CALL HCO_ERROR( 'ERROR 10', RC, THISLOC=LOC )
903 0 : RETURN
904 : ENDIF
905 :
906 : ! Determine dates to be used:
907 :
908 : ! To use start date
909 0 : IF ( OutTimeStamp == HcoDiagnStart ) THEN
910 0 : Yr = FLOOR( MOD(LastYMD*1.d0, 100000000.d0 ) / 1.0d4 )
911 0 : Mt = FLOOR( MOD(LastYMD*1.d0, 10000.d0 ) / 1.0d2 )
912 0 : Dy = FLOOR( MOD(LastYMD*1.d0, 100.d0 ) / 1.0d0 )
913 0 : Hr = FLOOR( MOD(LastHMS*1.d0, 1000000.d0 ) / 1.0d4 )
914 0 : Mn = FLOOR( MOD(LastHMS*1.d0, 10000.d0 ) / 1.0d2 )
915 :
916 : ! Use mid point
917 0 : ELSEIF ( OutTimeStamp == HcoDiagnMid ) THEN
918 :
919 : ! Julian day of start interval:
920 0 : Y1 = FLOOR( MOD(LastYMD*1.d0, 100000000.d0 ) / 1.0d4 )
921 0 : M1 = FLOOR( MOD(LastYMD*1.d0, 10000.d0 ) / 1.0d2 )
922 0 : D1 = FLOOR( MOD(LastYMD*1.d0, 100.d0 ) / 1.0d0 )
923 0 : h1 = FLOOR( MOD(LastHMS*1.d0, 1000000.d0 ) / 1.0d4 )
924 0 : n1 = FLOOR( MOD(LastHMS*1.d0, 10000.d0 ) / 1.0d2 )
925 0 : s1 = FLOOR( MOD(LastHMS*1.d0, 100.d0 ) / 1.0d0 )
926 :
927 : UTC = ( REAL(h1,dp) / 24.0_dp ) + &
928 : ( REAL(n1,dp) / 1440.0_dp ) + &
929 0 : ( REAL(s1,dp) / 86400.0_dp )
930 0 : DAY = REAL(D1,dp) + UTC
931 0 : JD1 = JULDAY( Y1, M1, DAY )
932 :
933 : ! Julian day of end interval:
934 : UTC = ( REAL(h2,dp) / 24.0_dp ) + &
935 : ( REAL(n2,dp) / 1440.0_dp ) + &
936 0 : ( REAL(s2,dp) / 86400.0_dp )
937 0 : DAY = REAL(D2,dp) + UTC
938 0 : JD2 = JULDAY( Y2, M2, DAY )
939 :
940 : ! Julian day in the middle
941 0 : JDMID = ( JD1 + JD2 ) / 2.0_dp
942 :
943 : ! Tranlate back into dates
944 0 : CALL CALDATE( JDMID, YYYYMMDD, HHMMSS )
945 0 : Yr = FLOOR ( MOD( YYYYMMDD, 100000000) / 1.0e4_dp )
946 0 : Mt = FLOOR ( MOD( YYYYMMDD, 10000 ) / 1.0e2_dp )
947 0 : Dy = FLOOR ( MOD( YYYYMMDD, 100 ) / 1.0e0_dp )
948 0 : Hr = FLOOR ( MOD( HHMMSS, 1000000 ) / 1.0e4_dp )
949 0 : Mn = FLOOR ( MOD( HHMMSS, 10000 ) / 1.0e2_dp )
950 :
951 : ! Otherwise, use end date
952 : ELSE
953 0 : Yr = Y2
954 0 : Mt = M2
955 0 : Dy = D2
956 0 : Hr = h2
957 0 : Mn = n2
958 : ENDIF
959 :
960 : ! Return w/ success
961 0 : RC = HCO_SUCCESS
962 :
963 : END SUBROUTINE ConstructTimeStamp
964 : !EOC
965 : END MODULE HCOIO_WRITE_MOD
966 : #endif
|