1module gridStateVectorFileIO_mod
2 ! MODULE gridStateVectorFileIO_mod (prefix='gio' category='4. Data Object transformations')
3 !
4 !:Purpose: The grid-point state vector I/O methods for reading from and writing to
5 ! files.
6 !
7 use mpi
8 use midasMpi_mod
9 use gridStateVector_mod
10 use interpolation_mod
11 use utilities_mod
12 use verticalCoord_mod
13 use horizontalCoord_mod
14 use oceanMask_mod
15 use varNameList_mod
16 use ramDisk_mod
17 use timeCoord_mod
18 use mathPhysConstants_mod
19 use codePrecision_mod
20 use Vgrid_Descriptors
21 implicit none
22 save
23 private
24
25 ! public subroutines and functions
26 public :: gio_readFromFile, gio_readTrials, gio_readFile
27 public :: gio_readMaskFromFile
28 public :: gio_getMaskLAM
29 public :: gio_writeToFile
30 public :: gio_fileUnitsToStateUnits
31
32 integer, external :: get_max_rss
33
34 ! Namelist variables
35 logical :: interpToPhysicsGrid ! for LAM grid, choose to keep physics variables on their original grid
36
37 contains
38 !--------------------------------------------------------------------------
39 ! gio_readFromFile
40 !--------------------------------------------------------------------------
41 subroutine gio_readFromFile(statevector_out, fileName, etiket_in, typvar_in, &
42 stepIndex_opt, unitConversion_opt, &
43 statevectorRef_opt, readHeightSfc_opt, &
44 containsFullField_opt, vcoFileIn_opt)
45 !
46 !:Purpose: Read an RPN standard file and put the contents into a
47 ! stateVector object. Main high level wrapper subroutine.
48 !
49 implicit none
50
51 ! Arguments:
52 type(struct_gsv), intent(inout) :: statevector_out
53 character(len=*), intent(in) :: fileName
54 character(len=*), intent(in) :: etiket_in
55 character(len=*), intent(in) :: typvar_in
56 integer, optional, intent(in) :: stepIndex_opt
57 logical, optional, intent(in) :: unitConversion_opt
58 type(struct_gsv), optional, intent(in) :: statevectorRef_opt ! Reference statevector providing optional fields (P0, TT, HU)
59 logical, optional, intent(in) :: readHeightSfc_opt
60 logical, optional, intent(in) :: containsFullField_opt
61 type(struct_vco), optional, pointer, intent(in) :: vcoFileIn_opt
62
63 ! Locals:
64 logical :: doHorizInterp, doVertInterp, unitConversion
65 logical :: readHeightSfc, containsFullField
66 logical :: foundVarNameInFile
67 integer :: stepIndex, varIndex
68 character(len=4) :: varName
69 type(struct_vco), pointer :: vco_file
70 type(struct_hco), pointer :: hco_file
71
72 nullify(vco_file, hco_file)
73
74 write(*,*)
75 write(*,*) 'gio_readFromFile: START'
76 call utl_tmg_start(160,'low-level--gsv_readFromFile')
77
78 if ( present(stepIndex_opt) ) then
79 stepIndex = stepIndex_opt
80 else
81 stepIndex = statevector_out%anltime
82 end if
83 if (stepIndex > stateVector_out%numStep .or. stepIndex < 1) then
84 write(*,*) 'stepIndex = ', stepIndex
85 call utl_abort('gio_readFromFile: invalid value for stepIndex')
86 end if
87
88 if ( present(unitConversion_opt) ) then
89 unitConversion = unitConversion_opt
90 else
91 unitConversion = .true.
92 end if
93
94 if (present(containsFullField_opt)) then
95 containsFullField = containsFullField_opt
96 else
97 containsFullField = .true.
98 end if
99 write(*,*) 'gio_readFromFile: containsFullField = ', containsFullField
100
101 if ( present(readHeightSfc_opt) ) then
102 readHeightSfc = readHeightSfc_opt
103 else
104 readHeightSfc = .false.
105 end if
106
107 ! set up vertical and horizontal coordinate for input file
108 if (present(vcoFileIn_opt)) then
109 vco_file => vcoFileIn_opt
110 write(*,*)
111 write(*,*) 'gio_readFromFile: vertical levels defined in user-supplied vco object will be read'
112 else
113 call vco_setupFromFile(vco_file,trim(fileName), beSilent_opt=.true.)
114 write(*,*)
115 write(*,*) 'gio_readFromFile: all the vertical levels will be read from ', trim(fileName)
116 end if
117
118 foundVarNameInFile = .false.
119
120 do varIndex = 1, vnl_numvarmax
121 varName = vnl_varNameList(varIndex)
122
123 if (.not. gsv_varExist(statevector_out,varName)) cycle
124
125 ! make sure variable is in the file
126 if (.not. utl_varNamePresentInFile(varName,fileName_opt=trim(fileName))) cycle
127
128 ! adopt a variable on the full/dynamic LAM grid
129 if (.not. statevector_out%hco%global .and. (trim(varName) == 'TM' .or. trim(varName) == 'MG')) cycle
130
131 foundVarNameInFile = .true.
132
133 exit
134
135 end do
136
137 ! special case when only TM (Surface Temperature) is in the file:
138 if (.not. foundVarNameInFile) then
139 varname = 'TM'
140 if (gsv_varExist( statevector_out, varname) .and. &
141 utl_varNamePresentInFile(varname, fileName_opt = trim(fileName))) &
142 foundVarNameInFile = .true.
143 end if
144
145 ! to be safe for situations where, e.g. someone wants to only read MG from a file
146 if (.not. foundVarNameInFile) then
147 varname = 'P0'
148 if (utl_varNamePresentInFile( varname, fileName_opt = trim( fileName))) &
149 foundVarNameInFile = .true.
150 end if
151
152 if (.not. foundVarNameInFile) call utl_abort('gio_readFromFile: NO variables found in the file!!!')
153
154 write(*,*) 'gio_readFromFile: defining hco by varname= ', varName
155
156 call hco_setupFromFile(hco_file, trim(fileName), etiket_in, gridName_opt='FILEGRID', varName_opt = varName)
157
158 ! test if horizontal and/or vertical interpolation needed for statevector grid
159 doVertInterp = .not.vco_equal(vco_file,statevector_out%vco)
160 doHorizInterp = .not.hco_equal(hco_file,statevector_out%hco)
161 write(*,*) 'gio_readFromFile: doVertInterp = ', doVertInterp, ', doHorizInterp = ', doHorizInterp
162
163 ! call appropriate subroutine to do actual work
164 if ((doVertInterp .or. doHorizInterp) .and. statevector_out%mpi_distribution=='Tiles') then
165 call readFromFileAndInterpToTiles(statevector_out, fileName, &
166 vco_file, hco_file, etiket_in, typvar_in, stepIndex, unitConversion, &
167 readHeightSfc, containsFullField, statevectorRef_opt=statevectorRef_opt)
168 else if ((doVertInterp .or. doHorizInterp) .and. .not.stateVector_out%mpi_local) then
169 call readFromFileAndInterp1Proc(statevector_out, fileName, &
170 vco_file, hco_file, etiket_in, typvar_in, stepIndex, unitConversion, &
171 readHeightSfc, containsFullField)
172 else if (.not.(doVertInterp .or. doHorizInterp) .and. stateVector_out%mpi_local) then
173 call readFromFileAndTransposeToTiles(statevector_out, fileName, &
174 etiket_in, typvar_in, stepIndex, unitConversion, &
175 readHeightSfc, containsFullField)
176 else
177 call readFromFileOnly(statevector_out, fileName, &
178 etiket_in, typvar_in, stepIndex, unitConversion, &
179 readHeightSfc, containsFullField)
180 end if
181
182 call utl_tmg_stop(160)
183 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
184 write(*,*) 'gio_readFromFile: END'
185
186 end subroutine gio_readFromFile
187
188 !--------------------------------------------------------------------------
189 ! readFromFileAndInterpToTiles
190 !--------------------------------------------------------------------------
191 subroutine readFromFileAndInterpToTiles(statevector_out, fileName, &
192 vco_file, hco_file, etiket_in, typvar_in, stepIndex, unitConversion, &
193 readHeightSfc, containsFullField, statevectorRef_opt)
194 !
195 !:Purpose: Read an RPN standard file and put the contents into a
196 ! stateVector object. Wrapper subroutine that also proceed with
197 ! distributed interpolation on MPI tiles.
198 !
199 !:Note: this routine currently only works correctly for reading FULL FIELDS,
200 ! not increments or perturbations... because of the HU -> LQ conversion
201 implicit none
202
203 ! Arguments:
204 type(struct_gsv), intent(inout) :: statevector_out
205 character(len=*), intent(in) :: fileName
206 type(struct_vco), pointer, intent(in) :: vco_file
207 type(struct_hco), pointer, intent(in) :: hco_file
208 character(len=*), intent(in) :: etiket_in
209 character(len=*), intent(in) :: typvar_in
210 integer, intent(in) :: stepIndex
211 logical, intent(in) :: unitConversion
212 logical, intent(in) :: readHeightSfc
213 logical, intent(in) :: containsFullField
214 type(struct_gsv), optional, intent(in) :: statevectorRef_opt ! Reference statevector providing optional fields (P0, TT, HU)
215
216 ! Locals:
217 real(4), pointer :: field3d_r4_ptr(:,:,:)
218 real(8), pointer :: field_in_ptr(:,:,:,:), field_out_ptr(:,:,:,:)
219 character(len=4), pointer :: varNamesToRead(:)
220 type(struct_gsv) :: statevector_file_r4, statevector_tiles
221 type(struct_gsv) :: statevector_hinterp_r4, statevector_vinterp
222
223 nullify(field3d_r4_ptr, field_in_ptr, field_out_ptr)
224
225 write(*,*) ''
226 write(*,*) 'readFromFileAndInterpToTiles: START'
227
228 nullify(varNamesToRead)
229 call gsv_varNamesList(varNamesToRead, statevector_out)
230
231 !-- 1.0 Read the file, distributed over mpi task with respect to variables/levels
232
233 ! initialize single precision 3D working copy of statevector for reading file
234 call gsv_allocate(statevector_file_r4, 1, hco_file, vco_file, &
235 dateStamp_opt=statevector_out%datestamplist(stepIndex), &
236 mpi_local_opt=.true., mpi_distribution_opt='VarsLevs', &
237 dataKind_opt=4, allocHeightSfc_opt=readHeightSfc, &
238 varNames_opt=varNamesToRead, &
239 hInterpolateDegree_opt=statevector_out%hInterpolateDegree, &
240 hExtrapolateDegree_opt=statevector_out%hExtrapolateDegree )
241
242 call gio_readFile(statevector_file_r4, filename, etiket_in, typvar_in, &
243 containsFullField, readHeightSfc_opt=readHeightSfc)
244
245 !-- 2.0 Horizontal Interpolation
246
247 ! initialize single precision 3D working copy of statevector for horizontal interpolation result
248 call gsv_allocate(statevector_hinterp_r4, 1, statevector_out%hco, vco_file, &
249 dateStamp_opt=statevector_out%datestamplist(stepIndex), &
250 mpi_local_opt=.true., mpi_distribution_opt='VarsLevs', &
251 dataKind_opt=4, allocHeightSfc_opt=readHeightSfc, &
252 varNames_opt=varNamesToRead, &
253 hInterpolateDegree_opt=statevector_out%hInterpolateDegree, &
254 hExtrapolateDegree_opt=statevector_out%hExtrapolateDegree )
255
256 call int_hInterp_gsv(statevector_file_r4, statevector_hinterp_r4)
257
258 call gsv_deallocate(statevector_file_r4)
259
260 !-- 3.0 Unit conversion
261
262 if ( unitConversion ) then
263 call gio_fileUnitsToStateUnits( statevector_hinterp_r4, containsFullField )
264 end if
265
266 !-- 4.0 MPI communication from vars/levels to lat/lon tiles
267
268 ! initialize double precision 3D working copy of statevector for mpi communication result
269 call gsv_allocate(statevector_tiles, 1, statevector_out%hco, vco_file, &
270 dateStamp_opt=statevector_out%datestamplist(stepIndex), &
271 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
272 dataKind_opt=8, allocHeightSfc_opt=readHeightSfc, &
273 varNames_opt=varNamesToRead )
274
275 call gsv_transposeVarsLevsToTiles(statevector_hinterp_r4, statevector_tiles)
276
277 call gsv_deallocate(statevector_hinterp_r4)
278
279 !-- 5.0 Vertical interpolation
280
281 ! initialize double precision 3D working copy of statevector for mpi communication result
282 call gsv_allocate(statevector_vinterp, 1, statevector_out%hco, statevector_out%vco, &
283 dateStamp_opt=statevector_out%datestamplist(stepIndex), &
284 mpi_local_opt=.true., mpi_distribution_opt='Tiles', dataKind_opt=8, &
285 allocHeightSfc_opt=readHeightSfc, varNames_opt=varNamesToRead )
286
287 call int_vInterp_gsv( statevector_tiles, statevector_vinterp, &
288 statevectorRef_opt=statevectorRef_opt)
289
290 call gsv_deallocate(statevector_tiles)
291
292 !-- 6.0 Copy result to output statevector
293 call gsv_copy(statevector_vinterp, statevector_out, stepIndexOut_opt=stepIndex)
294
295 call gsv_deallocate(statevector_vinterp)
296 deallocate(varNamesToRead)
297
298 write(*,*) 'readFromFileAndInterpToTiles: END'
299
300 end subroutine readFromFileAndInterpToTiles
301
302 !--------------------------------------------------------------------------
303 ! readFromFileAndTransposeToTiles
304 !--------------------------------------------------------------------------
305 subroutine readFromFileAndTransposeToTiles(statevector_out, fileName, &
306 etiket_in, typvar_in, stepIndex, unitConversion, &
307 readHeightSfc, containsFullField)
308 !
309 !:Purpose: Read an RPN standard file and put the contents into a
310 ! stateVector object. Wrapper subroutine that also proceed with
311 ! distributed transposition on MPI tiles.
312 !
313 !:Note: this routine currently only works correctly for reading FULL FIELDS,
314 ! not increments or perturbations... because of the HU -> LQ conversion
315 implicit none
316
317 ! Arguments:
318 type(struct_gsv), intent(inout) :: statevector_out
319 character(len=*), intent(in) :: fileName
320 character(len=*), intent(in) :: etiket_in
321 character(len=*), intent(in) :: typvar_in
322 integer, intent(in) :: stepIndex
323 logical, intent(in) :: unitConversion
324 logical, intent(in) :: readHeightSfc
325 logical, intent(in) :: containsFullField
326
327 ! Locals:
328 real(4), pointer :: field3d_r4_ptr(:,:,:)
329 real(8), pointer :: field_in_ptr(:,:,:,:), field_out_ptr(:,:,:,:)
330 character(len=4), pointer :: varNamesToRead(:)
331 type(struct_gsv) :: statevector_file_r4, statevector_tiles
332
333 nullify(field3d_r4_ptr, field_in_ptr, field_out_ptr)
334
335 write(*,*) ''
336 write(*,*) 'readFromFileAndTransposeToTiles: START'
337
338 nullify(varNamesToRead)
339 call gsv_varNamesList(varNamesToRead, statevector_out)
340
341 !-- 1.0 Read the file, distributed over mpi task with respect to variables/levels
342
343 ! initialize single precision 3D working copy of statevector for reading file
344 call gsv_allocate(statevector_file_r4, 1, statevector_out%hco, statevector_out%vco, &
345 dateStamp_opt=statevector_out%datestamplist(stepIndex), &
346 mpi_local_opt=.true., mpi_distribution_opt='VarsLevs', &
347 dataKind_opt=4, allocHeightSfc_opt=readHeightSfc, &
348 varNames_opt=varNamesToRead )
349
350 call gio_readFile(statevector_file_r4, filename, etiket_in, typvar_in, &
351 containsFullField, readHeightSfc_opt=readHeightSfc)
352
353 !-- 2.0 Unit conversion
354 if ( unitConversion ) then
355 call gio_fileUnitsToStateUnits( statevector_file_r4, containsFullField )
356 end if
357
358 !-- 3.0 MPI communication from vars/levels to lat/lon tiles
359 call gsv_allocate(statevector_tiles, 1, statevector_out%hco, statevector_out%vco, &
360 dateStamp_opt=statevector_out%datestamplist(stepIndex), &
361 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
362 dataKind_opt=8, allocHeightSfc_opt=readHeightSfc, &
363 varNames_opt=varNamesToRead)
364
365 call gsv_transposeVarsLevsToTiles(statevector_file_r4, statevector_tiles)
366
367 !-- 4.0 Copy result to output statevector
368 call gsv_copy(statevector_tiles, statevector_out, stepIndexOut_opt=stepIndex)
369
370 call gsv_deallocate(statevector_tiles)
371 call gsv_deallocate(statevector_file_r4)
372 deallocate(varNamesToRead)
373
374 write(*,*) 'readFromFileAndTransposeToTiles: END'
375
376 end subroutine readFromFileAndTransposeToTiles
377
378 !--------------------------------------------------------------------------
379 ! readFromFileAndInterp1Proc
380 !--------------------------------------------------------------------------
381 subroutine readFromFileAndInterp1Proc(statevector_out_r4, fileName, &
382 vco_file, hco_file, etiket_in, typvar_in, stepIndex, unitConversion, &
383 readHeightSfc, containsFullField)
384 !
385 !:Purpose: Read an RPN standard file and put the contents into a
386 ! stateVector object. Wrapper subroutine that also proceed with
387 ! (serial) interpolation.
388 !
389 implicit none
390
391 ! Arguments:
392 type(struct_gsv), intent(inout) :: statevector_out_r4
393 character(len=*), intent(in) :: fileName
394 type(struct_vco), pointer, intent(in) :: vco_file
395 type(struct_hco), pointer, intent(in) :: hco_file
396 character(len=*), intent(in) :: etiket_in
397 character(len=*), intent(in) :: typvar_in
398 integer, intent(in) :: stepIndex
399 logical, intent(in) :: unitConversion
400 logical, intent(in) :: readHeightSfc
401 logical, intent(in) :: containsFullField
402
403 ! Locals:
404 real(4), pointer :: field_in_ptr(:,:,:,:), field_out_ptr(:,:,:,:)
405 character(len=4), pointer :: varNamesToRead(:)
406 type(struct_gsv) :: statevector_file_r4, statevector_hinterp_r4
407 type(struct_gsv) :: statevector_vinterp_r4
408
409 nullify(field_in_ptr, field_out_ptr)
410
411 write(*,*) ''
412 write(*,*) 'readFromFileAndInterp1Proc: START'
413
414 nullify(varNamesToRead)
415 call gsv_varNamesList(varNamesToRead, statevector_out_r4)
416
417 !-- 1.0 Read the file
418
419 ! initialize single precision 3D working copy of statevector for reading file
420 call gsv_allocate(statevector_file_r4, 1, hco_file, vco_file, &
421 dateStamp_opt=statevector_out_r4%datestamplist(stepIndex), &
422 mpi_local_opt=.false., dataKind_opt=4, &
423 allocHeightSfc_opt=readHeightSfc, varNames_opt=varNamesToRead, &
424 hInterpolateDegree_opt=statevector_out_r4%hInterpolateDegree, &
425 hExtrapolateDegree_opt=statevector_out_r4%hExtrapolateDegree)
426
427 call gio_readFile(statevector_file_r4, filename, etiket_in, typvar_in, &
428 containsFullField, readHeightSfc_opt=readHeightSfc)
429
430 !-- 2.0 Horizontal Interpolation
431
432 ! initialize single precision 3D working copy of statevector for horizontal interpolation result
433 call gsv_allocate(statevector_hinterp_r4, 1, statevector_out_r4%hco, vco_file, &
434 dateStamp_opt=statevector_out_r4%datestamplist(stepIndex), &
435 mpi_local_opt=.false., dataKind_opt=4, &
436 allocHeightSfc_opt=readHeightSfc, varNames_opt=varNamesToRead, &
437 hInterpolateDegree_opt=statevector_out_r4%hInterpolateDegree, &
438 hExtrapolateDegree_opt=statevector_out_r4%hExtrapolateDegree)
439
440 call int_hInterp_gsv(statevector_file_r4, statevector_hinterp_r4)
441
442 call gsv_deallocate(statevector_file_r4)
443
444 !-- 3.0 Unit conversion (must come before vertical interp to get Psfc in Pascals)
445
446 if ( unitConversion ) then
447 call gio_fileUnitsToStateUnits( statevector_hinterp_r4, containsFullField )
448 end if
449
450 !-- 4.0 Vertical interpolation
451
452 ! initialize double precision 3D working copy of statevector for mpi communication result
453 call gsv_allocate(statevector_vinterp_r4, 1, statevector_out_r4%hco, statevector_out_r4%vco, &
454 dateStamp_opt=statevector_out_r4%datestamplist(stepIndex), &
455 mpi_local_opt=.false., dataKind_opt=4, &
456 allocHeightSfc_opt=readHeightSfc, varNames_opt=varNamesToRead)
457
458 call int_vInterp_gsv(statevector_hinterp_r4,statevector_vinterp_r4)
459
460 call gsv_deallocate(statevector_hinterp_r4)
461
462 !-- 5.0 Copy result to output statevector
463 call gsv_copy(statevector_vinterp_r4, statevector_out_r4, stepIndexOut_opt=stepIndex)
464
465 call gsv_deallocate(statevector_vinterp_r4)
466 deallocate(varNamesToRead)
467
468 write(*,*) 'readFromFileAndInterp1Proc: END'
469
470 end subroutine readFromFileAndInterp1Proc
471
472 !--------------------------------------------------------------------------
473 ! readFromFileOnly
474 !--------------------------------------------------------------------------
475 subroutine readFromFileOnly(statevector_out, fileName, etiket_in, typvar_in, &
476 stepIndex, unitConversion, readHeightSfc, &
477 containsFullField)
478 !
479 !:Purpose: Read an RPN standard file and put the contents into a
480 ! stateVector object. Wrapper subroutine
481 !
482 implicit none
483
484 ! Arguments:
485 type(struct_gsv), intent(inout) :: statevector_out
486 character(len=*), intent(in) :: fileName
487 character(len=*), intent(in) :: etiket_in
488 character(len=*), intent(in) :: typvar_in
489 integer , intent(in) :: stepIndex
490 logical , intent(in) :: unitConversion
491 logical , intent(in) :: readHeightSfc
492 logical , intent(in) :: containsFullField
493
494 write(*,*) ''
495 write(*,*) 'readFromFileOnly: Do simple reading with no interpolation and no mpi redistribution'
496
497 call gio_readFile(statevector_out, filename, etiket_in, typvar_in, &
498 containsFullField, readHeightSfc_opt = readHeightSfc, &
499 stepIndex_opt = stepIndex)
500
501 if (unitConversion) then
502 call gio_fileUnitsToStateUnits(statevector_out, containsFullField, stepIndex_opt = stepIndex)
503 end if
504
505 write(*,*) 'readFromFileOnly: END'
506
507 end subroutine readFromFileOnly
508
509 !--------------------------------------------------------------------------
510 ! gio_readFile
511 !--------------------------------------------------------------------------
512 subroutine gio_readFile(statevector, filename, etiket_in, typvar_in, &
513 containsFullField, readHeightSfc_opt, stepIndex_opt, &
514 ignoreDate_opt)
515 !
516 !:Purpose: Read an RPN standard file and put the contents into a
517 ! stateVector object. Low level subroutine that does the actual
518 ! file reading.
519 !
520 implicit none
521
522 ! Arguments:
523 type(struct_gsv), intent(inout) :: statevector
524 character(len=*), intent(in) :: fileName
525 character(len=*), intent(in) :: etiket_in
526 character(len=*), intent(in) :: typvar_in
527 logical, intent(in) :: containsFullField
528 logical, optional, intent(in) :: readHeightSfc_opt
529 integer, optional, intent(in) :: stepIndex_opt
530 logical, optional, intent(in) :: ignoreDate_opt
531
532 ! Locals:
533 integer :: nulfile, ierr, ip1, ni_file, nj_file, nk_file, kIndex, stepIndex
534 integer :: ikey, levIndex
535 integer :: stepIndexBeg, stepIndexEnd, ni_var, nj_var, nk_var
536 integer :: fnom, fstouv, fclos, fstfrm, fstlir, fstinf
537 integer :: fstprm, EZscintID_var, ezdefset, ezqkdef
538 integer :: dateo_var, deet_var, npas_var, nbits_var, datyp_var
539 integer :: ip1_var, ip2_var, ip3_var, swa_var, lng_var, dltf_var, ubc_var
540 integer :: extra1_var, extra2_var, extra3_var
541 integer :: ig1_var, ig2_var, ig3_var, ig4_var
542 integer :: varIndex, dateStampList(statevector%numStep)
543 character(len=4 ) :: nomvar_var
544 character(len=2 ) :: typvar_var
545 character(len=1 ) :: grtyp_var
546 character(len=12) :: etiket_var
547 real(4), pointer :: field_r4_ptr(:,:,:,:)
548 real(8), pointer :: field_r8_ptr(:,:,:,:)
549 real(4), pointer :: gd2d_file_r4(:,:), gd2d_r4_UV_ptr(:,:,:)
550 real(8), pointer :: gd2d_r8_UV_ptr(:,:,:)
551 real(8), pointer :: heightSfc_ptr(:,:)
552 real(4), allocatable :: gd2d_var_r4(:,:)
553 character(len=4) :: varName, varNameToRead
554 character(len=4) :: varLevel
555 type(struct_vco), pointer :: vco_file
556 type(struct_hco), pointer :: hco_file
557 logical :: foundVarNameInFile, ignoreDate
558
559 write(*,*) 'gio_readFile: starting'
560 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
561
562 call readNml()
563
564 vco_file => gsv_getVco(statevector)
565
566 if (statevector%mpi_distribution /= 'VarsLevs' .and. &
567 statevector%mpi_local) then
568 call utl_abort('gio_readFile: statevector must have ' // &
569 'complete horizontal fields on each mpi task.')
570 end if
571
572 if (present(stepIndex_opt)) then
573 stepIndexBeg = stepIndex_opt
574 stepIndexEnd = stepIndex_opt
575 else
576 stepIndexBeg = 1
577 stepIndexEnd = statevector%numStep
578 end if
579
580 if (present(ignoreDate_opt)) then
581 ignoreDate = ignoreDate_opt
582 else
583 ignoreDate = .false.
584 end if
585
586 if (.not. associated(statevector%dateStampList)) then
587 call utl_abort('gio_readFile: dateStampList of statevector is not associated with a target!')
588 else
589 dateStampList(:) = statevector%dateStampList(:)
590 if (ignoreDate) then
591 write(*,*) 'gio_readFile: as requested, ignoring the date when reading fields'
592 dateStampList(:) = -1
593 end if
594 end if
595
596 !- Open input field
597 nulfile = 0
598 write(*,*) 'gio_readFile: file name = ', trim(fileName)
599 ierr = fnom(nulfile, trim(fileName), 'RND+OLD+R/O', 0)
600
601 if (ierr >= 0) then
602 ierr = fstouv(nulfile,'RND+OLD')
603 else
604 call utl_abort('gio_readFile: problem opening input file')
605 end if
606
607 if (nulfile == 0) then
608 call utl_abort('gio_readFile: unit number for input file not valid')
609 end if
610
611 ! Read surface height if requested
612 if (present(readHeightSfc_opt)) then
613 if (readHeightSfc_opt .and. gsv_isAssocHeightSfc(statevector)) then
614 write(*,*) 'gio_readFile: reading the surface height'
615 varName = 'GZ'
616 ip1 = statevector%vco%ip1_sfc
617 typvar_var = typvar_in
618 ikey = fstinf(nulfile, ni_file, nj_file, nk_file, &
619 -1, etiket_in, &
620 -1, -1, -1, typvar_var, varName)
621
622 if (ikey < 0) then
623 if (trim(typvar_in) /= "") then
624 typvar_var(2:2) = '@'
625 ikey = fstinf(nulfile, ni_file, nj_file, nk_file, &
626 -1, etiket_in, &
627 -1, -1, -1, typvar_var, varName)
628 end if
629 if (ikey < 0) then
630 write(*,*) 'gio_readFile: etiket_in = ', etiket_in
631 write(*,*) 'gio_readFile: typvar_in = ', typvar_in
632 call utl_abort('gio_readFile: Problem with reading surface height from file')
633 end if
634 end if
635
636 if (ni_file /= statevector%hco%ni .or. nj_file /= statevector%hco%nj) then
637 write(*,*) 'ni, nj in file = ', ni_file, nj_file
638 write(*,*) 'ni, nj in statevector = ', statevector%hco%ni, statevector%hco%nj
639 call utl_abort('gio_readFile: Dimensions of surface height not consistent')
640 end if
641
642 allocate(gd2d_file_r4(ni_file, nj_file))
643 gd2d_file_r4(:,:) = 0.0d0
644 ierr = fstlir(gd2d_file_r4(:,:), nulfile, ni_file, nj_file, nk_file, &
645 -1, etiket_in, ip1, -1, -1, typvar_var,varName)
646 if (ierr < 0) then
647 write(*,*) 'ip1 = ', ip1
648 write(*,*) 'etiket_in = ', etiket_in
649 write(*,*) 'typvar_var = ', typvar_var
650 call utl_abort('gio_readFile: Problem with reading surface height from file')
651 end if
652 heightSfc_ptr => gsv_getHeightSfc(statevector)
653 heightSfc_ptr = real(gd2d_file_r4(1:gsv_getHco(statevector)%ni, &
654 1:gsv_getHco(statevector)%nj), 8) * 10.0d0
655 deallocate(gd2d_file_r4)
656 end if
657 end if
658
659 nullify(hco_file)
660 nullify(gd2d_file_r4)
661 if (statevector%mykCount > 0) then
662 if (statevector%hco%global) then
663
664 foundVarNameInFile = .false.
665 do varIndex = 1, vnl_numvarmax
666 varName = vnl_varNameList(varIndex)
667
668 if (.not. gsv_varExist(statevector, varName)) cycle
669
670 ! make sure variable is in the file
671 if (.not. utl_varNamePresentInFile(varName, fileName_opt = trim(fileName))) cycle
672
673 ! adopt a variable on the full/dynamic LAM grid
674 if ((trim(varName) == 'TM' .or. trim(varName) == 'MG')) cycle
675
676 foundVarNameInFile = .true.
677
678 exit
679 end do
680
681 ! special case when only TM (Surface Temperature) is in the file:
682 if (.not. foundVarNameInFile) then
683 varname = 'TM'
684 if (gsv_varExist(statevector, varname) .and. &
685 utl_varNamePresentInFile(varname, fileName_opt = trim(fileName))) &
686 foundVarNameInFile = .true.
687 end if
688
689 ! to be safe for situations where, e.g. someone wants to only read MG from a file
690 if (.not. foundVarNameInFile) then
691 varname = 'P0'
692 if (utl_varNamePresentInFile(varname, fileName_opt = trim( fileName))) &
693 foundVarNameInFile = .true.
694 end if
695
696 if (.not. foundVarNameInFile) call utl_abort('gio_readFile: NO variable is in the file')
697
698 call hco_setupFromFile(hco_file, filename, ' ', 'INPUTFILE', varName_opt = varName)
699
700 else
701 ! In LAM mode, force the input file dimensions to be always identical to the input statevector dimensions
702 hco_file => statevector%hco
703
704 ! Also attempt to set up the physics grid
705 if (interpToPhysicsGrid) then
706 var_loop: do varIndex = 1, vnl_numvarmax
707 varName = vnl_varNameList(varIndex)
708 if (.not. gsv_varExist(statevector, varName)) cycle var_loop
709 if (.not. vnl_isPhysicsVar(varName)) cycle var_loop
710 if (utl_varNamePresentInFile(varName, fileName_opt = filename) .and. &
711 .not. associated(statevector%hco_physics)) then
712 write(*,*) 'gio_readFile: set up physics grid using the variable:', varName
713 call hco_SetupFromFile(statevector%hco_physics, filename, ' ', &
714 'INPUTFILE', varName_opt = varName)
715 exit var_loop
716 end if
717 end do var_loop
718 end if
719
720 end if
721 allocate(gd2d_file_r4(hco_file%ni, hco_file%nj))
722 gd2d_file_r4(:,:) = 0.0
723 end if
724
725 ! Read all other fields needed for this MPI task
726 if (statevector%dataKind == 4) then
727 call gsv_getField(statevector, field_r4_ptr)
728 else
729 call gsv_getField(statevector, field_r8_ptr)
730 end if
731
732 do stepIndex = stepIndexBeg, stepIndexEnd
733 k_loop: do kIndex = statevector%mykBeg, statevector%mykEnd
734 varName = gsv_getVarNameFromK(statevector, kIndex)
735 levIndex = gsv_getLevFromK(statevector, kIndex)
736
737 if (.not.gsv_varExist(statevector, varName)) cycle k_loop
738
739 ! Check that the wanted field is present in the file
740 if (utl_varNamePresentInFile(varName, fileUnit_opt = nulfile)) then
741 varNameToRead = varName
742 else
743 select case (trim(varName))
744 case ('LVIS')
745 varNameToRead = 'VIS'
746 case ('Z_T', 'Z_M', 'P_T', 'P_M')
747 cycle k_loop
748 case ('LPR')
749 varNameToRead = 'PR'
750 case default
751 call utl_abort('gio_readFile: variable '//trim(varName)//' was not found in '//trim(fileName))
752 end select
753 end if
754
755 varLevel = vnl_varLevelFromVarname(varNameToRead)
756 if (varLevel == 'MM') then
757 ip1 = vco_file%ip1_M(levIndex)
758 else if (varLevel == 'TH') then
759 ip1 = vco_file%ip1_T(levIndex)
760 else if (varLevel == 'SF') then
761 ip1 = -1
762 else if (varLevel == 'SFTH') then
763 ip1 = vco_file%ip1_T_2m
764 else if (varLevel == 'SFMM') then
765 ip1 = vco_file%ip1_M_10m
766 else if (varLevel == 'OT') then
767 ip1 = vco_ip1_other(levIndex)
768 else if (varLevel == 'DP') then
769 ip1 = vco_file%ip1_depth(levIndex)
770 else if (varLevel == 'SS') then
771 ip1 = -1
772 else
773 write(*,*) 'varLevel =', varLevel
774 call utl_abort('gio_readFile: unknown varLevel')
775 end if
776
777 typvar_var = typvar_in
778
779 ! Make sure that the input variable has the same grid size than hco_file
780 ikey = fstinf(nulfile, ni_var, nj_var, nk_var, &
781 datestamplist(stepIndex), etiket_in, &
782 -1, -1, -1, typvar_var, varNameToRead)
783
784 if (ikey < 0) then
785 if (trim(typvar_in) /= '') then
786 typvar_var(2:2) = '@'
787 ikey = fstinf(nulfile, ni_var, nj_var, nk_var, &
788 datestamplist(stepIndex), etiket_in, &
789 -1, -1, -1, typvar_var, varNameToRead)
790 end if
791 if (ikey < 0) then
792 write(*,*) 'gio_readFile: looking for datestamp = ', datestamplist(stepIndex)
793 write(*,*) 'gio_readFile: etiket_in = ', etiket_in
794 write(*,*) 'gio_readFile: typvar_in = ', typvar_in
795 call utl_abort('gio_readFile: cannot find field ' // trim(varNameToRead) // ' in file ' // trim(fileName))
796 end if
797 end if
798
799 ierr = fstprm(ikey, & ! IN
800 dateo_var, deet_var, npas_var, ni_var, nj_var, nk_var, nbits_var, & ! OUT
801 datyp_var, ip1_var, ip2_var, ip3_var, typvar_var, nomvar_var, & ! OUT
802 etiket_var, grtyp_var, ig1_var, ig2_var, ig3_var, ig4_var, swa_var, & ! OUT
803 lng_var, dltf_var, ubc_var, extra1_var, extra2_var, extra3_var ) ! OUT
804 statevector%deet = deet_var
805 statevector%ip2List(stepIndex) = ip2_var
806 statevector%npasList(stepIndex) = npas_var
807 statevector%dateOriginList(stepIndex) = dateo_var
808 statevector%etiket = etiket_var
809
810 ! Check if we found a mask field by mistake - if yes, need to fix the code!
811 if (typvar_var == '@@') then
812 call utl_abort('gio_readFile: read a mask file by mistake - need to modify file or fix the code')
813 end if
814
815 if (ni_var == hco_file%ni .and. nj_var == hco_file%nj) then
816 ierr = fstlir(gd2d_file_r4(:,:), nulfile, ni_file, nj_file, nk_file, &
817 datestamplist(stepIndex), etiket_in, ip1, -1, -1, &
818 typvar_var, varNameToRead)
819 else
820 ! Special cases for variables that are on a different horizontal grid in LAM (e.g. TG)
821 write(*,*)
822 write(*,*) 'gio_readFile: variable on a different horizontal grid = ',trim(varNameToRead)
823 write(*,*) ni_var, hco_file%ni, nj_var, hco_file%nj
824 if (interpToPhysicsGrid) then
825 if (associated(statevector%hco_physics)) then
826 if (ni_var == statevector%hco_physics%ni .and. &
827 nj_var == statevector%hco_physics%nj) then
828 write(*,*) 'gio_readFile: this variable on same grid as other physics variables'
829 statevector%onPhysicsGrid(vnl_varListIndex(varName)) = .true.
830 else
831 call utl_abort('gio_readFile: this variable not on same grid as other physics variables')
832 end if
833 else
834 call utl_abort('gio_readFile: physics grid has not been set up')
835 end if
836 end if
837
838 if (statevector%hco%global) then
839 call utl_abort('gio_readFile: This is not allowed in global mode!')
840 end if
841
842 EZscintID_var = ezqkdef(ni_var, nj_var, grtyp_var, ig1_var, ig2_var, ig3_var, ig4_var, nulfile) ! IN
843
844 allocate(gd2d_var_r4(ni_var, nj_var))
845 gd2d_var_r4(:,:) = 0.0
846
847 ierr = fstlir(gd2d_var_r4(:,:), nulfile, ni_var, nj_var, nk_var, &
848 datestamplist(stepIndex), etiket_in, ip1, -1, -1, &
849 typvar_in, varNameToRead)
850
851 ierr = ezdefset(hco_file%EZscintID, EZscintID_var)
852 ierr = int_hInterpScalar(gd2d_file_r4, gd2d_var_r4, &
853 interpDegree = 'NEAREST', extrapDegree_opt = 'NEUTRAL')
854
855 ! read the corresponding mask if it exists
856 if (typvar_var(2:2) == '@') then
857 write(*,*) 'gio_readFile: read mask that needs interpolation for variable name: ', nomvar_var
858 call utl_abort('gio_readFile: not implemented yet')
859 end if
860
861 deallocate(gd2d_var_r4)
862 end if
863
864 if (varNameToRead == varName .or. .not. containsFullField) then
865
866 if (statevector%dataKind == 4) then
867 field_r4_ptr(:,:, kIndex, stepIndex) = gd2d_file_r4(1:statevector%hco%ni, 1:statevector%hco%nj)
868 else
869 field_r8_ptr(:,:, kIndex, stepIndex) = real(gd2d_file_r4(1:statevector%hco%ni, 1:statevector%hco%nj), 8)
870 end if
871
872 else
873
874 select case (trim(varName))
875 case ('LVIS')
876
877 if (statevector%dataKind == 4) then
878 field_r4_ptr(:,:, kIndex, stepIndex) = &
879 log(max(min(gd2d_file_r4(1:statevector%hco%ni, 1:statevector%hco%nj), &
880 mpc_maximum_vis_r4), mpc_minimum_vis_r4))
881 else
882 field_r8_ptr(:,:, kIndex, stepIndex) = real(&
883 log(max(min(gd2d_file_r4(1:statevector%hco%ni, 1:statevector%hco%nj), &
884 mpc_maximum_vis_r4), mpc_minimum_vis_r4)), 8)
885 end if
886
887 case ('LPR')
888
889 if (statevector%dataKind == 4) then
890 field_r4_ptr(:,:, kIndex, stepIndex) = &
891 log(mpc_minimum_pr_r4 + max(0.0, gd2d_file_r4(1:statevector%hco%ni, &
892 1:statevector%hco%nj)))
893 else
894 field_r8_ptr(:,:, kIndex, stepIndex) = real(&
895 log(mpc_minimum_pr_r4 + max(0.0, gd2d_file_r4(1:statevector%hco%ni, &
896 1:statevector%hco%nj))), 8)
897 end if
898
899 case default
900 call utl_abort('gio_readFile: Oups! This should not happen... Check the code.')
901 end select
902 endif
903
904 if (ierr < 0)then
905 write(*,*) varNameToRead, ip1, datestamplist(stepIndex)
906 call utl_abort('gio_readFile: Problem with reading file')
907 end if
908
909 ! When mpi distribution could put UU on a different mpi task than VV
910 ! or only one wind component present in statevector
911 ! then we re-read the corresponding UV component and store it
912 if (statevector%extraUVallocated) then
913 if (varName == 'UU') then
914 ierr = fstlir(gd2d_file_r4(:,:),nulfile, ni_file, nj_file, nk_file, &
915 datestamplist(stepIndex), etiket_in, ip1, -1, -1, &
916 typvar_in, 'VV')
917
918 if (statevector%dataKind == 4) then
919 call gsv_getFieldUV(statevector, gd2d_r4_UV_ptr, kIndex)
920 gd2d_r4_UV_ptr(1:gsv_getHco(statevector)%ni, &
921 1:gsv_getHco(statevector)%nj, stepIndex) &
922 = gd2d_file_r4(1:gsv_getHco(statevector)%ni, &
923 1:gsv_getHco(statevector)%nj)
924 else
925 call gsv_getFieldUV(statevector, gd2d_r8_UV_ptr, kIndex)
926 gd2d_r8_UV_ptr(1:gsv_getHco(statevector)%ni, &
927 1:gsv_getHco(statevector)%nj, stepIndex) &
928 = real(gd2d_file_r4(1:gsv_getHco(statevector)%ni, &
929 1:gsv_getHco(statevector)%nj), 8)
930 end if
931
932 else if (varName == 'VV') then
933 ierr = fstlir(gd2d_file_r4(:,:), nulfile, ni_file, nj_file, nk_file, &
934 datestamplist(stepIndex), etiket_in, ip1, -1, -1, &
935 typvar_in, 'UU')
936
937 if (statevector%dataKind == 4) then
938 call gsv_getFieldUV(statevector, gd2d_r4_UV_ptr, kIndex)
939 gd2d_r4_UV_ptr(1:gsv_getHco(statevector)%ni, &
940 1:gsv_getHco(statevector)%nj, stepIndex) &
941 = gd2d_file_r4(1:gsv_getHco(statevector)%ni, &
942 1:gsv_getHco(statevector)%nj)
943 else
944 call gsv_getFieldUV(statevector, gd2d_r8_UV_ptr, kIndex)
945 gd2d_r8_UV_ptr(1:gsv_getHco(statevector)%ni, &
946 1:gsv_getHco(statevector)%nj, stepIndex) &
947 = real(gd2d_file_r4(1:gsv_getHco(statevector)%ni, &
948 1:gsv_getHco(statevector)%nj), 8)
949 end if
950
951 end if
952 end if
953
954 end do k_loop
955 end do
956
957 if (statevector%hco%global .and. statevector%mykCount > 0) then
958 write(*,*) 'deallocating hco_file'
959 call hco_deallocate(hco_file)
960 end if
961
962 ierr = fstfrm(nulfile)
963 ierr = fclos(nulfile)
964 if (associated(gd2d_file_r4)) deallocate(gd2d_file_r4)
965
966 ! Read in an oceanMask if it is present in the file
967 call gio_readMaskFromFile(statevector, trim(filename))
968
969 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
970 write(*,*) 'gio_readFile: finished'
971
972 end subroutine gio_readFile
973
974 !--------------------------------------------------------------------------
975 ! gio_readMaskFromFile
976 !--------------------------------------------------------------------------
977 subroutine gio_readMaskFromFile(stateVector, filename)
978 !
979 !:Purpose: Check if any ocean mask fields exist. If so, read for the surface
980 ! or all ocean depth levels.
981 !
982 implicit none
983
984 ! Arguments:
985 type(struct_gsv), intent(inout) :: stateVector
986 character(len=*), intent(in) :: fileName
987
988 call ocm_readMaskFromFile(stateVector%oceanMask,gsv_getHco(statevector), &
989 gsv_getVco(statevector), filename)
990
991 end subroutine gio_readMaskFromFile
992
993 !--------------------------------------------------------------------------
994 ! gio_getMaskLAM
995 !--------------------------------------------------------------------------
996 subroutine gio_getMaskLAM(statevector_mask, hco_ptr, vco_ptr, hInterpolateDegree_opt)
997 !
998 !:Purpose: To read a LAM mask from a file (./analinc_mask by default).
999 !
1000 implicit none
1001
1002 ! Arguments:
1003 type(struct_gsv), intent(inout) :: statevector_mask
1004 type(struct_hco), pointer, intent(in) :: hco_ptr
1005 type(struct_vco), pointer, intent(in) :: vco_ptr
1006 character(len=*), optional, intent(in) :: hInterpolateDegree_opt
1007
1008 ! Locals:
1009 character(len=12) :: hInterpolationDegree
1010
1011 if (present(hInterpolateDegree_opt)) then
1012 hInterpolationDegree = hInterpolateDegree_opt
1013 else
1014 hInterpolationDegree = 'LINEAR'
1015 end if
1016
1017 call gsv_allocate(statevector_mask, 1, hco_ptr, vco_ptr, dateStamp_opt=-1, &
1018 dataKind_opt=pre_incrReal, &
1019 mpi_local_opt=.true., varNames_opt=(/'MSKC'/), &
1020 hInterpolateDegree_opt=hInterpolationDegree)
1021 call gio_readFromFile(statevector_mask, './analinc_mask', ' ', ' ', unitConversion_opt=.false., &
1022 vcoFileIn_opt=vco_ptr)
1023
1024 end subroutine gio_getMaskLAM
1025
1026 !--------------------------------------------------------------------------
1027 ! gio_readTrials
1028 !--------------------------------------------------------------------------
1029 subroutine gio_readTrials(stateVectorTrialIn)
1030 !
1031 !:Purpose: Reading trials
1032 !
1033 implicit none
1034
1035 ! Arguments:
1036 type(struct_gsv), target, intent(inout) :: stateVectorTrialIn
1037
1038 ! Locals:
1039 logical :: fileExists, allocHeightSfc
1040 logical :: useInputStateVectorTrial
1041 integer, parameter :: maxNumTrials = 100
1042 integer :: fnom, fstouv, fclos, fstfrm, fstinf
1043 integer :: ierr, ikey, stepIndex, stepIndexToRead, trialIndex, nulTrial
1044 integer :: ni_file, nj_file, nk_file, dateStamp, varNameIndex
1045 integer :: procToRead, numBatch, batchIndex, stepIndexBeg, stepIndexEnd
1046 character(len=2) :: fileNumber
1047 character(len=512) :: fileName
1048 character(len=4) :: varNameForDateStampSearch
1049 character(len=4), pointer :: varNamesToRead(:)
1050 type(struct_gsv), target :: stateVectorTrial
1051 type(struct_gsv), pointer :: stateVectorTrial_ptr
1052 type(struct_gsv) :: stateVector_1step_r4
1053
1054 call utl_tmg_start(1,'--ReadTrials')
1055
1056 if ( mmpi_myid == 0 ) then
1057 write(*,*) ''
1058 write(*,*) 'gio_readTrials: START'
1059 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1060 end if
1061
1062 if ( gsv_varExist(stateVectorTrialIn,'Z_T') .or. &
1063 gsv_varExist(stateVectorTrialIn,'Z_M') .or. &
1064 gsv_varExist(stateVectorTrialIn,'P_T') .or. &
1065 gsv_varExist(stateVectorTrialIn,'P_M') ) then
1066
1067 useInputStateVectorTrial = .false.
1068
1069 allocHeightSfc = ( stateVectorTrialIn%vco%Vcode /= 0 )
1070
1071 ! Allocate single-precision statevector without Z/P to read trials
1072 call gsv_allocate( stateVectorTrial, stateVectorTrialIn%numStep, &
1073 stateVectorTrialIn%hco, stateVectorTrialIn%vco, &
1074 dateStamp_opt=tim_getDateStamp(), &
1075 mpi_local_opt=stateVectorTrialIn%mpi_local, &
1076 mpi_distribution_opt='Tiles', dataKind_opt=4, &
1077 allocHeightSfc_opt=allocHeightSfc, &
1078 hInterpolateDegree_opt=stateVectorTrialIn%hInterpolateDegree, &
1079 allocHeight_opt=.false., allocPressure_opt=.false., &
1080 beSilent_opt=.false. )
1081 call gsv_zero( stateVectorTrial )
1082 stateVectorTrial_ptr => stateVectorTrial
1083 else
1084 useInputStateVectorTrial = .true.
1085
1086 stateVectorTrial_ptr => stateVectorTrialIn
1087 end if
1088
1089 nullify(varNamesToRead)
1090 call gsv_varNamesList(varNamesToRead, stateVectorTrial_ptr)
1091
1092 varNameForDateStampSearch = ' '
1093 do varNameIndex = 1, size(varNamesToRead)
1094 select case (trim(varNamesToRead(varNameIndex)))
1095 case ('Z_T','Z_M','P_T','P_M')
1096 cycle
1097 case default
1098 varNameForDateStampSearch = varNamesToRead(varNameIndex)
1099 exit
1100 end select
1101 end do
1102
1103 ! warn if not enough mpi tasks
1104 if ( mmpi_nprocs < stateVectorTrial_ptr%numStep ) then
1105 write(*,*) 'gio_readTrials: number of trial time steps, mpi tasks = ', stateVectorTrial_ptr%numStep, mmpi_nprocs
1106 write(*,*) 'gio_readTrials: for better efficiency, the number of mpi tasks should '
1107 write(*,*) ' be at least as large as number of trial time steps'
1108 end if
1109
1110 allocHeightSfc = stateVectorTrial_ptr%heightSfcPresent
1111
1112 ! figure out number of batches of time steps for reading
1113 numBatch = ceiling(real(stateVectorTrial_ptr%numStep) / real(mmpi_nprocs))
1114 write(*,*) 'gio_readTrials: reading will be done by number of batches = ', numBatch
1115
1116 BATCH: do batchIndex = 1, numBatch
1117
1118 stepIndexBeg = 1 + (batchIndex - 1) * mmpi_nprocs
1119 stepIndexEnd = min(stateVectorTrial_ptr%numStep, stepIndexBeg + mmpi_nprocs - 1)
1120 write(*,*) 'gio_readTrials: batchIndex, stepIndexBeg/End = ', batchIndex, stepIndexBeg, stepIndexEnd
1121
1122 ! figure out which time step I will read, if any (-1 if none)
1123 stepIndexToRead = -1
1124 do stepIndex = stepIndexBeg, stepIndexEnd
1125 procToRead = nint( real(stepIndex - stepIndexBeg) * real(mmpi_nprocs) / real(stepIndexEnd - stepIndexBeg + 1) )
1126 if ( procToRead == mmpi_myid ) stepIndexToRead = stepIndex
1127 if ( mmpi_myid == 0 ) write(*,*) 'gio_readTrials: stepIndex, procToRead = ', stepIndex, procToRead
1128 end do
1129
1130 ! loop over all times for which stateVector is allocated
1131 if ( stepIndexToRead /= -1 ) then
1132 dateStamp = stateVectorTrial_ptr%dateStampList(stepIndexToRead)
1133 write(*,*) 'gio_readTrials: reading background for time step: ',stepIndexToRead, dateStamp
1134 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1135
1136 ! identify which trial file corresponds with current datestamp
1137 ikey = 0
1138 do trialIndex = 1, maxNumTrials
1139 write(fileNumber,'(I2.2)') trialIndex
1140 fileName = 'trlm_' // trim(fileNumber)
1141 inquire(file=trim(fileName),exist=fileExists)
1142 if ( .not. fileExists ) exit
1143 nulTrial = 0
1144 ierr = fnom(nulTrial,trim(fileName),'RND+OLD+R/O',0)
1145 ierr = fstouv(nulTrial,'RND+OLD')
1146 ikey = fstinf(nulTrial, ni_file, nj_file, nk_file, &
1147 dateStamp, ' ', -1, -1, -1, ' ', varNameForDateStampSearch)
1148 ierr = fstfrm(nulTrial)
1149 ierr = fclos(nulTrial)
1150 if ( ikey > 0 ) exit
1151 end do
1152
1153 if ( ikey <= 0 .or. .not.fileExists ) then
1154 write(*,*) 'stepIndexToRead, dateStamp = ', stepIndexToRead, dateStamp
1155 call utl_abort('gio_readTrials: trial file not found for this increment timestep')
1156 end if
1157
1158 ! allocate stateVector for storing just 1 time step
1159 if ( batchIndex == 1 ) then
1160 call gsv_allocate( stateVector_1step_r4, 1, stateVectorTrial_ptr%hco, stateVectorTrial_ptr%vco, &
1161 dateStamp_opt=dateStamp, mpi_local_opt=.false., dataKind_opt=4, &
1162 allocHeightSfc_opt=allocHeightSfc, varNames_opt=varNamesToRead, &
1163 hInterpolateDegree_opt=stateVectorTrial_ptr%hInterpolateDegree, &
1164 hExtrapolateDegree_opt=stateVectorTrial_ptr%hExtrapolateDegree)
1165 call gsv_zero( stateVector_1step_r4 )
1166 else
1167 call gsv_modifyDate( stateVector_1step_r4, dateStamp )
1168 end if
1169
1170 ! read the trial file for this timestep
1171 fileName = ram_fullWorkingPath(fileName)
1172 call gio_readFromFile(stateVector_1step_r4, fileName, ' ', 'P', &
1173 readHeightSfc_opt=allocHeightSfc)
1174
1175 ! remove from ram disk to save some space
1176 ierr = ram_remove(fileName)
1177 else
1178
1179 if ( gsv_isAllocated(stateVector_1step_r4) ) call gsv_deallocate(stateVector_1step_r4)
1180
1181 end if ! I read a time step
1182
1183 if ( stateVectorTrial_ptr%mpi_distribution == 'VarsLevs' ) then
1184 call gsv_transposeStepToVarsLevs(stateVector_1step_r4, stateVectorTrial_ptr, stepIndexBeg)
1185 else if ( stateVectorTrial_ptr%mpi_distribution == 'Tiles' ) then
1186 call gsv_transposeStepToTiles(stateVector_1step_r4, stateVectorTrial_ptr, stepIndexBeg)
1187 else
1188 call utl_abort( 'gio_readTrials: not compatible with mpi_distribution = ' // &
1189 trim(stateVectorTrial_ptr%mpi_distribution) )
1190 end if
1191
1192 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1193
1194 if ( gsv_isAllocated(stateVector_1step_r4) .and. batchIndex == numBatch ) then
1195 call gsv_deallocate(stateVector_1step_r4)
1196 end if
1197
1198 end do BATCH
1199
1200 if ( .not. useInputStateVectorTrial ) then
1201 call gsv_copy( stateVectorTrial_ptr, stateVectorTrialIn, allowVarMismatch_opt=.true. )
1202 call gsv_deallocate( stateVectorTrial )
1203 end if
1204
1205 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1206 write(*,*) 'gio_readTrials: FINISHED'
1207 write(*,*) ''
1208
1209 call utl_tmg_stop(1)
1210
1211 end subroutine gio_readTrials
1212
1213 !--------------------------------------------------------------------------
1214 ! gio_writeToFile
1215 !--------------------------------------------------------------------------
1216 subroutine gio_writeToFile(statevector_in, fileName, etiket_in, &
1217 scaleFactor_opt, ip3_opt, stepIndex_opt, typvar_opt,&
1218 HUcontainsLQ_opt, unitConversion_opt, &
1219 writeHeightSfc_opt, numBits_opt, containsFullField_opt)
1220 !
1221 !:Purpose: Write a statevector object to an RPN standard file.
1222 !
1223 implicit none
1224
1225 ! Arguments:
1226 type(struct_gsv), target, intent(in) :: statevector_in
1227 character(len=*), intent(in) :: fileName
1228 character(len=*), intent(in) :: etiket_in
1229 real(8), optional, intent(in) :: scaleFactor_opt
1230 integer, optional, intent(in) :: ip3_opt, stepIndex_opt
1231 character(len=*), optional, intent(in) :: typvar_opt
1232 logical, optional, intent(in) :: HUcontainsLQ_opt
1233 logical, optional, intent(in) :: unitConversion_opt
1234 logical, optional, intent(in) :: writeHeightSfc_opt
1235 integer, optional, intent(in) :: numBits_opt
1236 logical, optional, intent(in) :: containsFullField_opt
1237
1238 ! Locals:
1239 logical :: iDoWriting, unitConversion, containsFullField
1240 integer :: fclos, fnom, fstouv, fstfrm
1241 integer :: nulfile, stepIndex
1242 integer :: ierr, fstecr, ezdefset
1243 integer :: ni, nj, nk
1244 integer :: dateo, npak, levIndex, nlev, varIndex, maskLevIndex
1245 integer :: ip1, ip2, ip3, deet, npas, datyp
1246 integer :: ig1 ,ig2 ,ig3 ,ig4
1247 integer :: yourid, nsize, youridy, youridx
1248 real(4) :: factor_r4, work_r4
1249 character(len=1) :: grtyp
1250 character(len=4) :: nomvar
1251 character(len=2) :: typvar
1252 character(len=12) :: etiket
1253 character(len=4) :: varLevel
1254 character(len=4), pointer :: varNamesToRead(:)
1255 integer, allocatable :: mask(:,:)
1256 real(4), allocatable :: work2d_r4(:,:), work2dFile_r4(:,:), gd_send_r4(:,:), gd_recv_r4(:,:,:)
1257 real(8), pointer :: field_r8(:,:,:,:), heightSfc_ptr(:,:)
1258 real(4), pointer :: field_r4(:,:,:,:)
1259 type(struct_gsv), pointer :: statevector
1260 type(struct_gsv), target :: statevector_tiles
1261
1262 write(*,*) 'gio_writeToFile: START'
1263
1264 call utl_tmg_start(161,'low-level--gsv_writeToFile')
1265
1266 call readNml()
1267
1268 !
1269 !- 1. Since this routine can only work with 'Tiles' distribution when mpi_local = .true.,
1270 ! transpose a statevector using 'VarsLevs' distribution
1271 !
1272 if ( stateVector_in%mpi_distribution == 'VarsLevs' .and. &
1273 stateVector_in%mpi_local ) then
1274 nullify(varNamesToRead)
1275 call gsv_varNamesList(varNamesToRead,statevector_in)
1276 call gsv_allocate(statevector_tiles, statevector_in%numStep, statevector_in%hco, &
1277 statevector_in%vco, dataKind_opt=statevector_in%dataKind, &
1278 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
1279 dateStampList_opt=statevector_in%dateStampList, &
1280 varNames_opt=varNamesToRead)
1281 call gsv_transposeVarsLevsToTiles(statevector_in, statevector_tiles)
1282 statevector => stateVector_tiles
1283 deallocate(varNamesToRead)
1284 else
1285 statevector => stateVector_in
1286 end if
1287
1288 !
1289 !- 2. Set some variables
1290 !
1291 if ( .not. statevector%mpi_local ) then
1292 write(*,*) 'gio_writeToFile: writing statevector that is already mpiglobal!'
1293 end if
1294
1295 if (present(ip3_opt)) then
1296 ip3 = ip3_opt
1297 else
1298 ip3 = 0
1299 end if
1300
1301 if (present(unitConversion_opt)) then
1302 unitConversion = unitConversion_opt
1303 else
1304 unitConversion = .true.
1305 end if
1306
1307 if (present(containsFullField_opt)) then
1308 containsFullField = containsFullField_opt
1309 else
1310 containsFullField = .false.
1311 end if
1312 write(*,*) 'gio_writeToFile: containsFullField = ', containsFullField
1313
1314 ! if step index not specified, choose anltime (usually center of window)
1315 if (present(stepIndex_opt)) then
1316 stepIndex = stepIndex_opt
1317 else
1318 stepIndex = statevector%anltime
1319 end if
1320
1321 if ( present(numBits_opt) ) then
1322 npak = -numBits_opt
1323 else
1324 npak = -32
1325 end if
1326
1327 ! initialization of parameters for writing to file
1328 if (statevector%dateOriginList(stepIndex) /= mpc_missingValue_int) then
1329 dateo = statevector%dateOriginList(stepIndex)
1330 else
1331 dateo = statevector%dateStampList(stepIndex)
1332 end if
1333
1334 if (statevector%deet /= mpc_missingValue_int) then
1335 deet = statevector%deet
1336 else
1337 deet = 0
1338 end if
1339
1340 if (statevector%npasList(stepIndex) /= mpc_missingValue_int) then
1341 npas = statevector%npasList(stepIndex)
1342 else
1343 npas = 0
1344 end if
1345
1346 if (statevector%ip2List(stepIndex) /= mpc_missingValue_int) then
1347 ip2 = statevector%ip2List(stepIndex)
1348 else
1349 ip2 = 0
1350 end if
1351
1352 if (statevector%etiket /= 'UNDEFINED') then
1353 etiket = statevector%etiket
1354 else
1355 etiket = trim(etiket_in)
1356 end if
1357
1358 ni = statevector%ni
1359 nj = statevector%nj
1360 nk = 1
1361 if ( present(typvar_opt) ) then
1362 typvar = trim(typvar_opt)
1363 else
1364 typvar = 'R'
1365 end if
1366 if ( statevector%oceanMask%maskPresent ) then
1367 typvar(2:2) = '@'
1368 end if
1369 grtyp = statevector%hco%grtyp
1370 ig1 = statevector%hco%ig1
1371 ig2 = statevector%hco%ig2
1372 ig3 = statevector%hco%ig3
1373 ig4 = statevector%hco%ig4
1374 datyp = 134
1375
1376 ! only proc 0 does writing or each proc when data is global
1377 ! (assuming only called for proc with global data)
1378 iDoWriting = (mmpi_myid == 0) .or. (.not. statevector%mpi_local)
1379
1380 !
1381 !- 3. Write the global StateVector
1382 !
1383 if (iDoWriting) then
1384
1385 !- Open output field
1386 nulfile = 0
1387 write(*,*) 'gio_writeToFile: file name = ',trim(fileName)
1388 ierr = fnom(nulfile,trim(fileName),'RND+APPEND',0)
1389
1390 if ( ierr >= 0 ) then
1391 ierr = fstouv(nulfile,'RND')
1392 else
1393 call utl_abort('gio_writeToFile: problem opening output file')
1394 end if
1395
1396 if (nulfile == 0 ) then
1397 call utl_abort('gio_writeToFile: unit number for output file not valid')
1398 end if
1399
1400 !- Write TicTacToc
1401 if ( (mmpi_myid == 0 .and. statevector%mpi_local) .or. .not.statevector%mpi_local ) then
1402 call writeTicTacToc(statevector,nulfile,etiket) ! IN
1403 endif
1404
1405 end if
1406
1407 allocate(gd_send_r4(statevector%lonPerPEmax,statevector%latPerPEmax))
1408 if ( mmpi_myid == 0 .or. (.not. statevector%mpi_local) ) then
1409 allocate(work2d_r4(statevector%ni,statevector%nj))
1410 if (statevector%mpi_local) then
1411 ! Receive tile data from all mpi tasks
1412 allocate(gd_recv_r4(statevector%lonPerPEmax,statevector%latPerPEmax,mmpi_nprocs))
1413 else
1414 ! Already have entire domain on mpi task (lat/lonPerPEmax == nj/ni)
1415 allocate(gd_recv_r4(statevector%lonPerPEmax,statevector%latPerPEmax,1))
1416 end if
1417 else
1418 allocate(gd_recv_r4(1,1,1))
1419 allocate(work2d_r4(1,1))
1420 end if
1421
1422 ! Write surface height, if requested
1423 if ( present(writeHeightSfc_opt) ) then
1424 if ( writeHeightSfc_opt .and. gsv_isAssocHeightSfc(statevector) ) then
1425 write(*,*) 'gio_writeToFile: writing surface height'
1426 heightSfc_ptr => gsv_getHeightSfc(statevector)
1427 ! MPI communication
1428 gd_send_r4(1:statevector%lonPerPE, &
1429 1:statevector%latPerPE) = &
1430 real(heightSfc_ptr(statevector%myLonBeg:statevector%myLonEnd, &
1431 statevector%myLatBeg:statevector%myLatEnd),4)
1432 if ( (mmpi_nprocs > 1) .and. statevector%mpi_local ) then
1433 nsize = statevector%lonPerPEmax * statevector%latPerPEmax
1434 call rpn_comm_gather(gd_send_r4, nsize, 'mpi_real4', &
1435 gd_recv_r4, nsize, 'mpi_real4', 0, 'grid', ierr )
1436 else
1437 ! just copy when either nprocs is 1 or data is global
1438 gd_recv_r4(:,:,1) = gd_send_r4(:,:)
1439 end if
1440 if ( mmpi_myid == 0 .and. statevector%mpi_local ) then
1441 !$OMP PARALLEL DO PRIVATE(youridy,youridx,yourid)
1442 do youridy = 0, (mmpi_npey-1)
1443 do youridx = 0, (mmpi_npex-1)
1444 yourid = youridx + youridy*mmpi_npex
1445 work2d_r4(statevector%allLonBeg(youridx+1):statevector%allLonEnd(youridx+1), &
1446 statevector%allLatBeg(youridy+1):statevector%allLatEnd(youridy+1)) = &
1447 gd_recv_r4(1:statevector%allLonPerPE(youridx+1), &
1448 1:statevector%allLatPerPE(youridy+1),yourid+1)
1449 end do
1450 end do
1451 !$OMP END PARALLEL DO
1452 else if ( .not. statevector%mpi_local ) then
1453 work2d_r4(:,:) = gd_recv_r4(:,:,1)
1454 end if
1455
1456 ! now do writing
1457 if (iDoWriting) then
1458 ip1 = statevector%vco%ip1_sfc
1459 nomvar = 'GZ'
1460
1461 !- Scale
1462 factor_r4 = real(1.0d0/10.0d0,4)
1463 work2d_r4(:,:) = factor_r4 * work2d_r4(:,:)
1464
1465 !- Writing to file
1466 ierr = fstecr(work2d_r4, work_r4, npak, nulfile, dateo, deet, npas, ni, nj, &
1467 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
1468 ig1, ig2, ig3, ig4, datyp, .false.)
1469 end if ! iDoWriting
1470
1471 end if
1472 end if
1473
1474 do varIndex = 1, vnl_numvarmax
1475
1476 if (gsv_varExist(statevector,vnl_varNameList(varIndex)) ) then
1477
1478 nlev = statevector%varNumLev(varIndex)
1479
1480 do levIndex = 1, nlev
1481
1482 if ( statevector%dataKind == 8 ) then
1483 call gsv_getField(statevector,field_r8,vnl_varNameList(varIndex))
1484 gd_send_r4(1:statevector%lonPerPE, &
1485 1:statevector%latPerPE) = &
1486 real(field_r8(statevector%myLonBeg:statevector%myLonEnd, &
1487 statevector%myLatBeg:statevector%myLatEnd,levIndex,stepIndex),4)
1488 else
1489 call gsv_getField(statevector,field_r4,vnl_varNameList(varIndex))
1490 gd_send_r4(1:statevector%lonPerPE, &
1491 1:statevector%latPerPE) = &
1492 field_r4(statevector%myLonBeg:statevector%myLonEnd, &
1493 statevector%myLatBeg:statevector%myLatEnd,levIndex,stepIndex)
1494 end if
1495
1496 nsize = statevector%lonPerPEmax*statevector%latPerPEmax
1497 if ( (mmpi_nprocs > 1) .and. (statevector%mpi_local) ) then
1498 call rpn_comm_gather(gd_send_r4, nsize, 'mpi_real4', &
1499 gd_recv_r4, nsize, 'mpi_real4', 0, 'grid', ierr )
1500 else
1501 ! just copy when either nprocs is 1 or data is global
1502 gd_recv_r4(:,:,1) = gd_send_r4(:,:)
1503 end if
1504
1505 if ( mmpi_myid == 0 .and. statevector%mpi_local ) then
1506 !$OMP PARALLEL DO PRIVATE(youridy,youridx,yourid)
1507 do youridy = 0, (mmpi_npey-1)
1508 do youridx = 0, (mmpi_npex-1)
1509 yourid = youridx + youridy*mmpi_npex
1510 work2d_r4(statevector%allLonBeg(youridx+1):statevector%allLonEnd(youridx+1), &
1511 statevector%allLatBeg(youridy+1):statevector%allLatEnd(youridy+1)) = &
1512 gd_recv_r4(1:statevector%allLonPerPE(youridx+1), &
1513 1:statevector%allLatPerPE(youridy+1), yourid+1)
1514 end do
1515 end do
1516 !$OMP END PARALLEL DO
1517 else if ( .not. statevector%mpi_local ) then
1518 work2d_r4(:,:) = gd_recv_r4(:,:,1)
1519 end if
1520
1521 ! now do writing
1522 if (iDoWriting) then
1523
1524 ! Set the ip1 value
1525 if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'MM') then
1526 ip1 = statevector%vco%ip1_M(levIndex)
1527 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'TH') then
1528 ip1 = statevector%vco%ip1_T(levIndex)
1529 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'SF') then
1530 ip1 = 0
1531 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'SFTH') then
1532 ip1 = statevector%vco%ip1_T_2m
1533 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'SFMM') then
1534 ip1 = statevector%vco%ip1_M_10m
1535 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'OT') then
1536 ip1 = vco_ip1_other(levIndex)
1537 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'DP') then
1538 ip1 = statevector%vco%ip1_depth(levIndex)
1539 else if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'SS') then
1540 ip1 = statevector%vco%ip1_seaLevel
1541 else
1542 varLevel = vnl_varLevelFromVarname(vnl_varNameList(varIndex))
1543 write(*,*) 'gio_writeToFile: unknown type of vertical level: ', varLevel
1544 call utl_abort('gio_writeToFile')
1545 end if
1546
1547 ! Set the level index for the mask (if present)
1548 if (vnl_varLevelFromVarname(vnl_varNameList(varIndex)) == 'DP') then
1549 maskLevIndex = levIndex
1550 else
1551 maskLevIndex = 1
1552 end if
1553
1554 ! Set the output variable name
1555 nomvar = trim(vnl_varNameList(varIndex))
1556 if ( trim(nomvar) == 'HU' .and. present(HUcontainsLQ_opt) ) then
1557 if ( HUcontainsLQ_opt ) nomvar = 'LQ'
1558 end if
1559
1560 if ( vnl_varKindFromVarname(trim(nomvar)) == 'CH' .and. containsFullField ) then
1561 ! Impose lower limits
1562 if ( gsv_minValVarKindCH(vnl_varListIndex(nomvar)) > 1.01*mpc_missingValue_r8 ) &
1563 work2d_r4(:,:) = max( work2d_r4(:,:), real(gsv_minValVarKindCH(vnl_varListIndex(trim(nomvar)))) )
1564 end if
1565
1566 ! Set the conversion factor
1567 if ( unitConversion ) then
1568
1569 if ( trim(nomvar) == 'UU' .or. trim(nomvar) == 'VV') then
1570 factor_r4 = mpc_knots_per_m_per_s_r4 ! m/s -> knots
1571 else if ( trim(nomvar) == 'P0' .or. trim(nomvar) == 'UP' .or. &
1572 trim(nomvar) == 'PB' .or. trim(nomvar) == 'P0LS' ) then
1573 factor_r4 = 0.01 ! Pa -> hPa
1574 else if ( vnl_varKindFromVarname(trim(nomvar)) == 'CH' ) then
1575 if ( gsv_conversionVarKindCHtoMicrograms ) then
1576 ! Apply inverse transform of unit conversion
1577 if ( trim(nomvar) == 'TO3' .or. trim(nomvar) == 'O3L' ) then
1578 factor_r4 = 1.0E-9*mpc_molar_mass_dry_air_r4 &
1579 /vnl_varMassFromVarName(trim(nomvar)) ! micrograms/kg -> vmr
1580 else
1581 factor_r4 = 1.0d0 ! no conversion
1582 end if
1583 else
1584 factor_r4 = 1.0d0 ! no conversion
1585 end if
1586 else
1587 factor_r4 = 1.0d0 ! no conversion
1588 end if
1589 else
1590 factor_r4 = 1.0
1591 end if
1592
1593 if (present(scaleFactor_opt)) factor_r4 = factor_r4 * real(scaleFactor_opt,4)
1594
1595 !- Scale
1596 work2d_r4(:,:) = factor_r4 * work2d_r4(:,:)
1597
1598 !- Convert Kelvin to Celcius only if full field
1599 if (containsFullField .and. (trim(nomvar) == 'TT' .or. trim(nomvar) == 'TM') ) then
1600 where (work2d_r4(:,:) > 100.0)
1601 work2d_r4(:,:) = work2d_r4(:,:) - mpc_k_c_degree_offset_r4
1602 end where
1603 end if
1604
1605 !- Do interpolation back to physics grid, if needed
1606 if ( interpToPhysicsGrid .and. statevector%onPhysicsGrid(varIndex) ) then
1607 write(*,*) 'writeToFile: interpolate this variable back to physics grid: ', &
1608 nomvar, associated(statevector%hco_physics)
1609 allocate(work2dFile_r4(statevector%hco_physics%ni,statevector%hco_physics%nj))
1610 work2dFile_r4(:,:) = 0.0
1611 ierr = ezdefset( statevector%hco_physics%EZscintID, statevector%hco%EZscintID )
1612 ierr = int_hInterpScalar( work2dFile_r4, work2d_r4, &
1613 interpDegree='NEAREST', extrapDegree_opt='NEUTRAL' )
1614
1615 !- Writing to file
1616 ierr = fstecr(work2dFile_r4, work_r4, npak, nulfile, dateo, deet, npas, &
1617 statevector%hco_physics%ni, statevector%hco_physics%nj, &
1618 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
1619 statevector%hco_physics%ig1, statevector%hco_physics%ig2, &
1620 statevector%hco_physics%ig3, statevector%hco_physics%ig4, &
1621 datyp, .false.)
1622 deallocate(work2dFile_r4)
1623
1624 else
1625
1626 !- Writing to file
1627 ierr = fstecr(work2d_r4, work_r4, npak, nulfile, dateo, deet, npas, ni, nj, &
1628 nk, ip1, ip2, ip3, typvar, nomvar, etiket, grtyp, &
1629 ig1, ig2, ig3, ig4, datyp, .false.)
1630
1631 end if
1632
1633 if ( statevector%oceanMask%maskPresent ) then
1634 if (.not.allocated(mask)) allocate(mask(ni,nj))
1635 call ocm_copyToInt(statevector%oceanMask,mask,maskLevIndex)
1636 ierr = fstecr(mask, work_r4, -1, nulfile, dateo, deet, npas, ni, nj, &
1637 nk, ip1, ip2, ip3, '@@', nomvar, etiket, grtyp, &
1638 ig1, ig2, ig3, ig4, 2, .false.)
1639 end if
1640
1641 end if ! iDoWriting
1642
1643 end do ! levIndex
1644
1645 end if ! varExist
1646
1647 end do ! varIndex
1648
1649 deallocate(work2d_r4)
1650 deallocate(gd_send_r4)
1651 deallocate(gd_recv_r4)
1652 if (allocated(mask)) deallocate(mask)
1653
1654 if (iDoWriting) then
1655 ierr = fstfrm(nulfile)
1656 ierr = fclos(nulfile)
1657 end if
1658
1659 !
1660 !- 4. Ending
1661 !
1662 if ( stateVector_in%mpi_distribution == 'VarsLevs' .and. &
1663 stateVector_in%mpi_local ) then
1664 call gsv_deallocate(statevector_tiles)
1665 end if
1666
1667 call utl_tmg_stop(161)
1668 write(*,*) 'gio_writeToFile: END'
1669
1670 end subroutine gio_writeToFile
1671
1672 !--------------------------------------------------------------------------
1673 ! writeTicTacToc
1674 !--------------------------------------------------------------------------
1675 subroutine writeTicTacToc(statevector,iun,etiket)
1676 !
1677 !:Purpose: Write a statevector object grid descriptors to an RPN standard
1678 ! file.
1679 !
1680 implicit none
1681
1682 ! Arguments:
1683 type(struct_gsv), intent(in) :: statevector
1684 integer, intent(in) :: iun
1685 character(len=*), intent(in) :: etiket
1686
1687 ! Locals:
1688 integer :: ier
1689 integer :: dateo, npak, status, fstecr
1690 integer :: ip1,ip2,ip3,deet,npas,datyp,ig1,ig2,ig3,ig4
1691 integer :: ig1_tictac,ig2_tictac,ig3_tictac,ig4_tictac
1692 character(len=1) :: grtyp
1693 character(len=2) :: typvar
1694
1695 !
1696 !- 1. Writing Tic-Tac
1697 !
1698 if ( statevector % hco % grtyp == 'Z' ) then
1699 npak = -32
1700 deet = 0
1701 ip1 = statevector%hco%ig1
1702 ip2 = statevector%hco%ig2
1703 ip3 = statevector%hco%ig3
1704 npas = 0
1705 datyp = 1
1706 grtyp = statevector%hco%grtypTicTac
1707 typvar = 'X'
1708 dateo = 0
1709
1710 call cxgaig ( grtyp, & ! IN
1711 ig1_tictac, ig2_tictac, ig3_tictac, ig4_tictac, & ! OUT
1712 real(statevector%hco%xlat1), real(statevector%hco%xlon1), & ! IN
1713 real(statevector%hco%xlat2), real(statevector%hco%xlon2) ) ! IN
1714
1715 ig1 = ig1_tictac
1716 ig2 = ig2_tictac
1717 ig3 = ig3_tictac
1718 ig4 = ig4_tictac
1719
1720 ier = utl_fstecr(statevector%hco%lon*mpc_degrees_per_radian_r8, npak, &
1721 iun, dateo, deet, npas, statevector%ni, 1, 1, ip1, &
1722 ip2, ip3, typvar, '>>', etiket, grtyp, ig1, &
1723 ig2, ig3, ig4, datyp, .true.)
1724
1725 ier = utl_fstecr(statevector%hco%lat*mpc_degrees_per_radian_r8, npak, &
1726 iun, dateo, deet, npas, 1, statevector%nj, 1, ip1, &
1727 ip2, ip3, typvar, '^^', etiket, grtyp, ig1, &
1728 ig2, ig3, ig4, datyp, .true.)
1729
1730 ! Also write the tic tac for the physics grid
1731 if ( any(statevector%onPhysicsGrid(:)) ) then
1732
1733 ip1 = statevector%hco_physics%ig1
1734 ip2 = statevector%hco_physics%ig2
1735 ip3 = statevector%hco_physics%ig3
1736 grtyp = statevector%hco_physics%grtypTicTac
1737
1738 call cxgaig ( grtyp, & ! IN
1739 ig1_tictac, ig2_tictac, ig3_tictac, ig4_tictac, & ! OUT
1740 real(statevector%hco_physics%xlat1), real(statevector%hco_physics%xlon1), & ! IN
1741 real(statevector%hco_physics%xlat2), real(statevector%hco_physics%xlon2) ) ! IN
1742
1743 ig1 = ig1_tictac
1744 ig2 = ig2_tictac
1745 ig3 = ig3_tictac
1746 ig4 = ig4_tictac
1747
1748 ier = utl_fstecr(statevector%hco_physics%lon*mpc_degrees_per_radian_r8, npak, &
1749 iun, dateo, deet, npas, statevector%hco_physics%ni, 1, 1, ip1, &
1750 ip2, ip3, typvar, '>>', etiket, grtyp, ig1, &
1751 ig2, ig3, ig4, datyp, .true.)
1752
1753 ier = utl_fstecr(statevector%hco_physics%lat*mpc_degrees_per_radian_r8, npak, &
1754 iun, dateo, deet, npas, 1, statevector%hco_physics%nj, 1, ip1, &
1755 ip2, ip3, typvar, '^^', etiket, grtyp, ig1, &
1756 ig2, ig3, ig4, datyp, .true.)
1757 end if
1758
1759 else if ( statevector % hco % grtyp == 'U' ) then
1760 npak = -32
1761 ier = fstecr(statevector%hco%tictacU, statevector%hco%tictacU, npak, iun, 0, 0, 0, size(statevector%hco%tictacU), 1, 1 , &
1762 statevector%hco%ig1, statevector%hco%ig2, statevector%hco%ig3, 'X', '^>', etiket, &
1763 'F', 1, 0, 0, 0, 5, .false.)
1764
1765 else if ( statevector % hco % grtyp == 'Y' ) then
1766 npak = -32
1767 deet = 0
1768 ip1 = statevector%hco%ig1
1769 ip2 = statevector%hco%ig2
1770 ip3 = statevector%hco%ig3
1771 npas = 0
1772 datyp = 1
1773 grtyp = statevector%hco%grtypTicTac
1774 typvar = 'X'
1775 dateo = 0
1776
1777 call cxgaig ( grtyp, & ! IN
1778 ig1_tictac, ig2_tictac, ig3_tictac, ig4_tictac, & ! OUT
1779 real(statevector%hco%xlat1), real(statevector%hco%xlon1), & ! IN
1780 real(statevector%hco%xlat2), real(statevector%hco%xlon2) ) ! IN
1781
1782 ig1 = ig1_tictac
1783 ig2 = ig2_tictac
1784 ig3 = ig3_tictac
1785 ig4 = ig4_tictac
1786
1787 ier = utl_fstecr(statevector%hco%lon2d_4*mpc_degrees_per_radian_r8, npak, &
1788 iun, dateo, deet, npas, statevector%ni, statevector%nj, 1, &
1789 ip1, ip2, ip3, typvar, '>>', etiket, grtyp, &
1790 ig1, ig2, ig3, ig4, datyp, .true.)
1791
1792 ier = utl_fstecr(statevector%hco%lat2d_4*mpc_degrees_per_radian_r8, npak, &
1793 iun, dateo, deet, npas, statevector%ni, statevector%nj, 1, &
1794 ip1, ip2, ip3, typvar, '^^', etiket, grtyp, &
1795 ig1, ig2, ig3, ig4, datyp, .true.)
1796
1797
1798 end if
1799
1800 !
1801 !- Writing Toc-Toc
1802 !
1803 if ( statevector%vco%vgridPresent ) then
1804 status = vgd_write(statevector%vco%vgrid,iun,'fst')
1805 if ( status /= VGD_OK ) then
1806 call utl_abort('writeTicTacToc: ERROR with vgd_write')
1807 end if
1808 end if
1809
1810 end subroutine writeTicTacToc
1811
1812 !--------------------------------------------------------------------------
1813 ! gio_fileUnitsToStateUnits
1814 !--------------------------------------------------------------------------
1815 subroutine gio_fileUnitsToStateUnits(statevector, containsFullField, stepIndex_opt)
1816 !
1817 !:Purpose: Unit conversion needed after reading RPN standard file
1818 !
1819 implicit none
1820
1821 ! Arguments:
1822 type(struct_gsv), intent(inout) :: statevector
1823 logical, intent(in) :: containsFullField
1824 integer, optional, intent(in) :: stepIndex_opt
1825
1826 ! Locals:
1827 integer :: stepIndex, stepIndexBeg, stepIndexEnd, kIndex
1828 real(4), pointer :: field_r4_ptr(:,:,:,:), fieldUV_r4_ptr(:,:,:)
1829 real(8), pointer :: field_r8_ptr(:,:,:,:), fieldUV_r8_ptr(:,:,:)
1830 real(8) :: multFactor
1831 character(len=4) :: varName
1832
1833 if (present(stepIndex_opt)) then
1834 stepIndexBeg = stepIndex_opt
1835 stepIndexEnd = stepIndex_opt
1836 else
1837 stepIndexBeg = 1
1838 stepIndexEnd = statevector%numStep
1839 end if
1840
1841 if (statevector%dataKind == 4) then
1842 call gsv_getField(statevector, field_r4_ptr)
1843 else
1844 call gsv_getField(statevector, field_r8_ptr)
1845 end if
1846
1847 step_loop: do stepIndex = stepIndexBeg, stepIndexEnd
1848
1849 ! Do unit conversion for all variables
1850 KINDEXCYCLE: do kIndex = statevector%mykBeg, statevector%mykEnd
1851 varName = gsv_getVarNameFromK(statevector, kIndex)
1852
1853 if (trim(varName) == 'UU' .or. trim(varName) == 'VV') then
1854 multFactor = mpc_m_per_s_per_knot_r8 ! knots -> m/s
1855 else if (trim(varName) == 'P0' .or. trim(varName) == 'P0LS') then
1856 multFactor = mpc_pa_per_mbar_r8 ! hPa -> Pa
1857 else if (vnl_varKindFromVarname(trim(varName)) == 'CH') then
1858 if (gsv_conversionVarKindCHtoMicrograms) then
1859 if (trim(varName) == 'TO3' .or. trim(varName) == 'O3L') then
1860 ! Convert from volume mixing ratio to micrograms/kg
1861 ! Standard ozone input would not require this conversion as it is already in micrograms/kg
1862 multFactor = 1.0d9 * vnl_varMassFromVarName(trim(varName)) / &
1863 mpc_molar_mass_dry_air_r8 ! vmr -> micrograms/kg
1864 else
1865 multFactor = 1.0d0 ! no conversion
1866 end if
1867 else
1868 multFactor = 1.0d0 ! no conversion
1869 end if
1870 else
1871 multFactor = 1.0d0 ! no conversion
1872 end if
1873
1874 if (multFactor /= 1.0d0) then
1875 if (statevector%dataKind == 4) then
1876 field_r4_ptr(:,:, kIndex, stepIndex) = real(multFactor * field_r4_ptr(:,:, kIndex, stepIndex), 4)
1877 else
1878 field_r8_ptr(:,:, kIndex, stepIndex) = real(multFactor * field_r8_ptr(:,:, kIndex, stepIndex), 8)
1879 end if
1880 end if
1881
1882 if (trim(varName) == 'TT' .and. containsFullField) then
1883 if (statevector%dataKind == 4) then
1884 field_r4_ptr(:,:, kIndex, stepIndex) = real(field_r4_ptr(:,:, kIndex, stepIndex) + &
1885 mpc_k_c_degree_offset_r8, 4)
1886 else
1887 field_r8_ptr(:,:, kIndex, stepIndex) = real(field_r8_ptr(:,:, kIndex, stepIndex) + &
1888 mpc_k_c_degree_offset_r8, 8)
1889 end if
1890 end if
1891
1892 if (trim(varName) == 'TM' .and. containsFullField) then
1893 if (statevector%dataKind == 4) then
1894 where (field_r4_ptr(:,:, kIndex, stepIndex) < 100.0)
1895 field_r4_ptr(:,:, kIndex, stepIndex) = real(field_r4_ptr(:,:, kIndex, stepIndex) + &
1896 mpc_k_c_degree_offset_r8, 4)
1897 end where
1898 else
1899 where (field_r8_ptr(:,:, kIndex, stepIndex) < 100.0)
1900 field_r8_ptr(:,:, kIndex, stepIndex) = real(field_r8_ptr(:,:, kIndex, stepIndex) + &
1901 mpc_k_c_degree_offset_r8, 8)
1902 end where
1903 end if
1904 end if
1905
1906 if (trim(varName) == 'VIS' .and. containsFullField) then
1907 if (statevector%dataKind == 4) then
1908 field_r4_ptr(:,:, kIndex, stepIndex) = min(field_r4_ptr(:,:, kIndex, stepIndex), mpc_maximum_vis_r4)
1909 else
1910 field_r8_ptr(:,:, kIndex, stepIndex) = min(field_r8_ptr(:,:, kIndex, stepIndex), mpc_maximum_vis_r8)
1911 end if
1912 end if
1913
1914 if (vnl_varKindFromVarname(trim(varName)) == 'CH' .and. containsFullField) then
1915 if (gsv_minValVarKindCH(vnl_varListIndex(varName)) > 1.01 * mpc_missingValue_r8) then
1916 if (statevector%dataKind == 4) then
1917 field_r4_ptr(:,:, kIndex, stepIndex) = max(field_r4_ptr(:,:, kIndex, stepIndex), &
1918 real(gsv_minValVarKindCH(vnl_varListIndex(trim(varName)))))
1919 else
1920 field_r8_ptr(:,:, kIndex, stepIndex) = max(field_r8_ptr(:,:, kIndex, stepIndex), &
1921 real(gsv_minValVarKindCH(vnl_varListIndex(trim(varName)))))
1922 end if
1923 end if
1924 end if
1925
1926 if (trim(varName) == 'PR' .and. containsFullField) then
1927 if (statevector%dataKind == 4) then
1928 field_r4_ptr(:,:, kIndex, stepIndex) = max(field_r4_ptr(:,:, kIndex, stepIndex), 0.0)
1929 else
1930 field_r8_ptr(:,:, kIndex, stepIndex) = max(field_r8_ptr(:,:, kIndex, stepIndex), 0.0)
1931 end if
1932 end if
1933 end do KINDEXCYCLE
1934
1935 ! Do unit conversion for extra copy of winds, if present
1936 IFWINDS: if (statevector%extraUVallocated) then
1937 multFactor = mpc_m_per_s_per_knot_r8 ! knots -> m/s
1938
1939 if (statevector%dataKind == 4) then
1940 !$OMP PARALLEL DO PRIVATE (kIndex, fieldUV_r4_ptr)
1941 do kIndex = statevector%myUVkBeg, statevector%myUVkEnd
1942 nullify(fieldUV_r4_ptr)
1943 call gsv_getFieldUV(statevector, fieldUV_r4_ptr, kIndex)
1944 fieldUV_r4_ptr(:,:, stepIndex) = real(multFactor * fieldUV_r4_ptr(:,:, stepIndex), 4)
1945 end do
1946 !$OMP END PARALLEL DO
1947 else
1948 !$OMP PARALLEL DO PRIVATE (kIndex, fieldUV_r8_ptr)
1949 do kIndex = statevector%myUVkBeg, statevector%myUVkEnd
1950 nullify(fieldUV_r8_ptr)
1951 call gsv_getFieldUV(statevector, fieldUV_r8_ptr, kIndex)
1952 fieldUV_r8_ptr(:,:, stepIndex) = real(multFactor * fieldUV_r8_ptr(:,:, stepIndex), 8)
1953 end do
1954 !$OMP END PARALLEL DO
1955 end if
1956 end if IFWINDS
1957
1958 end do step_loop
1959
1960 end subroutine gio_fileUnitsToStateUnits
1961
1962 !--------------------------------------------------------------------------
1963 ! readNml (private)
1964 !--------------------------------------------------------------------------
1965 subroutine readNml()
1966 !
1967 !:Purpose: Read the namelist NAMSTIO
1968 !
1969 implicit none
1970
1971 ! Locals:
1972 integer :: nulnam, ierr, fnom, fclos
1973 logical, save :: firstCall = .true.
1974
1975 NAMELIST /NAMSTIO/interpToPhysicsGrid
1976
1977 if (firstCall) then
1978
1979 ! set the default values
1980 interpToPhysicsGrid = .false.
1981
1982 ! read the namelist block if it exists
1983 if ( .not. utl_isNamelistPresent('NAMSTIO','./flnml') ) then
1984 if ( mmpi_myid == 0 ) then
1985 write(*,*) 'readNml (gio): namstio is missing in the namelist. The default values will be taken.'
1986 end if
1987 else
1988 ! Read namelist NAMSTIO
1989 nulnam = 0
1990 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
1991 read(nulnam,nml=namstio,iostat=ierr)
1992 if (ierr /= 0) call utl_abort('readNml (gio): Error reading namelist')
1993 if (mmpi_myid == 0) write(*,nml=namstio)
1994 ierr = fclos(nulnam)
1995 end if
1996
1997 firstCall = .false.
1998
1999 end if
2000
2001 end subroutine readNml
2002
2003end module gridStateVectorFileIO_mod
