1module increment_mod
2 ! MODULE increment_mod (prefix='inc' category='1. High-level functionality')
3 !
4 !:Purpose: To add a 4D increment to a given 4D background/reference state and
5 ! to output the results
6 !
7 use codePrecision_mod
8 use midasMpi_mod
9 use timeCoord_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use interpolation_mod
13 use horizontalCoord_mod
14 use verticalCoord_mod
15 use humidityLimits_mod
16 use utilities_mod
17 use message_mod
18 use gridVariableTransforms_mod
19 use BMatrix_mod
20 use varNamelist_mod
21 implicit none
22 save
23 private
24
25 ! public procedures
26 public :: inc_computeHighResAnalysis, inc_writeIncAndAnalHighRes, inc_getIncrement, inc_writeIncrement
27 public :: inc_writeAnalysis, inc_analPostProcessing
28
29 ! namelist variables
30 integer :: writeNumBits ! number of bits to use when writing analysis and high-res increment
31 logical :: writeHiresIncrement ! choose to write the high-res increment to a file
32 logical :: imposeRttovHuLimits ! choose to impose "rttov" humidity limits to analysis
33 logical :: useAnalIncMask ! for LAM only, choose to apply scale factor from a mask file to the increment
34 character(len=12) :: etiket_anlm ! 'etiket' used when writing the analysis
35 character(len=12) :: etiket_rehm ! 'etiket' used when writing the high-res increment
36 character(len=12) :: etiket_rebm ! 'etiket' used when writing the low-res increment
37 character(len=12) :: hInterpolationDegree ! type of interpolation to use: 'LINEAR' or 'CUBIC'
38 logical :: applyLiebmann ! choose to apply Liebmann extrapolation to SST and/or sea ice
39 logical :: SSTSpread ! choose to apply spatial spreading of the SST increment onto land
40 integer :: SSTSpreadMaxBoxSize ! control the amount of SST increment spreading
41 character(len=10) :: SSTSubgrid ! select subgrid on which to apply spreading: 'Yin' or 'Yan'
42
43CONTAINS
44
45 !--------------------------------------------------------------------------
46 ! readNameList
47 !--------------------------------------------------------------------------
48 subroutine readNameList
49 !
50 !:Purpose: Reading NAMINC namelist by any subroutines in increment_mod module.
51 !
52 implicit none
53
54 ! Locals:
55 integer :: nulnam, ierr
56 integer, external :: fnom, fclos
57 logical, save :: nmlAlreadyRead = .false.
58
59 NAMELIST /NAMINC/ writeHiresIncrement, etiket_rehm, etiket_anlm, &
60 etiket_rebm, writeNumBits, imposeRttovHuLimits, hInterpolationDegree, &
61 useAnalIncMask, applyLiebmann, SSTSpread, SSTSpreadMaxBoxSize, SSTSubgrid
62
63 if ( .not. nmlAlreadyRead ) then
64 nmlAlreadyRead = .true.
65
66 !- Setting default values
67 writeHiresIncrement = .true.
68 imposeRttovHuLimits = .true.
69 useAnalIncMask = .false.
70 etiket_rehm = 'INCREMENT'
71 etiket_rebm = 'INCREMENT'
72 etiket_anlm = 'ANALYSIS'
73 writeNumBits = 16
74 hInterpolationDegree = 'LINEAR'
75 applyLiebmann = .false.
76 SSTSpread = .false.
77 SSTSpreadMaxBoxSize = 0
78 SSTsubgrid = ' '
79
80 if ( .not. utl_isNamelistPresent('NAMINC','./flnml') ) then
81 call msg('readNameList (inc)', &
82 'NAMINC is missing in the namelist. The default values will be taken.', &
83 mpiAll_opt=.false.)
84 else
85 ! Reading the namelist
86 nulnam = 0
87 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
88 read(nulnam, nml=naminc, iostat=ierr)
89 if ( ierr /= 0) call utl_abort('readNameList: Error reading namelist')
90 ierr = fclos(nulnam)
91 end if
92 if ( mmpi_myid == 0 ) write(*,nml=naminc)
93 end if
94
95 end subroutine readNameList
96
97 !--------------------------------------------------------------------------
98 ! inc_computeHighResAnalysis
99 !--------------------------------------------------------------------------
100 subroutine inc_computeHighResAnalysis( statevectorIncLowRes, & ! IN
101 stateVectorUpdateHighRes, stateVectorPsfcHighRes) ! OUT
102 !
103 !:Purpose: Computing high-resolution analysis on the trial grid.
104 !
105 implicit none
106
107 ! Arguments:
108 type(struct_gsv), intent(in) :: statevectorIncLowRes
109 type(struct_gsv), intent(inout) :: stateVectorUpdateHighRes
110 type(struct_gsv), intent(out) :: stateVectorPsfcHighRes
111
112 ! Locals:
113 type(struct_gsv) :: statevectorRef, statevectorRefPsfc
114 type(struct_gsv) :: statevectorIncRefLowRes
115 type(struct_gsv) :: statevectorTrlRefVars
116 type(struct_gsv) :: statevectorTrlLowResTime
117 type(struct_gsv) :: statevectorTrlLowResVert
118 type(struct_gsv) :: statevector_mask
119 type(struct_gsv) :: stateVectorHighRes
120
121 type(struct_vco), pointer :: vco_trl, vco_inc
122 type(struct_hco), pointer :: hco_trl, hco_inc
123
124 real(pre_incrReal), pointer :: PsfcIncPrior(:,:,:,:), PsfcIncMasked(:,:,:,:)
125 real(pre_incrReal), pointer :: analIncMask(:,:,:)
126
127 integer :: stepIndex, numStep_inc, numStep_trl
128 integer, allocatable :: dateStampList(:)
129
130 character(len=4), allocatable :: varNamesRef(:), varNamesPsfc(:)
131 logical :: allocHeightSfc, refStateNeeded
132
133 call msg('inc_computeHighResAnalysis', 'START', verb_opt=2)
134 call msg_memUsage('inc_computeHighResAnalysis')
135
136 call utl_tmg_start(80,'--Increment')
137 call utl_tmg_start(81,'----ComputeHighResAnalysis')
138
139 ! Set/Read values for the namelist NAMINC
140 call readNameList
141
142 if ( gsv_isAllocated(stateVectorPsfcHighRes) ) then
143 call utl_abort('inc_computeHighResAnalysis: '&
144 //'stateVectorPsfcHighRes should not be allocated yet')
145 end if
146
147 ! If P0 present it is infered the statevector is a 3D atmospheric state
148 refStateNeeded = gsv_varExist(varName='P0')
149
150 ! Setup timeCoord module (date read from trial file)
151 numStep_inc = tim_nstepobsinc ! low-res time
152 numStep_trl = tim_nstepobs ! high-res time
153 allocate(dateStampList(numStep_inc))
154 call tim_getstamplist(dateStampList,numStep_inc,tim_getDatestamp())
155
156 hco_trl => gsv_getHco(stateVectorUpdateHighRes)
157 vco_trl => gsv_getVco(stateVectorUpdateHighRes)
158 hco_inc => gsv_getHco(statevectorIncLowRes)
159 vco_inc => gsv_getVco(statevectorIncLowRes)
160
161 if (vco_trl%Vcode == 0 .or. .not. gsv_varExist(varName='P0')) then
162 allocHeightSfc = .false.
163 else
164 allocHeightSfc = stateVectorUpdateHighRes%heightSfcPresent
165 end if
166
167 ! Read the analysis mask (in LAM mode only) - N.B. different from land/sea mask!!!
168 if (.not. hco_trl%global .and. useAnalIncMask) then
169 call gio_getMaskLAM(statevector_mask, hco_trl, vco_trl, hInterpolationDegree)
170 end if
171
172 ref_building: if ( refStateNeeded ) then
173 ! Build a reference variables
174 !
175 ! - Allocate the target reference statevector for vertical interpolation
176 ! that will be done in inc_interpolateAndAdd.
177 ! The reference needs to have the vertical structure of the input which is the
178 ! increment, but since horizontal interpolation is done first, it needs the
179 ! trial horizontal structure.
180 if (vco_trl%Vcode == 5100) then
181 allocate(varNamesRef(4))
182 varNamesRef = (/'P0', 'P0LS', 'TT', 'HU'/)
183 allocate(varNamesPsfc(2))
184 varNamesPsfc = (/'P0','P0LS'/)
185 else
186 allocate(varNamesRef(3))
187 varNamesRef = (/'P0', 'TT', 'HU'/)
188 allocate(varNamesPsfc(1))
189 varNamesPsfc = (/'P0'/)
190 end if
191 call gsv_allocate( statevectorRef, numStep_inc, hco_trl, vco_inc, &
192 dataKind_opt=pre_incrReal, &
193 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
194 varNames_opt=varNamesRef, allocHeightSfc_opt=allocHeightSfc, &
195 hInterpolateDegree_opt=hInterpolationDegree )
196
197 ! - Restriction of increment to only reference state variables
198 call gsv_allocate( statevectorIncRefLowRes, numStep_inc, hco_inc, vco_inc, &
199 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
200 dataKind_opt=pre_incrReal, varNames_opt=varNamesRef, &
201 allocHeightSfc_opt=allocHeightSfc)
202 call gsv_copy(statevectorIncLowRes, statevectorIncRefLowRes, &
203 allowVarMismatch_opt=.true.)
204
205 ! - Horizontal *only* interpolation of analysis increment (ref state variables)
206 ! using int_interp_gsv since it also does needed mpi transpositions
207 call msg('inc_computeHighResAnalysis', &
208 'horizontal interpolation of the Psfc increment', mpiAll_opt=.false.)
209 call int_interp_gsv(statevectorIncRefLowRes,statevectorRef)
210 call gsv_deallocate(statevectorIncRefLowRes)
211 !- Apply mask to increment in LAM
212 if (.not. hco_trl%global .and. useAnalIncMask) then
213 call gsv_getField(statevectorRef,PsfcIncPrior,'P0')
214 call gsv_getField(statevectorRef,PsfcIncMasked,'P0')
215 call gsv_getField(statevector_mask,analIncMask)
216 do stepIndex = 1, statevectorRef%numStep
217 PsfcIncMasked(:,:,1,stepIndex) = PsfcIncPrior(:,:,1,stepIndex) * analIncMask(:,:,1)
218 end do
219 if (gsv_varExist(statevectorRef,'P0LS')) then
220 call gsv_getField(statevectorRef,PsfcIncPrior,'P0LS')
221 call gsv_getField(statevectorRef,PsfcIncMasked,'P0LS')
222 do stepIndex = 1, statevectorRef%numStep
223 PsfcIncMasked(:,:,1,stepIndex) = PsfcIncPrior(:,:,1,stepIndex) * analIncMask(:,:,1)
224 end do
225 end if
226 end if
227
228 ! - Compute analysis of reference state variables to use for vertical interpolation
229 ! of increment
230 call msg('inc_computeHighResAnalysis', &
231 'Computing reference variables analysis to use for interpolation of increment', &
232 mpiAll_opt=.false.)
233
234 ! - build a restriction to reference state variables of trial
235 call gsv_allocate(statevectorTrlRefVars, numStep_trl, &
236 hco_trl, vco_trl, dateStamp_opt=tim_getDateStamp(), &
237 mpi_local_opt=.true., dataKind_opt=pre_incrReal, &
238 varNames_opt=varNamesRef, allocHeightSfc_opt=allocHeightSfc )
239 call gsv_copy(stateVectorUpdateHighRes, statevectorTrlRefVars,&
240 allowVarMismatch_opt=.true.)
241 ! - bring trial to low time res
242 call gsv_allocate(statevectorTrlLowResTime, numStep_inc, hco_trl, vco_trl, &
243 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
244 dataKind_opt=pre_incrReal, varNames_opt=varNamesRef, &
245 allocHeightSfc_opt=allocHeightSfc)
246 call gsv_copy(statevectorTrlRefVars, statevectorTrlLowResTime, &
247 allowTimeMismatch_opt=.true.)
248 ! - bring trial to low res (increment) vertical structure
249 ! (no vertical interpolation reference needed as it is a full state,
250 ! not a incremental state)
251 call gsv_allocate(statevectorTrlLowResVert, numStep_inc, hco_trl, vco_inc, &
252 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
253 dataKind_opt=pre_incrReal, varNames_opt=varNamesRef, &
254 allocHeightSfc_opt=allocHeightSfc)
255 call int_vInterp_gsv(statevectorTrlLowResTime, statevectorTrlLowResVert)
256 call gsv_deallocate(statevectorTrlLowResTime)
257 call gsv_deallocate(statevectorTrlRefVars)
258 ! - build the analysis reference state for the increment vertical interpolation.
259 ! At that stage, statevectorRef contains the increment on the trial
260 ! horizontal grid, but analysis vertical structure; consistent for addition.
261 call gsv_add(statevectorTrlLowResVert, statevectorRef)
262 call gsv_deallocate(statevectorTrlLowResVert)
263
264 ! - Time interpolation to get high-res Psfc analysis increment to output
265 call msg('inc_computeHighResAnalysis', &
266 'Time interpolation to get high-res Psfc analysis increment', &
267 mpiAll_opt=.false.)
268
269 ! - Restriction to P0 only + vco set on vco_trl for later consistency
270 call gsv_allocate(statevectorRefPsfc, numStep_inc, hco_trl, vco_trl, &
271 dataKind_opt=pre_incrReal, &
272 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
273 varNames_opt=varNamesPsfc, allocHeightSfc_opt=allocHeightSfc, &
274 hInterpolateDegree_opt=hInterpolationDegree )
275 call gsv_copy(statevectorRef, statevectorRefPsfc, allowVarMismatch_opt=.true., &
276 allowVcoMismatch_opt=.true.)
277 ! - high res time interpolation
278 call gsv_allocate(stateVectorPsfcHighRes, numStep_trl, hco_trl, vco_trl, &
279 dataKind_opt=pre_incrReal, &
280 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
281 varNames_opt=varNamesPsfc, allocHeightSfc_opt=allocHeightSfc, &
282 hInterpolateDegree_opt=hInterpolationDegree)
283 call int_tInterp_gsv(statevectorRefPsfc, stateVectorPsfcHighRes)
284 else ref_building
285 call msg('inc_computeHighResAnalysis','P0 not present, not computing reference', &
286 verb_opt=3)
287 end if ref_building
288
289 ! Compute the analysis
290 call msg('inc_computeHighResAnalysis', 'compute the analysis', mpiAll_opt=.false.)
291
292 ! Copy of trial to be passed to inc_interpolateAndAdd
293 call gsv_allocate( stateVectorHighRes, numStep_trl, hco_trl, vco_trl, &
294 dataKind_opt=stateVectorUpdateHighRes%dataKind, &
295 dateStamp_opt=tim_getDateStamp(), &
296 mpi_local_opt=stateVectorUpdateHighRes%mpi_local, &
297 allocHeightSfc_opt=stateVectorUpdateHighRes%heightSfcPresent, &
298 hInterpolateDegree_opt=hInterpolationDegree, &
299 allocHeight_opt=.false., allocPressure_opt=.false. )
300 call gsv_copy( stateVectorUpdateHighRes,stateVectorHighRes, &
301 allowVarMismatch_opt=.true.)
302
303 ! Interpolate low-res increments to high-res and add to the initial state
304 if (.not. hco_trl%global .and. useAnalIncMask) then
305 if ( refStateNeeded ) then
306 call inc_interpolateAndAdd(statevectorIncLowRes, stateVectorHighRes, &
307 statevectorRef_opt=statevectorRef, &
308 statevectorMaskLAM_opt=statevector_mask)
309 else
310 call inc_interpolateAndAdd(statevectorIncLowRes, stateVectorHighRes, &
311 statevectorMaskLAM_opt=statevector_mask)
312 end if
313 else
314 if ( refStateNeeded ) then
315 call inc_interpolateAndAdd(statevectorIncLowRes, stateVectorHighRes, &
316 statevectorRef_opt=statevectorRef)
317 else
318 call inc_interpolateAndAdd(statevectorIncLowRes, stateVectorHighRes)
319 end if
320 end if
321
322 call gsv_copy( stateVectorHighRes, stateVectorUpdateHighRes, &
323 allowVarMismatch_opt=.true.)
324 call gsv_deallocate(stateVectorHighRes)
325
326 if ( gsv_isAllocated(statevectorRef) ) call gsv_deallocate(statevectorRef)
327 if ( gsv_isAllocated(statevector_mask) ) call gsv_deallocate(statevector_mask)
328
329 call utl_tmg_stop(81)
330 call utl_tmg_stop(80)
331
332 call msg('inc_computeHighResAnalysis', 'END', verb_opt=2)
333
334 end subroutine inc_computeHighResAnalysis
335
336 !--------------------------------------------------------------------------
337 ! inc_analPostProcessing
338 !--------------------------------------------------------------------------
339 subroutine inc_analPostProcessing (stateVectorPsfcHighRes, stateVectorUpdateHighRes, & ! IN
340 stateVectorTrial, stateVectorPsfc, stateVectorAnal) ! OUT
341 !
342 !:Purpose: Post processing of the high resolution analysis including degrading
343 ! temporal resolution, variable transform, and humidity clipping.
344 !
345 implicit none
346
347 ! Arguments:
348 type(struct_gsv), intent(inout) :: stateVectorPsfcHighRes
349 type(struct_gsv), intent(in) :: stateVectorUpdateHighRes
350 type(struct_gsv), intent(out) :: stateVectorTrial
351 type(struct_gsv), intent(out) :: stateVectorPsfc
352 type(struct_gsv), intent(out) :: stateVectorAnal
353
354 ! Locals:
355 type(struct_vco), pointer :: vco_trl => null()
356 type(struct_hco), pointer :: hco_trl => null()
357 character(len=4), allocatable :: varNamesPsfc(:)
358
359 real(pre_incrReal), pointer :: oceanIce_ptr(:,:,:,:)
360
361 logical :: allocHeightSfc
362
363 call msg('inc_analPostProcessing', 'START', verb_opt=2)
364 call msg_memUsage('inc_analPostProcessing')
365
366 call utl_tmg_start(80,'--Increment')
367 call utl_tmg_start(82,'----AnalPostProcessing')
368
369 ! Re-read trials to make stateVectorTrial with degraded timesteps available
370 hco_trl => gsv_getHco(stateVectorUpdateHighRes)
371 vco_trl => gsv_getVco(stateVectorUpdateHighRes)
372 allocHeightSfc = ( vco_trl%Vcode /= 0 )
373
374 call gsv_allocate(stateVectorTrial, tim_nstepobsinc, hco_trl, vco_trl, &
375 dataKind_opt = pre_incrReal, &
376 dateStamp_opt = tim_getDateStamp(), mpi_local_opt = .true., &
377 allocHeightSfc_opt = allocHeightSfc, hInterpolateDegree_opt = 'LINEAR', &
378 allocHeight_opt = .false., allocPressure_opt = .false.)
379 call gsv_zero(stateVectorTrial)
380 call gio_readTrials(stateVectorTrial)
381 call msg_memUsage('inc_analPostProcessing')
382
383 if (vco_trl%Vcode == 0 .or. .not. gsv_varExist(varName='P0')) then
384 allocHeightSfc = .false.
385 else
386 allocHeightSfc = stateVectorTrial%heightSfcPresent
387 end if
388
389 if (vco_trl%Vcode == 5100) then
390 allocate(varNamesPsfc(2))
391 varNamesPsfc = (/'P0','P0LS'/)
392 else
393 allocate(varNamesPsfc(1))
394 varNamesPsfc = (/'P0'/)
395 end if
396
397 ! Degrade timesteps stateVector high-res Psfc, and high-res analysis
398 if (gsv_varExist(varName='P0')) then
399 call gsv_allocate(stateVectorPsfc, tim_nstepobsinc, hco_trl, vco_trl, &
400 dataKind_opt = pre_incrReal, &
401 dateStamp_opt = tim_getDateStamp(), mpi_local_opt = .true., &
402 varNames_opt = varNamesPsfc, allocHeightSfc_opt = allocHeightSfc, &
403 hInterpolateDegree_opt = hInterpolationDegree)
404 call gsv_copy(stateVectorPsfcHighRes, stateVectorPsfc, &
405 allowTimeMismatch_opt = .true., allowVarMismatch_opt=.true.)
406 call gsv_deallocate(stateVectorPsfcHighRes)
407 end if
408
409 call gsv_allocate(stateVectorAnal, tim_nstepobsinc, hco_trl, vco_trl, &
410 dataKind_opt = pre_incrReal, &
411 dateStamp_opt = tim_getDateStamp(), mpi_local_opt = .true., &
412 allocHeightSfc_opt = allocHeightSfc, hInterpolateDegree_opt = 'LINEAR', &
413 allocHeight_opt = .false., allocPressure_opt = .false.)
414 call gsv_copy(stateVectorUpdateHighRes, stateVectorAnal, &
415 allowVarMismatch_opt = .true., allowTimeMismatch_opt = .true.)
416
417 if (gsv_varExist(stateVectorAnal, 'GL')) then
418 ! Impose limits [0,1] on sea ice concentration analysis
419 call gsv_getField(stateVectorAnal, oceanIce_ptr, 'GL')
420 oceanIce_ptr(:,:,:,:) = min(oceanIce_ptr(:,:,:,:), 1.0d0)
421 oceanIce_ptr(:,:,:,:) = max(oceanIce_ptr(:,:,:,:), 0.0d0)
422 end if
423
424 if (applyLiebmann) then
425
426 ! Start the variable transformations
427 if (gsv_varExist(stateVectorAnal, 'GL')) then
428 if (gsv_varExist(stateVectorAnal, 'LG')) then
429 ! Compute the continuous sea ice concentration field (LG)
430 call gvt_transform(stateVectorAnal, 'oceanIceContinuous', stateVectorRef_opt = stateVectorTrial, varName_opt = 'LG')
431 end if
432 end if
433
434 ! Start the variable transformations
435 if(gsv_varExist(stateVectorAnal,'TM')) then
436 ! Compute the continuous SST field (TM)
437 call gvt_transform(stateVectorAnal, 'oceanIceContinuous', stateVectorRef_opt = stateVectorTrial, varName_opt = 'TM')
438 end if
439
440 end if
441
442 if (SSTSpread) then
443 if(gsv_varExist(stateVectorAnal,'TM')) then
444 ! Compute SST spreading (TM) on neigbouring land points
445 call gvt_transform(stateVectorAnal, 'SSTSpread', varName_opt = 'TM', maxBoxSize_opt = SSTSpreadMaxBoxSize, subgrid_opt = SSTSubgrid)
446 end if
447 end if
448
449 ! Convert all transformed variables into model variables (e.g. LVIS->VIS, LPR->PR) for analysis
450 call gvt_transform(stateVectorAnal, 'AllTransformedToModel', allowOverWrite_opt = .true.)
451
452 ! Impose limits on humidity analysis
453 call msg('inc_analPostProcessing', 'calling qlim_saturationLimit')
454 call qlim_saturationLimit(stateVectorAnal)
455 if (imposeRttovHuLimits) call qlim_rttovLimit(stateVectorAnal)
456
457 if (gsv_varKindExist('CH')) then
458 ! Apply boundaries to analysis of CH kind variables as needed.
459 call msg('inc_analPostProcessing', &
460 'applying minimum values to analysis for variables of CH kind')
461 call gvt_transform(stateVectorAnal, 'CH_bounds')
462 end if
463
464 call utl_tmg_stop(82)
465 call utl_tmg_stop(80)
466
467 call msg('inc_analPostProcessing', 'END', verb_opt=2)
468
469 end subroutine inc_analPostProcessing
470
471 !--------------------------------------------------------------------------
472 ! inc_writeIncAndAnalHighRes
473 !--------------------------------------------------------------------------
474 subroutine inc_writeIncAndAnalHighRes( stateVectorTrial, stateVectorPsfc, &
475 stateVectorAnal )
476 !
477 !:Purpose: Write the high-resolution analysis increments to the rehm file.
478 !
479 implicit none
480
481 ! Arguments:
482 type(struct_gsv), target, intent(inout) :: stateVectorTrial
483 type(struct_gsv), intent(inout) :: stateVectorPsfc
484 type(struct_gsv), intent(inout) :: stateVectorAnal
485
486 ! Locals:
487 type(struct_gsv) :: stateVectorIncHighRes
488 type(struct_gsv) :: stateVector_1step_r4, stateVectorPsfc_1step_r4
489 character(len=4), allocatable :: varNamesPsfc(:)
490
491 type(struct_vco), pointer :: vco_trl => null()
492 type(struct_hco), pointer :: hco_trl => null()
493
494 integer :: stepIndex, stepIndexBeg, stepIndexEnd, stepIndexToWrite, numStep
495 integer :: dateStamp, numBatch, batchIndex, procToWrite
496 integer, allocatable :: dateStampList(:)
497
498 character(len=256) :: incFileName, anlFileName
499 character(len=4) :: coffset
500 character(len=4), pointer :: varNames(:)
501
502 real(8) :: deltaHours
503 real(8), pointer :: HeightSfc_increment(:,:), HeightSfc_trial(:,:)
504
505 logical :: allocHeightSfc, writeHeightSfc
506
507 call msg('inc_writeIncAndAnalHighRes', 'START', verb_opt=2)
508 call msg_memUsage('inc_writeIncAndAnalHighRes')
509
510 call utl_tmg_start(80,'--Increment')
511 call utl_tmg_start(83,'----WriteIncAndAnalHighRes')
512
513 ! Setup timeCoord module (date read from trial file)
514 numStep = tim_nstepobsinc
515 allocate(dateStampList(numStep))
516 call tim_getstamplist(dateStampList,numStep,tim_getDatestamp())
517
518 hco_trl => gsv_getHco(stateVectorTrial)
519 vco_trl => gsv_getVco(stateVectorTrial)
520 if (vco_trl%Vcode == 0 .or. .not. gsv_varExist(varName='P0')) then
521 allocHeightSfc = .false.
522 else
523 allocHeightSfc = stateVectorTrial%heightSfcPresent
524 end if
525 writeHeightSfc = allocHeightSfc
526
527 if (vco_trl%Vcode == 5100) then
528 allocate(varNamesPsfc(2))
529 varNamesPsfc = (/'P0','P0LS'/)
530 else
531 allocate(varNamesPsfc(1))
532 varNamesPsfc = (/'P0'/)
533 end if
534
535 ! Convert all transformed variables into model variables (e.g. LVIS->VIS, LPR->PR) for original trial
536 call gvt_transform(stateVectorTrial, 'AllTransformedToModel',allowOverWrite_opt=.true.)
537
538 ! Copy the surface height from trial into stateVectorPsfc and stateVectorAnal
539 if ( allocHeightSfc ) then
540 HeightSfc_trial => gsv_getHeightSfc(stateVectorTrial)
541 HeightSfc_increment => gsv_getHeightSfc(stateVectorPsfc)
542 HeightSfc_increment(:,:) = HeightSfc_trial(:,:)
543
544 nullify(HeightSfc_increment)
545 HeightSfc_increment => gsv_getHeightSfc(stateVectorAnal)
546 HeightSfc_increment(:,:) = HeightSfc_trial(:,:)
547 end if
548
549 ! reAllocate incHighRes with the names of the model variables (e.g. VIS, PR)
550 nullify(varNames)
551 call gsv_varNamesList(varNames, stateVectorAnal)
552 call gsv_allocate(stateVectorIncHighRes, numStep, hco_trl, vco_trl, &
553 dataKind_opt=pre_incrReal, &
554 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
555 allocHeightSfc_opt=allocHeightSfc, hInterpolateDegree_opt=hInterpolationDegree, &
556 varNames_opt=varNames )
557
558 ! Recompute increments and write out the files in parallel with respect to time steps
559 call gsv_copy(stateVectorAnal, stateVectorIncHighRes)
560 call gsv_add(stateVectorTrial, stateVectorIncHighRes, -1.0d0)
561
562 ! figure out number of batches of time steps for writing
563 numBatch = ceiling(real(stateVectorAnal%numStep) / real(mmpi_nprocs))
564 call msg('inc_writeIncAndAnalHighRes', &
565 'writing will be done by number of batches = '//str(numBatch))
566
567 batch_loop: do batchIndex = 1, numBatch
568
569 stepIndexBeg = 1 + (batchIndex - 1) * mmpi_nprocs
570 stepIndexEnd = min(stateVectorAnal%numStep, stepIndexBeg + mmpi_nprocs - 1)
571 call msg('inc_writeIncAndAnalHighRes', 'batchIndex = '//str(batchIndex) &
572 //', stepIndexBeg = '//str(stepIndexBeg) &
573 //', stepIndexEnd = '//str(stepIndexEnd))
574
575 ! figure out which time step I will write, if any (-1 if none)
576 stepIndexToWrite = -1
577 do stepIndex = stepIndexBeg, stepIndexEnd
578 procToWrite = nint( real(stepIndex - stepIndexBeg) * real(mmpi_nprocs) / real(stepIndexEnd - stepIndexBeg + 1) )
579 if ( procToWrite == mmpi_myid ) stepIndexToWrite = stepIndex
580 call msg('inc_writeIncAndAnalHighRes', ' stepIndex = '//str(stepIndex) &
581 //', procToWrite = '//str(procToWrite), mpiAll_opt=.false.)
582 end do
583
584 ! determine date and allocate stateVector for storing just 1 time step, if I do writing
585 if ( stepIndexToWrite /= -1 ) then
586
587 dateStamp = stateVectorAnal%dateStampList(stepIndexToWrite)
588 call difdatr(dateStamp,tim_getDatestamp(),deltaHours)
589 if(nint(deltaHours*60.0d0).lt.0) then
590 write(coffset,'(I4.3)') nint(deltaHours*60.0d0)
591 else
592 write(coffset,'(I3.3)') nint(deltaHours*60.0d0)
593 endif
594
595 call gsv_allocate( stateVector_1step_r4, 1, stateVectorAnal%hco, stateVectorAnal%vco, &
596 dateStamp_opt=dateStamp, mpi_local_opt=.false., dataKind_opt=4, &
597 allocHeightSfc_opt=allocHeightSfc, &
598 varNames_opt=varNames )
599 call gsv_allocate( stateVectorPsfc_1step_r4, 1, stateVectorAnal%hco, stateVectorAnal%vco, &
600 dateStamp_opt=dateStamp, mpi_local_opt=.false., dataKind_opt=4, &
601 varNames_opt=varNamesPsfc, allocHeightSfc_opt=allocHeightSfc )
602 call gsv_zero(stateVectorPsfc_1step_r4)
603 end if
604
605 ! transpose ANALYSIS data from Tiles to Steps
606 call gsv_transposeTilesToStep(stateVector_1step_r4, stateVectorAnal, stepIndexBeg)
607
608 ! write the ANALYSIS file for one timestep on all tasks with data
609 if ( stepIndexToWrite /= -1 ) then
610 call msg_memUsage('inc_writeIncAndAnalHighRes')
611 anlFileName = './anlm_' // trim(coffset) // 'm'
612 call gio_writeToFile( stateVector_1step_r4, trim(anlFileName), etiket_anlm, &
613 typvar_opt='A', writeHeightSfc_opt=writeHeightSfc, &
614 numBits_opt=writeNumBits, containsFullField_opt=.true. )
615 end if
616
617 if (writeHiresIncrement) then
618
619 ! transpose INCREMENT data from Tiles to Steps
620 call gsv_transposeTilesToStep(stateVector_1step_r4, stateVectorIncHighRes, stepIndexBeg)
621
622 ! write the INCREMENT file for one timestep on all tasks with data
623 if ( stepIndexToWrite /= -1 ) then
624 call msg_memUsage('inc_writeIncAndAnalHighRes')
625 incFileName = './rehm_' // trim(coffset) // 'm'
626 call gio_writeToFile( stateVector_1step_r4, trim(incFileName), etiket_rehm, &
627 typvar_opt='R', &
628 numBits_opt=writeNumBits, containsFullField_opt=.false. )
629 end if
630
631 if (gsv_varExist(varName='P0')) then
632
633 ! transpose ANALYSIS PSFC AND height_SFC ONLY data from Tiles to Steps
634 call gsv_transposeTilesToStep(stateVectorPsfc_1step_r4, stateVectorPsfc, stepIndexBeg)
635
636 ! Also write analysis value of Psfc and surface height to increment file
637 if ( stepIndexToWrite /= -1 ) then
638 call msg_memUsage('inc_writeIncAndAnalHighRes')
639 incFileName = './rehm_' // trim(coffset) // 'm'
640 call gio_writeToFile( stateVectorPsfc_1step_r4, trim(incFileName), etiket_rehm, &
641 typvar_opt='A', writeHeightSfc_opt=writeHeightSfc, &
642 numBits_opt=writeNumBits, containsFullField_opt=.true. )
643 end if
644
645 end if
646
647 end if ! writeHiresIncrement
648
649 if ( stepIndexToWrite /= -1 ) then
650 call gsv_deallocate(stateVector_1step_r4)
651 call gsv_deallocate(stateVectorPsfc_1step_r4)
652 end if
653
654 end do batch_loop
655
656 call gsv_deallocate(stateVectorAnal)
657 if (gsv_varExist(varName='P0')) then
658 call gsv_deallocate(stateVectorPsfc)
659 end if
660 call gsv_deallocate(stateVectorIncHighRes)
661 call gsv_deallocate(stateVectorTrial)
662
663 call utl_tmg_stop(83)
664 call utl_tmg_stop(80)
665
666 call msg('inc_writeIncAndAnalHighRes', 'END', verb_opt=2)
667
668 end subroutine inc_writeIncAndAnalHighRes
669
670 !--------------------------------------------------------------------------
671 ! inc_getIncrement
672 !--------------------------------------------------------------------------
673 subroutine inc_getIncrement(incr_cv,statevector_incr,nvadim_mpilocal)
674 !
675 !:Purpose: Get true analysis increment from control vector.
676 !
677 implicit none
678
679 ! Arguments:
680 real(8), intent(in) :: incr_cv(:)
681 type(struct_gsv), intent(inout) :: statevector_incr
682 integer, intent(in) :: nvadim_mpilocal
683
684 call utl_tmg_start(80,'--Increment')
685 call utl_tmg_start(84,'----GetIncrement')
686
687 ! compute increment from control vector (multiply by B^1/2)
688 call bmat_sqrtB( incr_cv, nvadim_mpilocal, statevector_incr )
689
690 ! Compute new diagnotics based on NAMSTATE
691 if ( gsv_varExist(statevector_incr,'QR') .and. gsv_varExist(statevector_incr,'DD') ) then
692 call msg('inc_getIncrement', 'User is asking for Vort-Div analysis increment')
693 call gvt_transform( statevector_incr, & ! INOUT
694 'UVtoVortDiv' ) ! IN
695 if ( gsv_varExist(statevector_incr,'PP') .and. gsv_varExist(statevector_incr,'CC') ) then
696 call msg('inc_getIncrement', 'User is asking for Psi-Chi analysis increment')
697 call gvt_transform( statevector_incr, & ! INOUT
698 'VortDivToPsiChi') ! IN
699 end if
700 end if
701
702 call utl_tmg_stop(84)
703 call utl_tmg_stop(80)
704
705 end subroutine inc_getIncrement
706
707 !--------------------------------------------------------------------------
708 ! inc_writeIncrement
709 !--------------------------------------------------------------------------
710 subroutine inc_writeIncrement( stateVector_incr, &
711 ip3ForWriteToFile_opt )
712 !
713 !:Purpose: Write the low-resolution analysis increments to the rebm file.
714 !
715 implicit none
716
717 ! Arguments:
718 type(struct_gsv), intent(in) :: stateVector_incr
719 integer, optional, intent(in) :: ip3ForWriteToFile_opt
720
721 ! Locals:
722 integer :: stepIndex, dateStamp
723 real(8) :: deltaHours
724 character(len=4) :: coffset
725 character(len=30) :: fileName
726
727 call msg('inc_writeIncrement', 'START', verb_opt=2)
728
729 call utl_tmg_start(80,'--Increment')
730 call utl_tmg_start(85,'----WriteIncrement')
731
732 ! loop over times for which increment is computed
733 do stepIndex = 1, tim_nstepobsinc
734 if (gsv_isAllocated(statevector_incr)) then
735 dateStamp = gsv_getDateStamp(stateVector_incr,stepIndex)
736 call msg('inc_writeIncrement', 'writing increment for time step: ' &
737 //str(stepIndex)//', dateStamp = '//str(dateStamp), mpiAll_opt=.false.)
738
739 ! write the increment file for this time step
740 call difdatr(dateStamp,tim_getDatestamp(),deltaHours)
741 if(nint(deltaHours*60.0d0).lt.0) then
742 write(coffset,'(I4.3)') nint(deltaHours*60.0d0)
743 else
744 write(coffset,'(I3.3)') nint(deltaHours*60.0d0)
745 end if
746 fileName = './rebm_' // trim(coffset) // 'm'
747 call gio_writeToFile( stateVector_incr, fileName, etiket_rebm, scaleFactor_opt=1.0d0, &
748 ip3_opt=ip3ForWriteToFile_opt, stepIndex_opt=stepIndex, &
749 containsFullField_opt=.false. )
750 end if
751 end do
752
753 call utl_tmg_stop(85)
754 call utl_tmg_stop(80)
755
756 call msg('inc_writeIncrement', 'END', verb_opt=2)
757
758
759 end subroutine inc_writeIncrement
760
761 !--------------------------------------------------------------------------
762 ! inc_writeAnalysis
763 !--------------------------------------------------------------------------
764 subroutine inc_writeAnalysis(statevector_anal)
765 !:Purpose: To write to output standard file the analysid from statevector strucure (1Dvar case)
766 ! to output the results
767 !
768 implicit none
769
770 ! Arguments:
771 type(struct_gsv), intent(in) :: statevector_anal
772
773 ! Locals:
774 integer :: stepIndex, dateStamp
775 real(8) :: deltaHours
776 character(len=4) :: coffset
777 character(len=30) :: fileName
778
779 call msg('inc_writeAnalysis', 'START', verb_opt=2)
780
781 call utl_tmg_start(80,'--Increment')
782 call utl_tmg_start(86,'----WriteAnalysis')
783
784 !
785 !- Set/Read values for the namelist NAMINC
786 !
787 call readNameList()
788
789 ! loop over times for which increment is computed
790 do stepIndex = 1, tim_nstepobsinc
791 if (gsv_isAllocated(statevector_anal)) then
792 dateStamp = gsv_getDateStamp(statevector_anal,stepIndex)
793 call msg('inc_writeAnalysis', 'writing analysis for time step: ' &
794 //str(stepIndex)//', dateStamp = '//str(dateStamp), mpiAll_opt=.false.)
795
796 ! write the increment file for this time step
797 call difdatr(dateStamp,tim_getDatestamp(),deltaHours)
798 if(nint(deltaHours*60.0d0).lt.0) then
799 write(coffset,'(I4.3)') nint(deltaHours*60.0d0)
800 else
801 write(coffset,'(I3.3)') nint(deltaHours*60.0d0)
802 end if
803 fileName = './anlm_' // trim(coffset) // 'm'
804 call gio_writeToFile( statevector_anal, fileName, etiket_anlm, scaleFactor_opt = 1.0d0, &
805 ip3_opt = 0, stepIndex_opt = stepIndex, containsFullField_opt=.true. )
806 end if
807 end do
808
809 call utl_tmg_stop(86)
810 call utl_tmg_stop(80)
811
812 call msg('inc_writeAnalysis', 'END', verb_opt=2)
813 end subroutine inc_writeAnalysis
814
815 !--------------------------------------------------------------------------
816 ! inc_interpolateAndAdd
817 !--------------------------------------------------------------------------
818 subroutine inc_interpolateAndAdd(statevector_in,statevector_inout, &
819 statevectorRef_opt, statevectorMaskLAM_opt, &
820 scaleFactor_opt)
821 !
822 !:Purpose: Interpolate the low-resolution increments to trial grid and add to
823 ! get the high-resolution analysis.
824 !
825 implicit none
826
827 ! Arguments:
828 type(struct_gsv), target, intent(in) :: statevector_in
829 type(struct_gsv), intent(inout) :: statevector_inout
830 type(struct_gsv), optional, intent(in) :: statevectorRef_opt ! Reference statevector providing optional fields (P0, TT, HU)
831 type(struct_gsv), optional, intent(in) :: statevectorMaskLAM_opt
832 real(8), optional, intent(in) :: scaleFactor_opt
833
834 ! Locals:
835 type(struct_gsv) :: statevector_in_hvInterp
836 type(struct_gsv) :: statevector_in_hvtInterp
837 type(struct_gsv), target :: statevector_in_withP0LS
838 type(struct_gsv), pointer :: statevector_in_ptr
839
840 character(len=4), pointer :: varNamesToInterpolate(:), varNamesToInterpolate_withP0LS(:)
841
842 call msg('inc_interpolateAndAdd', 'START', verb_opt=2)
843
844 ! Error traps
845 if ( .not. gsv_isAllocated(statevector_in) ) then
846 call utl_abort('inc_interpolateAndAdd: gridStateVector_in not yet allocated! Aborting.')
847 end if
848 if ( .not. gsv_isAllocated(statevector_inout) ) then
849 call utl_abort('inc_interpolateAndAdd: gridStateVector_inout not yet allocated! Aborting.')
850 end if
851 if ( present(statevectorRef_opt) ) then
852 if ( statevector_in%numstep /= statevectorRef_opt%numstep) then
853 call utl_abort('inc_interpolateAndAdd: statevector_in and statevectorRef_opt numstep inconsistent')
854 end if
855 end if
856
857 nullify(varNamesToInterpolate)
858 call vnl_varNamesFromExistList(varNamesToInterpolate, statevector_in%varExistlist(:))
859
860 ! Check if P0LS needs to be added to the increment
861 if (statevector_inout%vco%vcode == 5100 .and. &
862 .not.gsv_varExist(statevector_in, 'P0LS')) then
863 varNamesToInterpolate_withP0LS => vnl_addToVarNames(varNamesToInterpolate,'P0LS')
864 call gsv_allocate(statevector_in_withP0LS, statevector_in%numstep, &
865 statevector_in%hco, statevector_in%vco, &
866 dateStamp_opt=tim_getDateStamp(), &
867 mpi_local_opt=statevector_in%mpi_local, mpi_distribution_opt='Tiles', &
868 dataKind_opt=statevector_in%dataKind, &
869 allocHeightSfc_opt=statevector_in%heightSfcPresent, &
870 varNames_opt=varNamesToInterpolate_withP0LS, &
871 hInterpolateDegree_opt=statevector_in%hInterpolateDegree)
872 call gsv_zero(statevector_in_withP0LS)
873 call gsv_copy(statevector_in, statevector_in_withP0LS, &
874 allowVarMismatch_opt=.true.)
875 statevector_in_ptr => statevector_in_withP0LS
876 deallocate(varNamesToInterpolate)
877 varNamesToInterpolate => varNamesToInterpolate_withP0LS
878 else
879 statevector_in_ptr => statevector_in
880 end if
881 ! Do the spatial interpolation of statevector_in onto the grid of statevector_inout
882 call gsv_allocate(statevector_in_hvInterp, statevector_in%numstep, &
883 statevector_inout%hco, statevector_inout%vco, &
884 dateStamp_opt=tim_getDateStamp(), &
885 mpi_local_opt=statevector_inout%mpi_local, mpi_distribution_opt='Tiles', &
886 dataKind_opt=statevector_inout%dataKind, &
887 allocHeightSfc_opt=statevector_inout%heightSfcPresent, &
888 varNames_opt=varNamesToInterpolate, &
889 hInterpolateDegree_opt=statevector_inout%hInterpolateDegree, &
890 hExtrapolateDegree_opt='VALUE' )
891
892 call int_interp_gsv(statevector_in, statevector_in_hvInterp, &
893 statevectorRef_opt=statevectorRef_opt)
894
895 ! Do the time interpolation
896 call gsv_allocate(statevector_in_hvtInterp, statevector_inout%numstep, &
897 statevector_inout%hco, statevector_inout%vco, &
898 dateStamp_opt=tim_getDateStamp(), &
899 mpi_local_opt=statevector_inout%mpi_local, mpi_distribution_opt='Tiles', &
900 dataKind_opt=statevector_inout%dataKind, &
901 allocHeightSfc_opt=statevector_inout%heightSfcPresent, &
902 varNames_opt=varNamesToInterpolate, &
903 hInterpolateDegree_opt=statevector_inout%hInterpolateDegree, &
904 hExtrapolateDegree_opt='VALUE' )
905 call int_tInterp_gsv(statevector_in_hvInterp, statevector_in_hvtInterp)
906 call gsv_deallocate(statevector_in_hvInterp)
907
908 ! Masking
909 if (present(statevectorMaskLAM_opt)) then
910 call gsv_applyMaskLAM(statevector_in_hvtInterp,statevectorMaskLAM_opt)
911 end if
912
913 ! Do the summation
914 call gsv_add(statevector_in_hvtInterp,statevector_inout,scaleFactor_opt=scaleFactor_opt)
915
916 call gsv_deallocate(statevector_in_hvtInterp)
917 deallocate(varNamesToInterpolate)
918
919 call msg('inc_interpolateAndAdd', 'END', verb_opt=2)
920
921 end subroutine inc_interpolateAndAdd
922
923END MODULE increment_mod
