1module gridVariableTransforms_mod
2 ! MODULE gridVariableTransforms_mod (prefix='gvt' category='4. Data Object transformations')
3 !
4 !:Purpose: To store various functions for variable transforms using inputs
5 ! from gridStateVector(s). Outputs are also placed in a
6 ! gridStateVector.
7 !
8 use midasMpi_mod
9 use codePrecision_mod
10 use mathPhysConstants_mod
11 use earthConstants_mod
12 use timeCoord_mod
13 use gridStateVector_mod
14 use gridStateVectorFileIO_mod
15 use interpolation_mod
16 use ensembleStateVector_mod
17 use lamSpectralTransform_mod
18 use globalSpectralTransform_mod
19 use lamAnalysisGridTransforms_mod
20 use horizontalCoord_mod
21 use verticalCoord_mod
22 use utilities_mod
23 use varNameList_mod
24 use calcHeightAndPressure_mod
25 use utilities_mod
26 use humiditylimits_mod
27 use getGridPosition_mod
28
29 implicit none
30 save
31 private
32
33 ! public procedures
34 public :: gvt_setup, gvt_transform, gvt_getStateVectorTrial
35 public :: gvt_setupRefFromTrialFiles, gvt_setupRefFromStateVector
36
37 logical :: varKindCHTrialsInitialized(vnl_numVarMax) = .false.
38
39 type(struct_hco), pointer :: hco_anl => null()
40 type(struct_vco), pointer :: vco_anl => null()
41 type(struct_hco), pointer :: hco_trl => null()
42 type(struct_vco), pointer :: vco_trl => null()
43
44 type(struct_gsv), target :: stateVectorRefHU
45 type(struct_gsv), target :: stateVectorTrialvarKindCH(vnl_numVarMax)
46 type(struct_gsv), target :: stateVectorRefHeight
47
48 ! module interfaces
49 interface gvt_transform
50 module procedure gvt_transform_gsv
51 module procedure gvt_transform_ens
52 end interface gvt_transform
53
54CONTAINS
55
56 !--------------------------------------------------------------------------
57 ! gvt_setup
58 !--------------------------------------------------------------------------
59 subroutine gvt_setup(hco_in,hco_core,vco_in)
60 !
61 !:Purpose: To set up a variable transformation object
62 !
63 implicit none
64
65 ! Arguments:
66 type(struct_hco), pointer, intent(inout) :: hco_in
67 type(struct_hco), pointer, intent(inout) :: hco_core
68 type(struct_vco), pointer, intent(inout) :: vco_in
69
70 if ( gsv_containsNonZeroValues(stateVectorRefHU) ) return
71 if ( gsv_containsNonZeroValues(stateVectorRefHeight) ) return
72 if ( any(varKindCHTrialsInitialized(:)) ) return
73
74 write(*,*) 'gvt_setup: starting'
75
76 hco_anl => hco_in
77 vco_anl => vco_in
78
79 call lgt_setupFromHco(hco_anl,hco_core)
80
81 write(*,*) 'gvt_setup: done'
82
83 end subroutine gvt_setup
84
85 !--------------------------------------------------------------------------
86 ! gvt_setupRefFromTrialFiles
87 !--------------------------------------------------------------------------
88 subroutine gvt_setupRefFromTrialFiles(varName, varKind_opt)
89 !
90 !:Purpose: Initialise reference statevector from file
91 !
92 !:Arguments:
93 ! :varKind_opt: optional variable "kind" argument presently used to
94 ! initialise the reference state in a chemical assimilation
95 ! context.
96 !
97 implicit none
98
99 ! Arguments:
100 character(len=*), intent(in) :: varName ! reference variable/type used
101 character(len=*), optional, intent(in) :: varKind_opt ! additional variable/type information mandatory for some initialization
102
103 ! Locals:
104 type(struct_gsv) :: statevector_noZnoP
105 integer :: varIndex
106
107 select case ( trim(varName) )
108 case ('HU')
109 ! initialize stateVectorRefHU on analysis grid
110 call gsv_allocate(stateVectorRefHU, tim_nstepobsinc, hco_anl, vco_anl, &
111 dateStamp_opt=tim_getDateStamp(), &
112 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
113 hInterpolateDegree_opt='LINEAR', &
114 varNames_opt=(/'HU','P0'/) )
115
116 ! read trial files using default horizontal interpolation degree
117 call gio_readTrials( stateVectorRefHU ) ! IN/OUT
118
119 case ('height')
120 if ( .not. gsv_isAllocated(stateVectorRefHeight) ) then
121 ! initialize stateVectorRefHeight on analysis grid
122 call gsv_allocate(stateVectorRefHeight, tim_nstepobsinc, hco_anl, &
123 vco_anl, dateStamp_opt=tim_getDateStamp(), &
124 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
125 hInterpolateDegree_opt='LINEAR', &
126 varNames_opt=&
127 (/'Z_T','Z_M','P_T','P_M','TT','HU','P0'/) )
128 call gsv_zero(stateVectorRefHeight)
129 end if
130
131 ! initialize statevector_noZnoP on analysis grid
132 call gsv_allocate(statevector_noZnoP, tim_nstepobsinc, hco_anl, vco_anl, &
133 dateStamp_opt=tim_getDateStamp(), &
134 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
135 hInterpolateDegree_opt='LINEAR', &
136 varNames_opt=(/'TT','HU','P0'/))
137 write(*,*) 'gvt_setupRefFromTrialFiles: statevector_noZnoP allocated'
138
139 ! read trial files using default horizontal interpolation degree
140 call gio_readTrials( statevector_noZnoP ) ! IN/OUT
141
142 ! copy the statevectors
143 call gsv_copy(statevector_noZnoP, stateVectorRefHeight, &
144 allowVarMismatch_opt=.true. )
145
146 call gsv_deallocate(statevector_noZnoP)
147
148 ! do height/P calculation of the grid
149 call czp_calcZandP_nl( stateVectorRefHeight )
150
151 case default
152 if ( present(varKind_opt) ) then
153 if (varKind_opt == 'CH' .and. &
154 vnl_varKindFromVarname(varName) == varKind_opt ) then
155
156 varIndex = vnl_varListIndex(varName)
157
158 ! initialize stateVectorTrialvarKindCH(varIndex) on analysis grid
159
160 call gsv_allocate(stateVectorTrialvarKindCH(varIndex), &
161 tim_nstepobsinc, hco_anl, vco_anl, &
162 dateStamp_opt=tim_getDateStamp(), &
163 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
164 hInterpolateDegree_opt='LINEAR', &
165 varNames_opt=(/trim(varName),'P0'/) )
166
167 ! read trial files using default horizontal interpolation degree
168 call gio_readTrials( stateVectorTrialvarKindCH(varIndex) ) ! IN/OUT
169
170 varKindCHTrialsInitialized(varIndex) = .true.
171
172 else
173 call utl_abort('gvt_setupRefFromTrialFiles: unknown variable ='//trim(varName))
174 end if
175 else
176 call utl_abort('gvt_setupRefFromTrialFiles: unknown variable ='//trim(varName))
177 end if
178 end select
179
180 end subroutine gvt_setupRefFromTrialFiles
181
182 !--------------------------------------------------------------------------
183 ! gvt_transform_gsv
184 !--------------------------------------------------------------------------
185 subroutine gvt_transform_gsv(statevector, transform, statevectorOut_opt, &
186 stateVectorRef_opt, varName_opt, &
187 allowOverWrite_opt, maxBoxSize_opt, subgrid_opt )
188 !
189 !:Purpose: Top-level switch routine for transformations on the grid.
190 !
191 implicit none
192
193 ! Arguments:
194 type(struct_gsv), intent(inout) :: statevector ! statevector operand of the transformation
195 character(len=*), intent(in) :: transform ! string identifying the requested transformation
196 type(struct_gsv), optional, intent(inout) :: statevectorOut_opt
197 type(struct_gsv), optional, intent(in) :: statevectorRef_opt ! reference statevector necessary for some transformation
198 logical, optional, intent(in) :: allowOverWrite_opt
199 character(len=*), optional, intent(in) :: varName_opt ! additional variable/type info mandatory for some transformation
200 integer, optional, intent(in) :: maxBoxSize_opt ! additional info required by SSTSpread
201 character(len=*), optional, intent(in) :: subgrid_opt ! additional info required by SSTSpread to spread SST values
202
203 ! check stateVector and statevectorOut_opt are on the same grid
204 if ( present(stateVectorRef_opt) ) then
205 if ( .not. hco_equal(stateVectorRef_opt%hco, stateVector%hco) .or. &
206 .not. vco_equal(stateVectorRef_opt%vco, stateVector%vco) ) then
207 call utl_abort('gvt_transform: stateVectorRef_opt not on same grid as stateVector')
208 end if
209 end if
210
211 select case(trim(transform))
212
213 case ('AllTransformedToModel') ! Do all transformed variables: LPRtoPR, LVIStoVIS
214 if ( gsv_varExist(statevector,'LPR') ) then
215 write(*,*) 'gvt_transform: calling LPRtoPR stateVector transformation'
216 call LPRtoPR_gsv(statevector, statevectorOut_opt=statevectorOut_opt, &
217 allowOverWrite_opt=allowOverWrite_opt)
218 end if
219 if ( gsv_varExist(statevector,'LVIS') ) then
220 write(*,*) 'gvt_transform: calling LVIStoVIS stateVector transformation'
221 call LVIStoVIS(statevector, statevectorOut_opt=statevectorOut_opt, &
222 allowOverWrite_opt=allowOverWrite_opt)
223 end if
224
225 case ('UVtoVortDiv')
226 if (present(statevectorOut_opt)) then
227 call utl_abort('gvt_transform: for UVtoVortDiv, the option statevectorOut_opt is not yet available')
228 end if
229 call UVtoVortDiv_gsv(statevector)
230
231 case ('VortDivToPsiChi')
232 if (.not. gsv_varExist(statevector,'QR') .or. &
233 .not. gsv_varExist(statevector,'DD') ) then
234 call utl_abort('gvt_transform: for VortDivToPsiChi, variables QR and DD must be allocated in gridstatevector')
235 end if
236 if (present(statevectorOut_opt)) then
237 call utl_abort('gvt_transform: for VortDivToPsiChi, the option statevectorOut_opt is not yet available')
238 end if
239 call VortDivToPsiChi_gsv(statevector)
240
241 case ('UVtoPsiChi')
242 if (.not. gsv_varExist(statevector,'PP') .or. &
243 .not. gsv_varExist(statevector,'CC') ) then
244 call utl_abort('gvt_transform: for UVToPsiChi, variables PP and CC must be allocated in gridstatevector')
245 end if
246 if (present(statevectorOut_opt)) then
247 call utl_abort('gvt_transform: for UVToPsiChi, the option statevectorOut_opt is not yet available')
248 end if
249 call UVtoPsiChi_gsv(statevector)
250
251 case ('LQtoHU')
252 if ( .not. gsv_varExist(statevector,'HU') ) then
253 call utl_abort('gvt_transform: for LQtoHU, variable HU must be allocated in gridstatevector')
254 end if
255 if (present(statevectorOut_opt)) then
256 call utl_abort('gvt_transform: for LQtoHU, the option statevectorOut_opt is not yet available')
257 end if
258 call LQtoHU(statevector)
259
260 case ('HUtoLQ')
261 if ( .not. gsv_varExist(statevector,'HU') ) then
262 call utl_abort('gvt_transform: for HUtoLQ, variable HU must be allocated in gridstatevector')
263 end if
264 if (present(statevectorOut_opt)) then
265 call utl_abort('gvt_transform: for HUtoLQ, the option statevectorOut_opt is not yet available')
266 end if
267 call HUtoLQ_gsv(statevector)
268
269 case ('LQtoHU_tlm')
270 if ( .not. gsv_varExist(statevector,'HU') .and. &
271 .not. gsv_varExist(statevector,'LQ') ) then
272 call utl_abort('gvt_transform: for LQtoHU_tlm, variable HU or LQ must be allocated in gridstatevector')
273 end if
274 if ( gsv_varExist(statevector,'HU') .and. &
275 gsv_varExist(statevector,'LQ') ) then
276 call utl_abort('gvt_transform: for LQtoHU_tlm, only one of HU and LQ must be allocated in gridstatevector')
277 end if
278 if (present(statevectorOut_opt)) then
279 call utl_abort('gvt_transform: for LQtoHU_tlm, the option statevectorOut_opt is not yet available')
280 end if
281 call LQtoHU_tlm(statevector, stateVectorRef_opt)
282
283 case ('HUtoLQ_tlm')
284 if ( .not. gsv_varExist(statevector,'HU') ) then
285 call utl_abort('gvt_transform: for HUtoLQ_tlm, variable HU must be allocated in gridstatevector')
286 end if
287 if (present(statevectorOut_opt)) then
288 call utl_abort('gvt_transform: for HUtoLQ_tlm, the option statevectorOut_opt is not yet available')
289 end if
290 call HUtoLQ_tlm(statevector, stateVectorRef_opt)
291
292 case ('LQtoHU_ad')
293 if ( .not. gsv_varExist(statevector,'HU') ) then
294 call utl_abort('gvt_transform: for LQtoHU_ad, variable HU must be allocated in gridstatevector')
295 end if
296 if (present(statevectorOut_opt)) then
297 call utl_abort('gvt_transform: for LQtoHU_ad, the option statevectorOut_opt is not yet available')
298 end if
299 call LQtoHU_tlm(statevector, stateVectorRef_opt) ! self-adjoint
300
301 case ('HUtoLQ_ad')
302 if ( .not. gsv_varExist(statevector,'HU') ) then
303 call utl_abort('gvt_transform: for HUtoLQ_ad, variable HU must be allocated in gridstatevector')
304 end if
305 if (present(statevectorOut_opt)) then
306 call utl_abort('gvt_transform: for HUtoLQ_ad, the option statevectorOut_opt is not yet available')
307 end if
308 call HUtoLQ_tlm(statevector, stateVectorRef_opt) ! self-adjoint
309
310 case ('LPRoPR')
311 if ( .not. gsv_varExist(statevector,'PR') ) then
312 call utl_abort('gvt_transform: for LPRtoPR, variable PR must be allocated in gridstatevector')
313 end if
314 call LPRtoPR_gsv(statevector, statevectorOut_opt, &
315 allowOverWrite_opt=allowOverWrite_opt)
316
317 case ('PRtoLPR')
318 if ( .not. gsv_varExist(statevector,'PR') ) then
319 call utl_abort('gvt_transform: for PRtoLPR, variable PR must be allocated in gridstatevector')
320 end if
321 call PRtoLPR_gsv(statevector, statevectorOut_opt)
322
323 case ('LVIStoVIS')
324 if (present(statevectorOut_opt)) then
325 if ( .not. gsv_varExist(statevector,'LVIS')) then
326 call utl_abort('gvt_transform: for LVIStoVIS, variable LVIS must be allocated in gridstatevector IN')
327 end if
328 if ( .not. gsv_varExist(statevectorOut_opt,'VIS')) then
329 call utl_abort('gvt_transform: for LVIStoVIS, variable VIS must be allocated in gridstatevector OUT')
330 end if
331 call LVIStoVIS(statevector, statevectorOut_opt=statevectorOut_opt)
332 else
333 if (.not. gsv_varExist(statevector,'VIS') .or. &
334 .not. gsv_varExist(statevector,'LVIS') ) then
335 call utl_abort('gvt_transform: for LVIStoVIS, variables LVIS and VIS must be allocated in gridstatevector')
336 end if
337 call LVIStoVIS(statevector)
338 end if
339
340 case ('ZandP_nl')
341 if ( present(statevectorOut_opt) ) then
342 call utl_abort('gvt_transform: for ZandP_nl, the option statevectorOut_opt is not yet available')
343 end if
344 call czp_calcZandP_nl(stateVector)
345
346 case ('ZandP_tl')
347 if ( present(statevectorOut_opt) ) then
348 call utl_abort('gvt_transform: for ZandP_tl, the option statevectorOut_opt is not yet available')
349 end if
350 call ZandP_tl(stateVector)
351
352 case ('ZandP_ad')
353 if ( present(statevectorOut_opt) ) then
354 call utl_abort('gvt_transform: for ZandP_ad, the option statevectorOut_opt is not yet available')
355 end if
356 call ZandP_ad(stateVector)
357
358 case ('expCH_tlm')
359 if ( .not. gsv_varKindExist('CH') ) then
360 call utl_abort('gvt_transform: for expCH_tlm, variables of CH kind must be allocated in gridstatevector')
361 else if ( .not.present(varName_opt) ) then
362 call utl_abort('gvt_transform: for expCH_tlm, missing variable name')
363 else if ( .not. gsv_varExist(statevector,trim(varName_opt)) ) then
364 call utl_abort('gvt_transform: for expCH_tlm, variable '//trim(varName_opt)//' must be allocated in gridstatevector')
365 else if ( vnl_varKindFromVarname(trim(varName_opt)) /= 'CH' ) then
366 call utl_abort('gvt_transform: Invalid kind of varName for expCH_tlm')
367 else if ( .not. present(statevectorRef_opt)) then
368 call utl_abort('gvt_transform: for expCH_tlm, the option statevectorRef_opt must be present')
369 else if (present(statevectorOut_opt)) then
370 call utl_abort('gvt_transform: for expCH_tlm, the option statevectorOut_opt is not yet available')
371 end if
372 call expCH_tlm(statevector, varName_opt, stateVectorRef_opt)
373
374 case ('expCH_ad')
375 if ( .not. gsv_varKindExist('CH') ) then
376 call utl_abort('gvt_transform: for expCH_ad, variables of CH kind must be allocated in gridstatevector')
377 else if ( .not.present(varName_opt) ) then
378 call utl_abort('gvt_transform: for expCH_ad missing variable name')
379 else if ( .not. gsv_varExist(statevector,trim(varName_opt)) ) then
380 call utl_abort('gvt_transform: for expCH_ad, variable '//trim(varName_opt)//' must be allocated in gridstatevector')
381 else if ( vnl_varKindFromVarname(trim(varName_opt)) /= 'CH' ) then
382 call utl_abort('gvt_transform: Invalid kind of varName for expCH_ad')
383 else if ( .not. present(statevectorRef_opt)) then
384 call utl_abort('gvt_transform: for expCH_ad, the option statevectorRef_opt must be present')
385 else if (present(statevectorOut_opt)) then
386 call utl_abort('gvt_transform: for expCH_ad, the option statevectorOut_opt is not yet available')
387 end if
388 call expCH_tlm(statevector, varName_opt, stateVectorRef_opt) ! self-adjoint
389
390 case ('CH_bounds')
391 if ( .not. gsv_varKindExist('CH') ) then
392 call utl_abort('gvt_transform: for CH_bounds, variables of CH kind must be allocated in gridstatevector')
393 end if
394 call CH_bounds(statevector)
395
396 case ('oceanIceContinuous')
397 if (.not. present(statevectorRef_opt)) then
398 call utl_abort('gvt_transform: for gvt_oceanIceContinuous, the option statevectorRef_opt must be present')
399 end if
400 if (.not.present(varName_opt)) then
401 call utl_abort('gvt_transform: for gvt_oceanIceContinuous, missing variable name')
402 end if
403 call gvt_oceanIceContinuous(statevector, stateVectorRef_opt, varName_opt)
404
405 case ('SSTSpread')
406 if (.not.present(varName_opt)) then
407 call utl_abort('gvt_transform: for gvt_SSTSpread, missing variable name')
408 end if
409 if (.not.present(maxBoxSize_opt)) then
410 call utl_abort('gvt_transform: for gvt_SSTSpread, missing max box size')
411 end if
412 if (.not.present(subgrid_opt)) then
413 call utl_abort('gvt_transform: for gvt_SSTSpread, missing subgrid')
414 end if
415 call gvt_SSTSpread(statevector, varName_opt, maxBoxSize_opt, subgrid_opt)
416
417 case default
418 write(*,*)
419 write(*,*) 'Unsupported function : ', trim(transform)
420 call utl_abort('gvt_transform')
421 end select
422
423 end subroutine gvt_transform_gsv
424
425 !--------------------------------------------------------------------------
426 ! gvt_transform_ens
427 !--------------------------------------------------------------------------
428 subroutine gvt_transform_ens( ens,transform, allowOverWrite_opt, &
429 varName_opt, huMinValue_opt)
430 !
431 !:Purpose: Top-level switch routine for ensemble transformations on the
432 ! grid
433 !
434 implicit none
435
436 ! Arguments:
437 type(struct_ens), intent(inout) :: ens ! operand (ensemble of statevector)
438 character(len=*), intent(in) :: transform ! string identifying the requested transformation
439 logical, optional, intent(in) :: allowOverWrite_opt
440 character(len=*), optional, intent(in) :: varName_opt ! additional variable/type information mandatory for some transformation
441 real(8), optional, intent(in) :: huMinValue_opt ! HU min value for 'HUtoLQ' transformation
442
443 select case(trim(transform))
444 case ('AllTransformedToModel') ! Do all transformed variables: LPRtoPR
445 if ( ens_varExist(ens,'LPR') ) then
446 write(*,*) 'gvt_transform: calling LPRtoPR ensemble transformation'
447 call LPRtoPR_ens(ens, allowOverWrite_opt=allowOverWrite_opt)
448 end if
449 case ('HUtoLQ')
450 call HUtoLQ_ens(ens, huMinValue_opt=huMinValue_opt)
451 case ('LPRtoPR')
452 call LPRtoPR_ens(ens, allowOverWrite_opt=allowOverWrite_opt)
453 case ('UVtoPsiChi')
454 call UVtoPsiChi_ens(ens)
455 case ('UVtoVortDiv')
456 call UVtoVortDiv_ens(ens)
457 case ('logCH')
458 if ( .not.present(varName_opt) ) then
459 call utl_abort('gvt_transform: for logCH missing variable name')
460 else if ( vnl_varKindFromVarname(trim(varName_opt)) /= 'CH' ) then
461 call utl_abort('gvt_transform: Invalid kind of varName for logCH')
462 end if
463 call logCH_ens(ens,varName_opt)
464 case default
465 call utl_abort('gvt_transform_ens: Unsupported function '//trim(transform))
466 end select
467
468 end subroutine gvt_transform_ens
469
470 !--------------------------------------------------------------------------
471 ! gvt_getStateVectorTrial
472 !--------------------------------------------------------------------------
473 function gvt_getStateVectorTrial(varName) result(statevector_ptr)
474 !
475 !:Purpose: Returns a pointer to requested reference statevector
476 !
477 implicit none
478
479 ! Arguments:
480 character(len=*), intent(in) :: varName ! reference variable/type requested
481 ! Result:
482 type(struct_gsv), pointer :: statevector_ptr
483
484 select case ( trim(varName) )
485 case ('HU')
486 if ( .not. gsv_containsNonZeroValues(stateVectorRefHU) ) then
487 call utl_abort('gvt_getStateVectorTrial: do trials to stateVectorRefHU transform at higher level')
488 end if
489 statevector_ptr => stateVectorRefHU
490
491 case ('height')
492 if ( .not. gsv_containsNonZeroValues(stateVectorRefHeight) ) then
493 call utl_abort('gvt_getStateVectorTrial: do trials to stateVectorRefHeight transform at higher level')
494 end if
495 statevector_ptr => stateVectorRefHeight
496
497 case default
498 call utl_abort('gvt_getStateVectorTrial: unknown variable ='//trim(varName))
499 end select
500
501 end function gvt_getStateVectorTrial
502
503 !--------------------------------------------------------------------------
504 ! gvt_setupRefFromStateVector
505 !--------------------------------------------------------------------------
506 subroutine gvt_setupRefFromStateVector( stateVectorOnTrlGrid, varName, &
507 applyLimitOnHU_opt )
508 !
509 !:Purpose: computing the reference stateVector on the analysis grid at each
510 ! outer-loop iteration. The calculation is skipped if stateVectorRef* is
511 ! initialized (gsv_containsNonZeroValue(stateVectorRef*)=.true.).
512 ! The input stateVector is the high spatial/temporal resolution
513 ! statevector used for reading the trials and should contain
514 ! TT/HU/P0 if stateVectorRefHeight is asked for.
515 !
516 implicit none
517
518 ! Arguments:
519 type(struct_gsv), intent(in) :: stateVectorOnTrlGrid ! high spatial/temporal resolution statevector
520 character(len=*), intent(in) :: varName ! reference variable/type used
521 logical, optional, intent(in) :: applyLimitOnHU_opt
522
523 ! Locals:
524 type(struct_gsv) :: stateVectorLowResTime
525 type(struct_gsv) :: stateVectorLowResTimeSpace
526 type(struct_gsv), target :: stateVectorRefHUTT
527 logical :: allocHeightSfc
528 character(len=4), pointer :: varNames(:)
529
530 if ( mmpi_myid == 0 ) write(*,*) 'gvt_setupRefFromStateVector: START ', trim(varName)
531
532 if ( .not. associated(hco_trl) ) hco_trl => gsv_getHco(stateVectorOnTrlGrid)
533 if ( .not. associated(vco_trl) ) vco_trl => gsv_getVco(stateVectorOnTrlGrid)
534
535 select case ( trim(varName) )
536 case ('HU')
537 if ( .not. present(applyLimitOnHU_opt) ) then
538 call utl_abort('gvt_setupRefFromStateVector: applyLimitOnHU_opt for RefHU missing')
539 end if
540
541 if ( .not. gsv_isAllocated(stateVectorRefHU) ) then
542 call gsv_allocate(stateVectorRefHU, tim_nstepobsinc, hco_anl, vco_anl,&
543 dateStamp_opt=tim_getDateStamp(), &
544 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
545 hInterpolateDegree_opt='LINEAR', &
546 varNames_opt=(/'HU','P0'/) )
547 else
548 call gsv_zero( stateVectorRefHU )
549 end if
550
551 allocHeightSfc = ( vco_trl%Vcode /= 0 )
552
553 call gsv_allocate(stateVectorRefHUTT, tim_nstepobsinc, hco_anl, vco_anl,&
554 dateStamp_opt=tim_getDateStamp(),&
555 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
556 hInterpolateDegree_opt='LINEAR', &
557 varNames_opt=(/'HU','TT','P0'/) )
558
559 ! First, degrade the time steps
560 call gsv_allocate(stateVectorLowResTime, tim_nstepobsinc, hco_trl, &
561 vco_trl, dataKind_opt=pre_incrReal, &
562 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true.,&
563 allocHeightSfc_opt=allocHeightSfc, &
564 hInterpolateDegree_opt='LINEAR', &
565 allocHeight_opt=.false., allocPressure_opt=.false. )
566 call gsv_copy( stateVectorOnTrlGrid, stateVectorLowResTime, &
567 allowTimeMismatch_opt=.true., allowVarMismatch_opt=.true. )
568
569 ! Second, interpolate to the low-resolution spatial grid.
570 nullify(varNames)
571 call gsv_varNamesList(varNames, stateVectorLowResTime)
572 call gsv_allocate(stateVectorLowResTimeSpace, tim_nstepobsinc, hco_anl, &
573 vco_anl, dataKind_opt=pre_incrReal, &
574 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true.,&
575 allocHeightSfc_opt=.true., &
576 hInterpolateDegree_opt='LINEAR', varNames_opt=varNames)
577 call int_interp_gsv(stateVectorLowResTime, stateVectorLowResTimeSpace)
578
579 ! Now copy only P0, HU, and TT to create reference stateVector.
580 call gsv_copy(stateVectorLowResTimeSpace, stateVectorRefHUTT, &
581 allowTimeMismatch_opt=.false., allowVarMismatch_opt=.true.)
582 call gsv_deallocate(stateVectorLowResTimeSpace)
583 call gsv_deallocate(stateVectorLowResTime)
584
585 ! Impose humidity limits on stateVectorRefHUTT
586 if ( applyLimitOnHU_opt ) then
587 write(*,*) 'var: impose limits on stateVectorRefHUTT'
588 call qlim_saturationLimit(stateVectorRefHUTT)
589 call qlim_rttovLimit(stateVectorRefHUTT)
590 end if
591
592 call gsv_copy(stateVectorRefHUTT, stateVectorRefHU, &
593 allowTimeMismatch_opt=.false., allowVarMismatch_opt=.true.)
594 call gsv_deallocate(stateVectorRefHUTT)
595
596 case ('height')
597 if ( .not. gsv_isAllocated(stateVectorRefHeight) ) then
598 call gsv_allocate( stateVectorRefHeight, tim_nstepobsinc, hco_anl, &
599 vco_anl, dateStamp_opt=tim_getDateStamp(), &
600 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
601 hInterpolateDegree_opt='LINEAR', &
602 varNames_opt=&
603 (/'Z_T','Z_M','P_T','P_M','TT','HU','P0','P0LS'/) )
604 else
605 call gsv_zero( stateVectorRefHeight )
606 end if
607
608 ! First, degrade the time steps
609 nullify(varNames)
610 call gsv_varNamesList(varNames, stateVectorRefHeight)
611 call gsv_allocate(stateVectorLowResTime, tim_nstepobsinc, hco_trl, &
612 vco_trl, dateStamp_opt=tim_getDateStamp(), &
613 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
614 hInterpolateDegree_opt='LINEAR', &
615 varNames_opt=varNames )
616 call gsv_copy(stateVectorOnTrlGrid, stateVectorLowResTime, &
617 allowTimeMismatch_opt=.true., allowVarMismatch_opt=.true.)
618
619 ! Second, interpolate to the low-resolution spatial grid.
620 nullify(varNames)
621 call gsv_varNamesList(varNames, stateVectorLowResTime)
622 call gsv_allocate(stateVectorLowResTimeSpace, tim_nstepobsinc, hco_anl,&
623 vco_anl, dateStamp_opt=tim_getDateStamp(), &
624 mpi_local_opt=.true., allocHeightSfc_opt=.true., &
625 hInterpolateDegree_opt='LINEAR', varNames_opt=varNames)
626 call int_interp_gsv(stateVectorLowResTime, stateVectorLowResTimeSpace)
627
628 ! Now copy to create final stateVector height.
629 call gsv_copy(stateVectorLowResTimeSpace, stateVectorRefHeight, &
630 allowTimeMismatch_opt=.false., allowVarMismatch_opt=.true.)
631 call gsv_deallocate(stateVectorLowResTimeSpace)
632 call gsv_deallocate(stateVectorLowResTime)
633
634 ! do height/P calculation of the grid
635 call czp_calcZandP_nl(stateVectorRefHeight)
636
637 end select
638
639 if ( mmpi_myid == 0 ) write(*,*) 'gvt_setupRefFromStateVector: END'
640
641 end subroutine gvt_setupRefFromStateVector
642
643 !--------------------------------------------------------------------------
644 ! LQtoHU
645 !--------------------------------------------------------------------------
646 subroutine LQtoHU(statevector)
647 !
648 !:Purpose: Specific humidity logarithm exponentiation.
649 !
650 implicit none
651
652 ! Arguments:
653 type(struct_gsv), intent(inout) :: statevector
654
655 ! Locals:
656 integer :: lonIndex,latIndex,levIndex,stepIndex
657 real(8), pointer :: hu_ptr(:,:,:,:), lq_ptr(:,:,:,:)
658
659 call gsv_getField(statevector,hu_ptr,'HU')
660 call gsv_getField(statevector,lq_ptr,'HU')
661
662 do stepIndex = 1, statevector%numStep
663 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
664 do levIndex = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
665 do latIndex = statevector%myLatBeg, statevector%myLatEnd
666 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
667 hu_ptr(lonIndex,latIndex,levIndex,stepIndex) = &
668 exp(lq_ptr(lonIndex,latIndex,levIndex,stepIndex))
669 end do
670 end do
671 end do
672 !$OMP END PARALLEL DO
673 end do
674
675 end subroutine LQtoHU
676
677 !--------------------------------------------------------------------------
678 ! HUtoLQ_gsv
679 !--------------------------------------------------------------------------
680 subroutine HUtoLQ_gsv(statevector)
681 !
682 !:Purpose: Logarithmic transformation of specific humidity
683 !
684 implicit none
685
686 ! Arguments:
687 type(struct_gsv), intent(inout) :: statevector
688
689 ! Locals:
690 integer :: lonIndex,latIndex,levIndex,stepIndex
691 real(4), pointer :: hu_ptr_r4(:,:,:,:), lq_ptr_r4(:,:,:,:)
692 real(8), pointer :: hu_ptr_r8(:,:,:,:), lq_ptr_r8(:,:,:,:)
693
694 if ( gsv_getDataKind(statevector) == 8 ) then
695
696 call gsv_getField(statevector,hu_ptr_r8,'HU')
697 call gsv_getField(statevector,lq_ptr_r8,'HU')
698
699 do stepIndex = 1, statevector%numStep
700 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
701 do levIndex = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
702 do latIndex = statevector%myLatBeg, statevector%myLatEnd
703 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
704 lq_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
705 log(max(hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex),&
706 gsv_rhumin) )
707 end do
708 end do
709 end do
710 !$OMP END PARALLEL DO
711 end do
712
713 else
714
715 call gsv_getField(statevector,hu_ptr_r4,'HU')
716 call gsv_getField(statevector,lq_ptr_r4,'HU')
717
718 do stepIndex = 1, statevector%numStep
719 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
720 do levIndex = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
721 do latIndex = statevector%myLatBeg, statevector%myLatEnd
722 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
723 lq_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
724 log(max(hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex),&
725 real(gsv_rhumin,4)) )
726 end do
727 end do
728 end do
729 !$OMP END PARALLEL DO
730 end do
731
732 end if
733
734 end subroutine HUtoLQ_gsv
735
736 !--------------------------------------------------------------------------
737 ! HUtoLQ_ens
738 !--------------------------------------------------------------------------
739 subroutine HUtoLQ_ens(ens, huMinValue_opt)
740 !
741 !:Purpose: Specific humidity logarithm exponentiation (ensemble processing)
742 !
743 implicit none
744
745 ! Arguments:
746 type(struct_ens), intent(inout) :: ens
747 real(8), optional, intent(in) :: huMinValue_opt
748
749 ! Locals:
750 integer :: lonIndex, latIndex, levIndex, stepIndex, memberIndex
751 integer :: myLatBeg, myLatEnd, myLonBeg, myLonEnd
752 character(len=4) :: varName
753 real(4), pointer :: ptr4d_r4(:,:,:,:)
754 real(4) :: huMinValue
755
756 if (present(huMinValue_opt)) then
757 huMinValue = huMinValue_opt
758 else
759 huMinValue = MPC_MINIMUM_HU_R4
760 end if
761
762 call ens_getLatLonBounds(ens, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
763
764 do levIndex = 1, ens_getNumK(ens)
765
766 varName = ens_getVarNameFromK(ens,levIndex)
767 if (varName /= 'HU') cycle
768
769 ptr4d_r4 => ens_getOneLev_r4(ens,levIndex)
770
771 do latIndex = myLatBeg, myLatEnd
772 do lonIndex = myLonBeg, myLonEnd
773 do stepIndex = 1, ens_getNumStep(ens)
774 do memberIndex = 1, ens_getNumMembers(ens)
775 ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
776 log(max( ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex),&
777 huMinValue) )
778 end do
779 end do
780 end do
781 end do
782
783 end do
784
785 end subroutine HUtoLQ_ens
786
787 !--------------------------------------------------------------------------
788 ! LQtoHU_tlm
789 !--------------------------------------------------------------------------
790 subroutine LQtoHU_tlm(statevector, stateVectorRef_opt)
791 !
792 !:Purpose: Tangent linear of exponentiation transformation of specific
793 ! humidity in logarithmic representation
794 !
795 implicit none
796
797 ! Arguments:
798 type(struct_gsv), intent(inout) :: statevector
799 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
800
801 ! Locals:
802 integer :: lonIndex,latIndex,levIndex,stepIndex
803 real(8), pointer :: hu_ptr_r8(:,:,:,:), lq_ptr_r8(:,:,:,:)
804 real(8), pointer :: hu_trial(:,:,:,:)
805 real(4), pointer :: hu_ptr_r4(:,:,:,:), lq_ptr_r4(:,:,:,:)
806 character(len=4) :: varName
807
808 if ( gsv_varExist(statevector,'HU') ) then
809 varName = 'HU'
810 else if ( gsv_varExist(statevector,'LQ') ) then
811 varName = 'LQ'
812 else
813 call utl_abort('LQtoHU_tlm: either HU or LQ must be part of stateVector')
814 end if
815
816 if ( present(statevectorRef_opt) ) then
817 call gsv_getField(stateVectorRef_opt,hu_trial,'HU')
818 else
819 if ( .not. gsv_containsNonZeroValues(stateVectorRefHU) ) then
820 call utl_abort('LQtoHU_tlm: do trials to stateVectorRefHU transform at higher level')
821 end if
822 call gsv_getField(stateVectorRefHU,hu_trial,'HU')
823 end if
824
825 if ( gsv_getDataKind(statevector) == 4 ) then
826 call gsv_getField(statevector,hu_ptr_r4,varName)
827 call gsv_getField(statevector,lq_ptr_r4,varName)
828
829 do stepIndex = 1, statevector%numStep
830 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
831 do levIndex = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
832 do latIndex = statevector%myLatBeg, statevector%myLatEnd
833 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
834 hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
835 lq_ptr_r4(lonIndex,latIndex,levIndex,stepIndex)* &
836 max( hu_trial(lonIndex,latIndex,levIndex,stepIndex),&
837 MPC_MINIMUM_HU_R4)
838 end do
839 end do
840 end do
841 !$OMP END PARALLEL DO
842 end do
843 else
844 call gsv_getField(statevector,hu_ptr_r8,varName)
845 call gsv_getField(statevector,lq_ptr_r8,varName)
846
847 do stepIndex = 1, statevector%numStep
848 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
849 do levIndex = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
850 do latIndex = statevector%myLatBeg, statevector%myLatEnd
851 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
852 hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
853 lq_ptr_r8(lonIndex,latIndex,levIndex,stepIndex)* &
854 max( hu_trial(lonIndex,latIndex,levIndex,stepIndex),&
855 MPC_MINIMUM_HU_R8)
856 end do
857 end do
858 end do
859 !$OMP END PARALLEL DO
860 end do
861 end if
862
863 end subroutine LQtoHU_tlm
864
865 !--------------------------------------------------------------------------
866 ! HUtoLQ_tlm
867 !--------------------------------------------------------------------------
868 subroutine HUtoLQ_tlm(statevector, stateVectorRef_opt)
869 !
870 !:Purpose: Tangent linear of logarithmic transformation of specific
871 ! humidity
872 !
873 implicit none
874
875 ! Arguments:
876 type(struct_gsv), intent(inout) :: statevector
877 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
878
879 ! Locals:
880 integer :: lonIndex,latIndex,levIndex,stepIndex
881 real(8), pointer :: hu_ptr_r8(:,:,:,:), lq_ptr_r8(:,:,:,:)
882 real(8), pointer :: hu_trial(:,:,:,:)
883 real(4), pointer :: hu_ptr_r4(:,:,:,:), lq_ptr_r4(:,:,:,:)
884
885 if ( present(statevectorRef_opt) ) then
886 call gsv_getField(stateVectorRef_opt,hu_trial,'HU')
887 else
888 if ( .not. gsv_containsNonZeroValues(stateVectorRefHU) ) then
889 call utl_abort('HUtoLQ_tlm: do trials to stateVectorRefHU transform at higher level')
890 end if
891 call gsv_getField(stateVectorRefHU,hu_trial,'HU')
892 end if
893
894 if ( gsv_getDataKind(statevector) == 4 ) then
895 call gsv_getField(statevector,hu_ptr_r4,'HU')
896 call gsv_getField(statevector,lq_ptr_r4,'HU')
897
898 do stepIndex = 1, statevector%numStep
899 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
900 do levIndex = 1,gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
901 do latIndex = statevector%myLatBeg, statevector%myLatEnd
902 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
903 lq_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
904 hu_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) / &
905 max( hu_trial(lonIndex,latIndex,levIndex,stepIndex),&
906 MPC_MINIMUM_HU_R8)
907 end do
908 end do
909 end do
910 !$OMP END PARALLEL DO
911 end do
912 else
913 call gsv_getField(statevector,hu_ptr_r8,'HU')
914 call gsv_getField(statevector,lq_ptr_r8,'HU')
915
916 do stepIndex = 1, statevector%numStep
917 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
918 do levIndex = 1,gsv_getNumLev(statevector,vnl_varLevelFromVarname('HU'))
919 do latIndex = statevector%myLatBeg, statevector%myLatEnd
920 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
921 lq_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
922 hu_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) / &
923 max( hu_trial(lonIndex,latIndex,levIndex,stepIndex),&
924 MPC_MINIMUM_HU_R8)
925 end do
926 end do
927 end do
928 !$OMP END PARALLEL DO
929 end do
930 end if
931
932 end subroutine HUtoLQ_tlm
933
934 !--------------------------------------------------------------------------
935 ! LPRtoPR_gsv
936 !--------------------------------------------------------------------------
937 subroutine LPRtoPR_gsv(stateVector, statevectorOut_opt, allowOverWrite_opt)
938 !
939 !:Purpose: Quantity of precipitation logarithm exponentiation
940 !
941 implicit none
942
943 ! Arguments:
944 type(struct_gsv), intent(inout) :: stateVector
945 type(struct_gsv), optional, intent(inout) :: statevectorOut_opt
946 logical, optional, intent(in) :: allowOverWrite_opt
947
948 ! Locals:
949 integer :: lonIndex,latIndex,levIndex,stepIndex
950 logical :: overWriteNeeded
951 real(4), pointer :: pr_ptr_r4(:,:,:,:), lpr_ptr_r4(:,:,:,:)
952 real(8), pointer :: pr_ptr_r8(:,:,:,:), lpr_ptr_r8(:,:,:,:)
953
954 ! Check if overWrite of PR is needed, but not allowed
955 overWriteNeeded = .false.
956 if (present(statevectorOut_opt)) then
957 if (.not. gsv_varExist(stateVectorOut_opt,'PR')) then
958 overWriteNeeded = .true.
959 end if
960 else
961 if (.not. gsv_varExist(stateVector,'PR')) then
962 overWriteNeeded = .true.
963 end if
964 end if
965 if (overWriteNeeded) then
966 if(present(allowOverWrite_opt)) then
967 if (.not.allowOverWrite_opt) then
968 call utl_abort('LPRtoPR_gsv: allowOverWrite_opt is false, but PR not present in stateVector')
969 end if
970 else
971 call utl_abort('LPRtoPR_gsv: allowOverWrite_opt not specified, but PR not present in stateVector')
972 end if
973 end if
974
975 if ( gsv_getDataKind(stateVector) == 8 ) then
976
977 if (present(statevectorOut_opt)) then
978 if (gsv_varExist(stateVectorOut_opt,'PR')) then
979 call gsv_getField(statevectorOut_opt,pr_ptr_r8,'PR')
980 else
981 call gsv_getField(statevectorOut_opt,pr_ptr_r8,'LPR')
982 end if
983 else
984 if (gsv_varExist(stateVector,'PR')) then
985 call gsv_getField(stateVector,pr_ptr_r8,'PR')
986 else
987 call gsv_getField(stateVector,pr_ptr_r8,'LPR')
988 end if
989 end if
990 call gsv_getField(stateVector,lpr_ptr_r8,'LPR')
991
992 do stepIndex = 1, stateVector%numStep
993 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
994 do levIndex = 1, gsv_getNumLev(stateVector,vnl_varLevelFromVarname('LPR'))
995 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
996 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
997 if (lpr_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) > &
998 0.99D0*MPC_missingValue_R8) then
999
1000 pr_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
1001 max( 0.0D0, exp( lpr_ptr_r8(lonIndex,latIndex,levIndex,&
1002 stepIndex)) - &
1003 MPC_MINIMUM_PR_R8)
1004 else
1005 pr_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
1006 MPC_missingValue_R8
1007 end if
1008 end do
1009 end do
1010 end do
1011 !$OMP END PARALLEL DO
1012 end do
1013
1014 else
1015
1016 if (present(statevectorOut_opt)) then
1017 if (gsv_varExist(stateVectorOut_opt,'PR')) then
1018 call gsv_getField(statevectorOut_opt,pr_ptr_r4,'PR')
1019 else
1020 call gsv_getField(statevectorOut_opt,pr_ptr_r4,'LPR')
1021 end if
1022 else
1023 if (gsv_varExist(stateVector,'PR')) then
1024 call gsv_getField(stateVector,pr_ptr_r4,'PR')
1025 else
1026 call gsv_getField(stateVector,pr_ptr_r4,'LPR')
1027 end if
1028 end if
1029 call gsv_getField(stateVector,lpr_ptr_r4,'LPR')
1030
1031 do stepIndex = 1, stateVector%numStep
1032 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1033 do levIndex = 1,gsv_getNumLev(stateVector,vnl_varLevelFromVarname('LPR'))
1034 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
1035 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
1036 if (lpr_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) > &
1037 0.99*MPC_missingValue_R4) then
1038
1039 pr_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
1040 max( 0.0, exp( lpr_ptr_r4(lonIndex,latIndex,levIndex,&
1041 stepIndex)) - &
1042 MPC_MINIMUM_PR_R4)
1043 else
1044 pr_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
1045 MPC_missingValue_R4
1046 end if
1047 end do
1048 end do
1049 end do
1050 !$OMP END PARALLEL DO
1051 end do
1052
1053 end if
1054
1055 ! Change the variable name from LPR to PR if PR does not exist
1056 if (present(statevectorOut_opt)) then
1057 if (.not. gsv_varExist(stateVectorOut_opt,'PR')) then
1058 call gsv_modifyVarName(stateVectorOut_opt,'LPR','PR')
1059 end if
1060 else
1061 if (.not. gsv_varExist(stateVector,'PR')) then
1062 call gsv_modifyVarName(stateVector,'LPR','PR')
1063 end if
1064 end if
1065
1066 end subroutine LPRtoPR_gsv
1067
1068 !--------------------------------------------------------------------------
1069 ! LPRtoPR_ens
1070 !--------------------------------------------------------------------------
1071 subroutine LPRtoPR_ens(ens, allowOverWrite_opt)
1072 !
1073 !:Purpose: Quantity of precipitation logarithm exponentiation
1074 ! (ensemble processing)
1075 !
1076 implicit none
1077
1078 ! Arguments:
1079 type(struct_ens), intent(inout) :: ens
1080 logical, optional, intent(in) :: allowOverWrite_opt
1081
1082 ! Locals:
1083 integer :: lonIndex, latIndex, levIndex, stepIndex, memberIndex
1084 integer :: myLatBeg, myLatEnd, myLonBeg, myLonEnd, kIndexLPR, kIndexPR
1085 logical :: overWriteNeeded
1086 real(4), pointer :: PR_ptr_r4(:,:,:,:), LPR_ptr_r4(:,:,:,:)
1087
1088 ! Check if overWrite of PR is needed, but not allowed
1089 overWriteNeeded = .false.
1090 if (.not. ens_varExist(ens,'PR')) then
1091 overWriteNeeded = .true.
1092 end if
1093 if (overWriteNeeded) then
1094 if(present(allowOverWrite_opt)) then
1095 if (.not.allowOverWrite_opt) then
1096 call utl_abort('LPRtoPR_ens: allowOverWrite_opt is false, but PR not present in ensemble')
1097 end if
1098 else
1099 call utl_abort('LPRtoPR_ens: allowOverWrite_opt not specified, but PR not present in ensemble')
1100 end if
1101 end if
1102
1103 call ens_getLatLonBounds(ens, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
1104
1105 levIndex = 1
1106 kIndexLPR = ens_getKFromLevVarName(ens, levIndex, 'LPR')
1107 if (ens_varExist(ens,'PR')) then
1108 kIndexPR = ens_getKFromLevVarName(ens, levIndex, 'PR')
1109 else
1110 kIndexPR = ens_getKFromLevVarName(ens, levIndex, 'LPR')
1111 end if
1112
1113 LPR_ptr_r4 => ens_getOneLev_r4(ens,kIndexLPR)
1114 PR_ptr_r4 => ens_getOneLev_r4(ens,kIndexPR)
1115
1116 do latIndex = myLatBeg, myLatEnd
1117 do lonIndex = myLonBeg, myLonEnd
1118 do stepIndex = 1, ens_getNumStep(ens)
1119 do memberIndex = 1, ens_getNumMembers(ens)
1120 PR_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
1121 max( 0.0, exp( LPR_ptr_r4(memberIndex,stepIndex,lonIndex,&
1122 latIndex)) - &
1123 MPC_MINIMUM_PR_R4)
1124 ! maybe also impose minimum value of LPR
1125 !LPR_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
1126 ! max(log(MPC_MINIMUM_PR_R4), LPR_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex))
1127 end do
1128 end do
1129 end do
1130 end do
1131
1132 ! Change the variable name from LPR to PR if PR does not exist
1133 if (.not. ens_varExist(ens,'PR')) then
1134 call ens_modifyVarName(ens,'LPR','PR')
1135 end if
1136
1137 end subroutine LPRtoPR_ens
1138
1139 !--------------------------------------------------------------------------
1140 ! PRtoLPR_gsv
1141 !--------------------------------------------------------------------------
1142 subroutine PRtoLPR_gsv(stateVector, statevectorOut_opt)
1143 !
1144 !:Purpose: Logarithmic transformation of quantity of precipitation
1145 !
1146 implicit none
1147
1148 ! Arguments:
1149 type(struct_gsv), intent(inout) :: stateVector
1150 type(struct_gsv), optional, intent(inout) :: statevectorOut_opt
1151
1152 ! Locals:
1153 integer :: lonIndex,latIndex,levIndex,stepIndex
1154 real(4), pointer :: pr_ptr_r4(:,:,:,:), lpr_ptr_r4(:,:,:,:)
1155 real(8), pointer :: pr_ptr_r8(:,:,:,:), lpr_ptr_r8(:,:,:,:)
1156
1157 if ( gsv_getDataKind(stateVector) == 8 ) then
1158
1159 if (present(statevectorOut_opt)) then
1160 call gsv_getField(statevectorOut_opt,lpr_ptr_r8,'LPR')
1161 else
1162 call gsv_getField(stateVector,lpr_ptr_r8,'LPR')
1163 end if
1164 call gsv_getField(stateVector,pr_ptr_r8,'PR')
1165
1166 do stepIndex = 1, stateVector%numStep
1167 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1168 do levIndex = 1, gsv_getNumLev(stateVector,vnl_varLevelFromVarname('PR'))
1169 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
1170 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
1171 if (pr_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) > -1.0D0) then
1172 lpr_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
1173 log( MPC_MINIMUM_PR_R8 + &
1174 max(0.0d0,pr_ptr_r8(lonIndex,latIndex,levIndex,&
1175 stepIndex)))
1176 else
1177 lpr_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
1178 MPC_missingValue_R8
1179 end if
1180 end do
1181 end do
1182 end do
1183 !$OMP END PARALLEL DO
1184 end do
1185
1186 else
1187
1188 if (present(statevectorOut_opt)) then
1189 call gsv_getField(statevectorOut_opt,lpr_ptr_r4,'LPR')
1190 else
1191 call gsv_getField(stateVector,lpr_ptr_r4,'LPR')
1192 end if
1193 call gsv_getField(stateVector,pr_ptr_r4,'PR')
1194
1195 do stepIndex = 1, stateVector%numStep
1196 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1197 do levIndex = 1,gsv_getNumLev(stateVector,vnl_varLevelFromVarname('PR'))
1198 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
1199 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
1200 if (pr_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) > -1.0) then
1201 lpr_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
1202 log( MPC_MINIMUM_PR_R4 + &
1203 max(0.0,pr_ptr_r4(lonIndex,latIndex,levIndex,&
1204 stepIndex)))
1205 else
1206 lpr_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
1207 MPC_missingValue_R4
1208 end if
1209 end do
1210 end do
1211 end do
1212 !$OMP END PARALLEL DO
1213 end do
1214
1215 end if
1216
1217 end subroutine PRtoLPR_gsv
1218
1219 !--------------------------------------------------------------------------
1220 ! LVIStoVIS
1221 !--------------------------------------------------------------------------
1222 subroutine LVIStoVIS(stateVector, statevectorOut_opt, allowOverWrite_opt)
1223 !
1224 !:Purpose: Visibility logarithm exponentiation
1225 !
1226 implicit none
1227
1228 ! Arguments:
1229 type(struct_gsv), intent(inout) :: stateVector
1230 type(struct_gsv), optional, intent(inout) :: statevectorOut_opt
1231 logical, optional, intent(in) :: allowOverWrite_opt
1232
1233 ! Locals:
1234 integer :: lonIndex,latIndex,levIndex,stepIndex
1235 logical :: overWriteNeeded
1236 real(4), pointer :: vis_ptr_r4(:,:,:,:), lvis_ptr_r4(:,:,:,:)
1237 real(8), pointer :: vis_ptr_r8(:,:,:,:), lvis_ptr_r8(:,:,:,:)
1238
1239 ! Check if overWrite of PR is needed, but not allowed
1240 overWriteNeeded = .false.
1241 if (present(statevectorOut_opt)) then
1242 if (.not. gsv_varExist(stateVectorOut_opt,'PR')) then
1243 overWriteNeeded = .true.
1244 end if
1245 else
1246 if (.not. gsv_varExist(stateVector,'PR')) then
1247 overWriteNeeded = .true.
1248 end if
1249 end if
1250 if (overWriteNeeded) then
1251 if(present(allowOverWrite_opt)) then
1252 if (.not.allowOverWrite_opt) then
1253 call utl_abort('LVIStoVIS_gsv: allowOverWrite_opt is false, but PR not present in stateVector')
1254 end if
1255 else
1256 call utl_abort('LVIStoVIS_gsv: allowOverWrite_opt not specified, but PR not present in stateVector')
1257 end if
1258 end if
1259
1260 if ( gsv_getDataKind(stateVector) == 8 ) then
1261
1262 if (present(statevectorOut_opt)) then
1263 if (gsv_varExist(stateVectorOut_opt,'VIS')) then
1264 call gsv_getField(statevectorOut_opt,vis_ptr_r8,'VIS')
1265 else
1266 call gsv_getField(statevectorOut_opt,vis_ptr_r8,'LVIS')
1267 end if
1268 else
1269 if (gsv_varExist(stateVector,'VIS')) then
1270 call gsv_getField(stateVector,vis_ptr_r8,'VIS')
1271 else
1272 call gsv_getField(stateVector,vis_ptr_r8,'LVIS')
1273 end if
1274 end if
1275 call gsv_getField(stateVector,lvis_ptr_r8,'LVIS')
1276
1277 do stepIndex = 1, stateVector%numStep
1278 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1279 do levIndex = 1, gsv_getNumLev( stateVector,&
1280 vnl_varLevelFromVarname('LVIS'))
1281 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
1282 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
1283 vis_ptr_r8(lonIndex,latIndex,levIndex,stepIndex) = &
1284 exp(lvis_ptr_r8(lonIndex,latIndex,levIndex,stepIndex))
1285 end do
1286 end do
1287 end do
1288 !$OMP END PARALLEL DO
1289 end do
1290
1291 else
1292
1293 if (present(statevectorOut_opt)) then
1294 if (gsv_varExist(stateVectorOut_opt,'VIS')) then
1295 call gsv_getField(statevectorOut_opt,vis_ptr_r4,'VIS')
1296 else
1297 call gsv_getField(statevectorOut_opt,vis_ptr_r4,'LVIS')
1298 end if
1299 else
1300 if (gsv_varExist(stateVector,'VIS')) then
1301 call gsv_getField(stateVector,vis_ptr_r4,'VIS')
1302 else
1303 call gsv_getField(stateVector,vis_ptr_r4,'LVIS')
1304 end if
1305 end if
1306 call gsv_getField(stateVector,lvis_ptr_r4,'LVIS')
1307
1308 do stepIndex = 1, stateVector%numStep
1309 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1310 do levIndex = 1, gsv_getNumLev( stateVector,&
1311 vnl_varLevelFromVarname('LVIS'))
1312 do latIndex = stateVector%myLatBeg, stateVector%myLatEnd
1313 do lonIndex = stateVector%myLonBeg, stateVector%myLonEnd
1314 vis_ptr_r4(lonIndex,latIndex,levIndex,stepIndex) = &
1315 exp(lvis_ptr_r4(lonIndex,latIndex,levIndex,stepIndex))
1316 end do
1317 end do
1318 end do
1319 !$OMP END PARALLEL DO
1320 end do
1321
1322 end if
1323
1324 ! Change the variable name from LVIS to VIS if VIS does not exist
1325 if (present(statevectorOut_opt)) then
1326 if (.not. gsv_varExist(stateVectorOut_opt,'VIS')) then
1327 call gsv_modifyVarName(stateVectorOut_opt,'LVIS','VIS')
1328 end if
1329 else
1330 if (.not. gsv_varExist(stateVector,'VIS')) then
1331 call gsv_modifyVarName(stateVector,'LVIS','VIS')
1332 end if
1333 end if
1334
1335 end subroutine LVIStoVIS
1336
1337
1338 !--------------------------------------------------------------------------
1339 ! ZandP_tl
1340 !--------------------------------------------------------------------------
1341 subroutine ZandP_tl(stateVector)
1342 !
1343 !:Purpose: Tangeant linear of height and pressure computation.
1344 ! The computation order depends on the native model
1345 ! representation (height or pressure based).
1346 !
1347 implicit none
1348
1349 ! Arguments:
1350 type(struct_gsv), intent(inout) :: stateVector
1351
1352 call czp_calcZandP_tl(statevector, stateVectorRefHeight)
1353
1354 end subroutine ZandP_tl
1355
1356 !--------------------------------------------------------------------------
1357 ! ZandP_ad
1358 !--------------------------------------------------------------------------
1359 subroutine ZandP_ad(stateVector)
1360 !
1361 !:Purpose: Adjoint of the tangeant linear computation of both height and
1362 ! pressure. The computation order depends on the native model
1363 ! representation (height or pressure based).
1364 !
1365 implicit none
1366
1367 ! Arguments:
1368 type(struct_gsv), intent(inout) :: stateVector
1369
1370 call czp_calcZandP_ad(statevector, stateVectorRefHeight)
1371
1372 end subroutine ZandP_ad
1373
1374 !--------------------------------------------------------------------------
1375 ! UVtoVortDiv_gsv
1376 !--------------------------------------------------------------------------
1377 subroutine UVtoVortDiv_gsv(statevector)
1378 !
1379 !:Purpose: Wind components to relative vorticity and divergence transformation.
1380 !
1381 implicit none
1382
1383 ! Arguments:
1384 type(struct_gsv), intent(inout) :: statevector
1385
1386 ! Locals:
1387 real(8), pointer :: uu_ptr(:,:,:,:), vv_ptr(:,:,:,:)
1388 real(8), pointer :: qr_ptr(:,:,:,:), dd_ptr(:,:,:,:)
1389 integer :: stepIndex
1390 integer :: nlev_M
1391
1392 call gsv_getField(statevector,uu_ptr,'UU')
1393 call gsv_getField(statevector,vv_ptr,'VV')
1394 if (gsv_varExist(statevector,'QR') .and. &
1395 gsv_varExist(statevector,'DD')) then
1396 call gsv_getField(statevector,qr_ptr,'QR')
1397 call gsv_getField(statevector,dd_ptr,'DD')
1398 else
1399 call gsv_getField(statevector,qr_ptr,'UU')
1400 call gsv_getField(statevector,dd_ptr,'VV')
1401 end if
1402
1403 nlev_M = gsv_getNumLev (statevector,'MM')
1404
1405 if ( statevector%hco%global ) then
1406 call utl_abort('UVtoVortDiv_gsv: global mode not available')
1407 else
1408 do stepIndex = 1, statevector%numStep
1409 call lgt_UVToVortDiv(qr_ptr(:,:,:,stepIndex), dd_ptr(:,:,:,stepIndex), & ! OUT
1410 uu_ptr(:,:,:,stepIndex), vv_ptr(:,:,:,stepIndex), & ! IN
1411 nlev_M ) ! IN
1412 end do
1413 end if
1414
1415 end subroutine UVtoVortDiv_gsv
1416
1417 !--------------------------------------------------------------------------
1418 ! vortDivtoPsiChi_gsv
1419 !--------------------------------------------------------------------------
1420 subroutine vortDivToPsiChi_gsv(statevector)
1421 !
1422 !:Purpose: Relative vorticity and divergence to stream function and
1423 ! velocity potential transformation.
1424 !
1425 implicit none
1426
1427 ! Arguments:
1428 type(struct_gsv), intent(inout) :: statevector
1429
1430 ! Locals:
1431 integer :: stepIndex
1432 logical, save :: firstTime = .true.
1433 real(8), pointer :: qr_ptr(:,:,:,:), dd_ptr(:,:,:,:)
1434 real(8), pointer :: pp_ptr(:,:,:,:), cc_ptr(:,:,:,:)
1435 type(struct_lst), save :: lst_lapl ! Spectral transform Parameters for Vort/Div -> Psi/Chi
1436 integer :: nlev_M
1437 integer :: nTrunc
1438
1439 nTrunc = max(statevector%hco%ni,statevector%hco%nj)/4 + 1
1440
1441 if (gsv_varExist(statevector,'QR') .and. &
1442 gsv_varExist(statevector,'DD')) then
1443 call gsv_getField(statevector,qr_ptr,'QR')
1444 call gsv_getField(statevector,dd_ptr,'DD')
1445 else
1446 call gsv_getField(statevector,qr_ptr,'UU')
1447 call gsv_getField(statevector,dd_ptr,'VV')
1448 end if
1449
1450 if (gsv_varExist(statevector,'PP') .and. &
1451 gsv_varExist(statevector,'CC')) then
1452 call gsv_getField(statevector,pp_ptr,'PP')
1453 call gsv_getField(statevector,cc_ptr,'CC')
1454 else
1455 call gsv_getField(statevector,pp_ptr,'UU')
1456 call gsv_getField(statevector,cc_ptr,'VV')
1457 end if
1458
1459 nlev_M = gsv_getNumLev (statevector,'MM')
1460
1461 if ( statevector%hco%global ) then
1462 call utl_abort('vortDivToPsiChi_gsv: global mode not available')
1463 else
1464 if (firstTime) then
1465 call lst_Setup(lst_lapl, & ! OUT
1466 statevector%hco%ni, statevector%hco%nj, & ! IN
1467 statevector%hco%dlon, nTrunc, & ! IN
1468 'LatLonMN', nlev_M ) ! IN
1469 firstTime = .false.
1470 end if
1471
1472 do stepIndex = 1, statevector%numStep
1473
1474 pp_ptr(:,:,:,stepIndex) = qr_ptr(:,:,:,stepIndex)
1475 cc_ptr(:,:,:,stepIndex) = dd_ptr(:,:,:,stepIndex)
1476
1477 ! Vort -> Psi
1478 call lst_Laplacian(lst_lapl, & ! IN
1479 pp_ptr(:,:,:,stepIndex), & ! INOUT
1480 'Inverse', nlev_M ) ! IN
1481
1482 ! Div -> Chi
1483 call lst_Laplacian(lst_lapl, & ! IN
1484 cc_ptr(:,:,:,stepIndex), & ! INOUT
1485 'Inverse', nlev_M ) ! IN
1486
1487 end do
1488
1489 end if
1490
1491 end subroutine VortDivToPsiChi_gsv
1492
1493 !--------------------------------------------------------------------------
1494 ! UVtoPsiChi_gsv
1495 !--------------------------------------------------------------------------
1496 subroutine UVtoPsiChi_gsv(statevector)
1497 !
1498 !:Purpose: Wind components to stream function and velocity potential
1499 ! transformation.
1500 !
1501 implicit none
1502
1503 ! Arguments:
1504 type(struct_gsv), intent(inout) :: statevector
1505
1506 ! Locals:
1507 integer :: stepIndex
1508 real(8), pointer :: uu_ptr(:,:,:,:), vv_ptr(:,:,:,:)
1509 real(8), pointer :: psi_ptr(:,:,:,:), chi_ptr(:,:,:,:)
1510 real(8), allocatable :: gridState(:,:,:), spectralState(:,:,:)
1511 real(8) :: dla2
1512 integer :: nlev_M, levIndex
1513 integer :: ila_mpiglobal, ila_mpilocal
1514 ! spectral transform configuration (saved)
1515 integer, save :: gstID = -1
1516 integer, save :: nla_mpilocal, maxMyNla, ntrunc
1517 integer, save :: mymBeg,mymEnd,mymSkip,mymCount
1518 integer, save :: mynBeg,mynEnd,mynSkip,mynCount
1519 integer, save, pointer :: ilaList_mpiglobal(:), ilaList_mpilocal(:)
1520
1521 write(*,*) 'UVtoPsiChi_gsv: starting'
1522 call flush(6)
1523
1524 if ( .not. statevector%hco%global ) then
1525
1526 call UVtoVortDiv_gsv (statevector)
1527 call vortDivToPsiChi_gsv(statevector)
1528
1529 else
1530
1531 call gsv_getField(statevector,uu_ptr,'UU')
1532 call gsv_getField(statevector,vv_ptr,'VV')
1533 if (gsv_varExist(statevector,'PP') .and. &
1534 gsv_varExist(statevector,'CC')) then
1535 call gsv_getField(statevector,psi_ptr,'PP')
1536 call gsv_getField(statevector,chi_ptr,'CC')
1537 else
1538 call gsv_getField(statevector,psi_ptr,'UU')
1539 call gsv_getField(statevector,chi_ptr,'VV')
1540 end if
1541 nlev_M = gsv_getNumLev(statevector,'MM')
1542
1543 if ( gstID < 0 ) then
1544 !ntrunc = statevector%nj
1545 ntrunc = 180
1546 gstID = gst_setup(statevector%ni,statevector%nj,ntrunc,2*nlev_M)
1547 call mmpi_setup_m(ntrunc,mymBeg,mymEnd,mymSkip,mymCount)
1548 call mmpi_setup_n(ntrunc,mynBeg,mynEnd,mynSkip,mynCount)
1549 call gst_ilaList_mpiglobal( ilaList_mpiglobal,nla_mpilocal,maxMyNla,&
1550 gstID,mymBeg,mymEnd,mymSkip,mynBeg,mynEnd,&
1551 mynSkip)
1552 call gst_ilaList_mpilocal(ilaList_mpilocal,gstID,mymBeg,mymEnd,mymSkip,&
1553 mynBeg,mynEnd,mynSkip)
1554 end if
1555
1556 dla2 = ec_ra * ec_ra
1557 allocate(gridState(statevector%lonPerPE,statevector%latPerPE,2*nlev_M))
1558 allocate(spectralState(nla_mpilocal,2,2*nlev_M))
1559
1560 do stepIndex = 1, statevector%numStep
1561
1562 gridState(:,:,1:nlev_M) = uu_ptr(:,:,:,stepIndex)
1563 gridState(:,:,(nlev_M+1):2*nlev_M) = vv_ptr(:,:,:,stepIndex)
1564
1565 call gst_setID(gstID)
1566 call gst_gdsp(spectralState,gridState,nlev_M)
1567
1568 do levIndex = 1, 2*nlev_M
1569 do ila_mpilocal = 1, nla_mpilocal
1570 ila_mpiglobal = ilaList_mpiglobal(ila_mpilocal)
1571 spectralState(ila_mpilocal,:,levIndex) = &
1572 spectralState(ila_mpilocal,:,levIndex) * &
1573 dla2*gst_getR1snp1(ila_mpiglobal)
1574 enddo
1575 enddo
1576
1577 call gst_speree(spectralState,gridState)
1578
1579 psi_ptr(:,:,:,stepIndex) = gridState(:,:,1:nlev_M)
1580 chi_ptr(:,:,:,stepIndex) = gridState(:,:,(nlev_M+1):2*nlev_M)
1581
1582 end do
1583
1584 write(*,*) 'deallocate'; call flush(6)
1585 deallocate(gridState)
1586 deallocate(spectralState)
1587
1588 end if
1589
1590 write(*,*) 'UVtoPsiChi_gsv: finished'
1591 call flush(6)
1592
1593 end subroutine UVtoPsiChi_gsv
1594
1595 !--------------------------------------------------------------------------
1596 ! UVtoPsiChi_ens
1597 !--------------------------------------------------------------------------
1598 subroutine UVtoPsiChi_ens(ens)
1599 !
1600 !:Purpose: Wind components to stream function and velocity potential
1601 ! transformation (ensemble processing).
1602 !
1603 implicit none
1604
1605 ! Arguments:
1606 type(struct_ens), intent(inout) :: ens
1607
1608 ! Locals:
1609 type(struct_hco), pointer :: hco_ens => null()
1610 type(struct_gsv) :: gridStateVector_oneMember
1611 integer :: memberIndex
1612
1613 write(*,*)
1614 write(*,*) 'gvt_UVtoPsiChi_ens: starting'
1615
1616 hco_ens => ens_getHco(ens)
1617
1618 !
1619 !- 1. Create a working stateVector
1620 !
1621 call gsv_allocate(gridStateVector_oneMember, 1, hco_ens, ens_getVco(ens), &
1622 varNames_opt=(/'UU','VV'/), &
1623 datestamp_opt=tim_getDatestamp(), &
1624 mpi_local_opt=.true., dataKind_opt=8)
1625
1626 !
1627 !- 2. Loop on members
1628 !
1629 do memberIndex = 1, ens_getNumMembers(ens)
1630
1631 !- 2.1 Copy to a stateVector
1632 call ens_copyMember(ens, gridStateVector_oneMember, memberIndex)
1633
1634 !- 2.2 Do the transform
1635 call UVtoPsiChi_gsv(gridStateVector_oneMember)
1636
1637 !- 2.3 Put the result back in the input ensembleStateVector
1638 call ens_insertMember(ens, gridStateVector_oneMember, memberIndex)
1639 end do
1640
1641 !
1642 !- 3. Cleaning
1643 !
1644 call gsv_deallocate(gridStateVector_oneMember)
1645
1646 write(*,*) 'gvt_UVtoPsiChi_ens: finished'
1647
1648 end subroutine UVtoPsiChi_ens
1649
1650 !--------------------------------------------------------------------------
1651 ! UVtoVorDiv_ens
1652 !--------------------------------------------------------------------------
1653 subroutine UVtoVortDiv_ens(ens)
1654 !
1655 !:Purpose: Wind components to relative vorticity and divergence transformation
1656 ! (ensemble processing)
1657 !
1658 implicit none
1659
1660 ! Arguments:
1661 type(struct_ens), intent(inout) :: ens
1662
1663 ! Locals:
1664 type(struct_hco), pointer :: hco_ens => null()
1665 type(struct_gsv) :: gridStateVector_oneMember
1666 integer :: memberIndex
1667
1668 write(*,*)
1669 write(*,*) 'gvt_UVtoVortDiv_ens: starting'
1670
1671 hco_ens => ens_getHco(ens)
1672
1673 if (hco_ens%global ) then
1674 call utl_abort('gvt_UVtoVortDiv_ens: global mode not yet available')
1675 end if
1676
1677 !
1678 !- 1. Create a working stateVector
1679 !
1680 call gsv_allocate(gridStateVector_oneMember, 1, hco_ens, ens_getVco(ens), &
1681 varNames_opt=(/'UU','VV'/), &
1682 datestamp_opt=tim_getDatestamp(), &
1683 mpi_local_opt=.true., dataKind_opt=8)
1684
1685 !
1686 !- 2. Loop on members
1687 !
1688 do memberIndex = 1, ens_getNumMembers(ens)
1689
1690 !- 2.1 Copy to a stateVector
1691 call ens_copyMember(ens, gridStateVector_oneMember, memberIndex)
1692
1693 !- 2.2 Do the transform
1694 call UVtoVortDiv_gsv(gridStateVector_oneMember)
1695
1696 !- 2.3 Put the result back in the input ensembleStateVector
1697 call ens_insertMember(ens, gridStateVector_oneMember, memberIndex)
1698 end do
1699
1700 !
1701 !- 3. Cleaning
1702 !
1703 call gsv_deallocate(gridStateVector_oneMember)
1704
1705 write(*,*) 'gvt_UVtoVortDiv_ens: finished'
1706
1707 end subroutine UVtoVortDiv_ens
1708
1709 !--------------------------------------------------------------------------
1710 ! logCH_ens
1711 !--------------------------------------------------------------------------
1712 subroutine logCH_ens(ens,varName)
1713 !
1714 !:Purpose: Logarithmic transformation of a chemical species concentration
1715 ! (ensemble processing)
1716 !
1717 implicit none
1718
1719 ! Arguments:
1720 type(struct_ens), intent(inout) :: ens
1721 character(len=4), intent(in) :: varName
1722
1723 ! Locals:
1724 integer :: lonIndex, latIndex, levIndex, stepIndex, memberIndex
1725 integer :: myLatBeg, myLatEnd, myLonBeg, myLonEnd
1726 character(len=4) :: varName_ens
1727 real(4) :: minVal
1728 real(4), pointer :: ptr4d_r4(:,:,:,:)
1729
1730 call ens_getLatLonBounds(ens, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
1731
1732 do levIndex = 1, ens_getNumK(ens)
1733
1734 varName_ens = ens_getVarNameFromK(ens,levIndex)
1735 if ( trim(varName_ens) /= trim(varName) ) cycle
1736
1737 ptr4d_r4 => ens_getOneLev_r4(ens,levIndex)
1738 minVal=real(gsv_minValVarKindCH(vnl_varListIndex(varName)),4)
1739
1740 do latIndex = myLatBeg, myLatEnd
1741 do lonIndex = myLonBeg, myLonEnd
1742 do stepIndex = 1, ens_getNumStep(ens)
1743 do memberIndex = 1, ens_getNumMembers(ens)
1744 ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex) = &
1745 log(max(ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex),&
1746 minVal))
1747 end do
1748 end do
1749 end do
1750 end do
1751
1752 end do
1753
1754 end subroutine logCH_ens
1755
1756 !--------------------------------------------------------------------------
1757 ! expCH_tlm
1758 !--------------------------------------------------------------------------
1759 subroutine expCH_tlm(statevector, varName, stateVectorRef_opt)
1760 !
1761 !:Purpose: Tangent linear of exponentiation of chemical species
1762 ! concentration.
1763 ! Transform d[log(x)] to dx where x = 'stateVectorRef_opt',
1764 ! the input 'statevector' component is d[log(x)] and
1765 ! the output 'statevector' component is dx.
1766 !
1767 implicit none
1768
1769 ! Arguments:
1770 type(struct_gsv), intent(inout) :: statevector
1771 character(len=*), intent(in) :: varName
1772 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
1773
1774 ! Locals:
1775 integer :: lonIndex,latIndex,levIndex,stepIndex,varIndex
1776 real(8), pointer :: var_ptr(:,:,:,:), logVar_ptr(:,:,:,:)
1777 real(8), pointer :: var_trial(:,:,:,:)
1778 real(8) :: minVal
1779
1780 if ( present(statevectorRef_opt) ) then
1781 call gsv_getField(stateVectorRef_opt,var_trial,trim(varName))
1782 else
1783 varIndex = vnl_varListIndex(varName)
1784 if ( .not. varKindCHTrialsInitialized(varIndex) ) then
1785 call utl_abort('expCH_tlm: do trials to stateVectorRefChem transform at higher level')
1786 end if
1787 call gsv_getField(stateVectorTrialvarKindCH(varIndex),var_trial,&
1788 trim(varName))
1789 end if
1790
1791 call gsv_getField(statevector,var_ptr,trim(varName))
1792 call gsv_getField(statevector,logVar_ptr,trim(varName))
1793
1794 minVal=gsv_minValVarKindCH(vnl_varListIndex(varName))
1795
1796 do stepIndex = 1, statevector%numStep
1797 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1798 do levIndex = 1, gsv_getNumLev( statevector,&
1799 vnl_varLevelFromVarname(trim(varName)))
1800 do latIndex = statevector%myLatBeg, statevector%myLatEnd
1801 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
1802 var_ptr(lonIndex,latIndex,levIndex,stepIndex) = &
1803 logVar_ptr(lonIndex,latIndex,levIndex,stepIndex) &
1804 * max(var_trial(lonIndex,latIndex,levIndex,stepIndex),minVal)
1805 end do
1806 end do
1807 end do
1808 !$OMP END PARALLEL DO
1809 end do
1810
1811 end subroutine expCH_tlm
1812
1813 !--------------------------------------------------------------------------
1814 ! CH_bounds
1815 !--------------------------------------------------------------------------
1816 subroutine CH_bounds(statevector)
1817 !
1818 !:Purpose: Impose boundary values to variables of CH kind.
1819 !
1820 implicit none
1821
1822 ! Arguments:
1823 type(struct_gsv), intent(inout) :: statevector
1824
1825 ! Locals:
1826 integer :: varIndex,lonIndex,latIndex,levIndex,stepIndex
1827 real(8), pointer :: var_ptr(:,:,:,:)
1828 real(8) :: minVal
1829 character(len=4) :: varName
1830
1831 do varIndex = 1, vnl_numvarmax
1832 varName = vnl_varNameList(varIndex)
1833 if ( .not.gsv_varExist(varName=trim(varName)) ) cycle
1834 if ( vnl_varKindFromVarname(trim(varName)) /= 'CH' ) cycle
1835
1836 call gsv_getField(statevector,var_ptr,trim(varName))
1837
1838 minVal=gsv_minValVarKindCH(vnl_varListIndex(varName))
1839
1840 do stepIndex = 1, statevector%numStep
1841 !$OMP PARALLEL DO PRIVATE(lonIndex,latIndex,levIndex)
1842 do levIndex = 1, gsv_getNumLev( statevector,&
1843 vnl_varLevelFromVarname(trim(varName)))
1844 do latIndex = statevector%myLatBeg, statevector%myLatEnd
1845 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
1846 var_ptr(lonIndex,latIndex,levIndex,stepIndex) = &
1847 max(var_ptr(lonIndex,latIndex,levIndex,stepIndex),minVal)
1848 end do
1849 end do
1850 end do
1851 !$OMP END PARALLEL DO
1852 end do
1853
1854 end do
1855
1856 end subroutine CH_bounds
1857
1858 !--------------------------------------------------------------------------
1859 ! gvt_oceanIceContinuous
1860 !--------------------------------------------------------------------------
1861 subroutine gvt_oceanIceContinuous(stateVector, stateVectorRef, outputVarName)
1862 !
1863 !:Purpose: Solve laplaces equation at a subset of gridpoints
1864 ! subject to the boundary conditions imposed by the
1865 ! surrounding points. Uses the method of sequential
1866 ! relaxation or Liebmann relaxation, which converges
1867 ! more rapidly than the simultaneous relaxation (see
1868 ! numerical weather analysis and prediction by P. D.
1869 ! Thompson, 1961, pp92-98). NOTE: this subroutine
1870 ! currently uses the oceanMask only for the first level.
1871 ! Therefore, if it is applied to 3D masked fields that
1872 ! have level-dependent masks, the code will abort so
1873 ! that the user can make the necessary adjustments to
1874 ! ensure the correct behaviour.
1875 !
1876 implicit none
1877
1878 ! Arguments:
1879 type(struct_gsv), intent(inout) :: statevector
1880 type(struct_gsv), intent(in) :: stateVectorRef
1881 character(len=*), intent(in) :: outputVarName
1882
1883 ! Locals:
1884 type(struct_gsv) :: statevector_analysis_1step_r8
1885 type(struct_gsv) :: statevector_trial_1step_r8
1886 real(8), pointer :: analysis_ptr(:,:,:,:), trial_ptr(:,:,:,:)
1887 real(8), pointer :: input_ptr(:,:,:,:)
1888 real(8) :: alpha, factor, correc, rms, maxAbsCorr, basic
1889 integer :: lonIndex, latIndex, levIndex, stepIndex
1890 integer :: ipass, numPass, numCorrect
1891 logical :: orca12
1892 character(len=2) :: inputVarName
1893
1894 ! abort if 3D mask is present, since we may not handle this situation correctly
1895 if (stateVector%oceanMask%nLev > 1) then
1896 call utl_abort('gvt_oceanIceContinuous: 3D mask present - this case not properly handled')
1897 end if
1898
1899 ! allocate statevector for single time steps
1900 if (mmpi_myid < stateVector%numStep) then
1901 if (outputVarName == 'LG') then
1902 call gsv_allocate(stateVector_analysis_1step_r8, 1, stateVector%hco, &
1903 stateVector%vco, mpi_local_opt = .false., &
1904 dataKind_opt = 8, varNames_opt = (/'GL','LG'/))
1905 call gsv_allocate(stateVector_trial_1step_r8, 1, stateVectorRef%hco, &
1906 stateVectorRef%vco, mpi_local_opt=.false., &
1907 dataKind_opt=8, varNames_opt = (/'GL','LG'/))
1908 else if(outputVarName == 'TM') then
1909 call gsv_allocate(stateVector_analysis_1step_r8, 1, stateVector%hco, &
1910 stateVector%vco, mpi_local_opt = .false., &
1911 dataKind_opt = 8, varNames_opt = (/outputVarName/))
1912 call gsv_allocate(stateVector_trial_1step_r8, 1, stateVectorRef%hco, &
1913 stateVectorRef%vco, mpi_local_opt = .false., &
1914 dataKind_opt = 8, varNames_opt = (/outputVarName/))
1915 else
1916 call utl_abort('gvt_oceanIceContinuous: unrecognized variable name: '//trim(outputVarName))
1917 end if
1918 end if
1919
1920 call gsv_transposeTilesToStep(stateVector_analysis_1step_r8, stateVector, 1)
1921 call gsv_transposeTilesToStep(stateVector_trial_1step_r8, stateVectorRef, 1)
1922
1923 if (gsv_isAllocated(stateVector_analysis_1step_r8)) then
1924
1925 if (outputVarName == 'LG') then
1926 inputVarName = 'GL'
1927 else
1928 inputVarName = outputVarName
1929 end if
1930 call gsv_getField(stateVector_analysis_1step_r8, input_ptr, inputVarName)
1931 call gsv_getField(stateVector_analysis_1step_r8, analysis_ptr, outputVarName)
1932 call gsv_getField(stateVector_trial_1step_r8, trial_ptr, outputVarName)
1933
1934 alpha = 1.975d0
1935 if( gpos_gridIsOrca(stateVector%ni, stateVector%nj, real(stateVector%hco%lat2d_4,8), real(stateVector%hco%lon2d_4,8)) ) then
1936 orca12 = .true.
1937 numPass = 1000
1938 factor = alpha*0.88d0
1939 else
1940 orca12 = .false.
1941 numPass = 500
1942 factor = alpha
1943 end if
1944
1945 write(*,*) 'gvt_oceanIceContinuous: Liebmann relaxation'
1946 write(*,*) 'gvt_oceanIceContinuous: Number of free points: ', count(.not. stateVector%oceanMask%mask)
1947 write(*,*) 'gvt_oceanIceContinuous: Number of fixed points: ', count( stateVector%oceanMask%mask)
1948 write(*,*) 'gvt_oceanIceContinuous: Total number of grid points: ', stateVector%ni*stateVector%nj
1949 write(*,*) 'gvt_oceanIceContinuous: Total number of iterations: ', numPass
1950
1951 do stepIndex = 1, statevector%numStep
1952 write(*,*) 'gvt_oceanIceContinuous: stepIndex = ', stepIndex
1953 do levIndex = 1, gsv_getNumLev(statevector,vnl_varLevelFromVarname(outputVarName))
1954
1955 write(*,*) 'gvt_oceanIceContinuous: levIndex = ',levIndex
1956
1957 ! Initialisation
1958 do latIndex = 1, stateVector%nj
1959 do lonIndex = 1, stateVector%ni
1960 if (stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
1961 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = input_ptr(lonIndex, latIndex, levIndex, stepIndex)
1962 else
1963 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = trial_ptr(lonIndex, latIndex, levIndex, stepIndex)
1964 end if
1965 end do
1966 end do
1967
1968 do ipass = 1, numPass
1969
1970 rms = 0.0d0
1971 numCorrect = 0
1972 maxAbsCorr = 0.0d0
1973
1974 if(stateVector%hco%grtyp == 'U') then
1975
1976 do latIndex = 2, stateVector%nj/2-1
1977 do lonIndex = 2, stateVector%ni-1
1978 if (.not. stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
1979 basic = (analysis_ptr(lonIndex + 1, latIndex , levIndex, stepIndex) + &
1980 analysis_ptr(lonIndex - 1, latIndex , levIndex, stepIndex) + &
1981 analysis_ptr(lonIndex , latIndex + 1, levIndex, stepIndex) + &
1982 analysis_ptr(lonIndex , latIndex - 1, levIndex, stepIndex)) / 4.0d0
1983 correc = factor * (basic - analysis_ptr(lonIndex, latIndex, levIndex, stepIndex))
1984 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) + correc
1985 rms = rms + correc * correc
1986 numCorrect = numCorrect + 1
1987 if(abs(correc) > maxAbsCorr) maxAbsCorr = abs(correc)
1988 end if
1989 end do
1990 end do
1991
1992 do latIndex = stateVector%nj/2 + 2, stateVector%nj-1
1993 do lonIndex = 2, stateVector%ni-1
1994 if (.not. stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
1995 basic = (analysis_ptr(lonIndex + 1, latIndex , levIndex, stepIndex) + &
1996 analysis_ptr(lonIndex - 1, latIndex , levIndex, stepIndex) + &
1997 analysis_ptr(lonIndex , latIndex + 1, levIndex, stepIndex) + &
1998 analysis_ptr(lonIndex , latIndex - 1, levIndex, stepIndex)) / 4.0d0
1999 correc = factor * (basic - analysis_ptr(lonIndex, latIndex, levIndex, stepIndex))
2000 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) + correc
2001 rms = rms + correc * correc
2002 numCorrect = numCorrect + 1
2003 if(abs(correc) > maxAbsCorr) maxAbsCorr = abs(correc)
2004 end if
2005 end do
2006 end do
2007
2008 else
2009
2010 do latIndex = 2, stateVector%nj-1
2011 do lonIndex = 2, stateVector%ni-1
2012 if (.not. stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
2013 basic = (analysis_ptr(lonIndex + 1, latIndex , levIndex, stepIndex) + &
2014 analysis_ptr(lonIndex - 1, latIndex , levIndex, stepIndex) + &
2015 analysis_ptr(lonIndex , latIndex + 1, levIndex, stepIndex) + &
2016 analysis_ptr(lonIndex , latIndex - 1, levIndex, stepIndex)) / 4.0d0
2017 correc = factor * (basic - analysis_ptr(lonIndex, latIndex, levIndex, stepIndex))
2018 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) + correc
2019 rms = rms + correc * correc
2020 numCorrect = numCorrect + 1
2021 if(abs(correc) > maxAbsCorr) maxAbsCorr = abs(correc)
2022 end if
2023 end do
2024 end do
2025
2026 end if
2027
2028 if( orca12 ) then
2029 ! Periodicity in the X direction
2030 lonIndex = 1
2031 do latIndex = 2, stateVector%nj-1
2032 if (.not. stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
2033 basic = (analysis_ptr(lonIndex + 1, latIndex, levIndex, stepIndex) + &
2034 analysis_ptr(stateVector%ni-2, latIndex, levIndex, stepIndex) + &
2035 analysis_ptr(lonIndex, latIndex + 1, levIndex, stepIndex) + &
2036 analysis_ptr(lonIndex, latIndex - 1, levIndex, stepIndex)) / 4.0d0
2037 correc = factor * (basic - analysis_ptr(lonIndex, latIndex, levIndex, stepIndex))
2038 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) + correc
2039 rms = rms + correc * correc
2040 numCorrect = numCorrect + 1
2041 if(abs(correc) > maxAbsCorr ) maxAbsCorr = abs(correc)
2042 end if
2043 end do
2044 lonIndex = stateVector%ni
2045 do latIndex = 2, stateVector%nj-1
2046 if (.not. stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
2047 basic = (analysis_ptr(3, latIndex, levIndex, stepIndex) + &
2048 analysis_ptr(lonIndex - 1, latIndex, levIndex, stepIndex) + &
2049 analysis_ptr(lonIndex, latIndex + 1, levIndex, stepIndex) + &
2050 analysis_ptr(lonIndex, latIndex - 1, levIndex, stepIndex)) / 4.0d0
2051 correc = factor * (basic - analysis_ptr(lonIndex, latIndex, levIndex, stepIndex))
2052 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) + correc
2053 rms = rms + correc * correc
2054 numCorrect = numCorrect + 1
2055 if( abs(correc) > maxAbsCorr ) maxAbsCorr = abs(correc)
2056 end if
2057 end do
2058 ! North fold
2059 latIndex = stateVector%nj
2060 do lonIndex = 2, stateVector%ni-1
2061 if (.not. stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
2062 basic = (analysis_ptr(lonIndex + 1, latIndex, levIndex, stepIndex)+&
2063 analysis_ptr(lonIndex - 1, latIndex, levIndex, stepIndex)+&
2064 analysis_ptr(stateVector%ni + 2 - lonIndex, latIndex - 2, levIndex, stepIndex) + &
2065 analysis_ptr(stateVector%ni + 2 - lonIndex, latIndex - 1, levIndex, stepIndex) ) / 4.0d0
2066 correc = factor*(basic - analysis_ptr(lonIndex, latIndex, levIndex, stepIndex))
2067 analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) = analysis_ptr(lonIndex, latIndex, levIndex, stepIndex) + correc
2068 rms = rms + correc * correc
2069 numCorrect = numCorrect + 1
2070 if(abs(correc) > maxAbsCorr) maxAbsCorr = abs(correc)
2071 end if
2072 end do
2073 end if
2074
2075 end do
2076
2077 if(numCorrect > 0) rms = sqrt(rms / real(numCorrect))
2078
2079 write(*,*) 'gvt_oceanIceContinuous: number of points corrected = ', numCorrect
2080 write(*,*) 'gvt_oceanIceContinuous: RMS correction during last iteration: ', rms
2081 write(*,*) 'gvt_oceanIceContinuous: MAX absolute correction during last iteration: ', maxAbsCorr
2082 write(*,*) 'gvt_oceanIceContinuous: Field min value = ', minval(analysis_ptr(:,:,levIndex,stepIndex))
2083 write(*,*) 'gvt_oceanIceContinuous: Field max value = ', maxval(analysis_ptr(:,:,levIndex,stepIndex))
2084
2085 if(maxAbsCorr > 1.0) then
2086 call utl_abort('gvt_oceanIceContinuous: Unstable algorithm !')
2087 end if
2088
2089 end do
2090 end do
2091
2092 if (outputVarName == 'LG') then
2093 !- Impose limits [0,1] on sea ice concentration analysis
2094 analysis_ptr(:,:,:,:) = min(analysis_ptr(:,:,:,:), 1.0d0)
2095 analysis_ptr(:,:,:,:) = max(analysis_ptr(:,:,:,:), 0.0d0)
2096 end if
2097
2098 end if
2099
2100 call gsv_transposeStepToTiles(stateVector_analysis_1step_r8, stateVector, 1)
2101
2102 if (mmpi_myid < stateVector%numStep) then
2103 call gsv_deallocate(stateVector_analysis_1step_r8)
2104 call gsv_deallocate(stateVector_trial_1step_r8)
2105 end if
2106
2107 end subroutine gvt_oceanIceContinuous
2108
2109 !--------------------------------------------------------------------------
2110 ! gvt_SSTSpread
2111 !--------------------------------------------------------------------------
2112 subroutine gvt_SSTSpread(stateVector, variableName, maxBoxSize, subgrid)
2113 !
2114 !:Purpose: Spread SST values on neigbouring land surface points
2115 !
2116 implicit none
2117
2118 ! Arguments:
2119 type(struct_gsv), intent(inout) :: stateVector ! state vector of SST analysis
2120 character(len=*), intent(in) :: variableName ! variable name
2121 integer , intent(in) :: maxBoxSize ! maximum box size of SST values spreading
2122 character(len=*), intent(in) :: subgrid ! spread SST values on neighbouring land points of "Yin" or "Yan" subgrid
2123
2124 ! Locals:
2125 logical, allocatable :: isWaterValue(:,:) ! .True. for water points, .False. for land points
2126 logical, allocatable :: updatedIsWaterValue(:,:) ! If the current value is already updated, set it to .True.
2127 real(4), allocatable :: updatedField(:,:) ! updated surface temperature on land
2128 real(4) :: updatedValueSum
2129 integer :: top, bottom, left, right, np
2130 integer :: boxSize, k, l, m
2131 integer :: in(100), jn(100), ngp
2132 integer :: lonIndex, latIndex
2133 integer :: latIndexBeg, latIndexEnd
2134 type(struct_gsv) :: stateVector_1step
2135 real(8), pointer :: field_ptr(:,:,:,:)
2136
2137 write(*,'(a,i2,a)') 'gvt_SSTSpread: spread SST values on ', maxBoxSize,' neighbouring land points on '//trim(subgrid)//' subgrid...'
2138 if (trim(subgrid) == 'Yin') then
2139 latIndexBeg = 1
2140 latIndexEnd = stateVector%hco%nj / 2
2141 else if(trim(subgrid) == 'Yan') then
2142 latIndexBeg = stateVector%hco%nj / 2 + 1
2143 latIndexEnd = stateVector%hco%nj
2144 else
2145 call utl_abort('gvt_SSTSpread: unknown subgrid: '//trim(subgrid))
2146 end if
2147
2148 ! abort if 3D mask is present, since we may not handle this situation correctly
2149 if (stateVector%oceanMask%nLev > 1) then
2150 call utl_abort('gvt_SSTSpread: 3D mask present - this case not properly handled')
2151 end if
2152
2153 ! abort if 3D variable is present
2154 if( gsv_getNumLev(stateVector, vnl_varLevelFromVarname(variableName)) > 1) then
2155 call utl_abort('gvt_SSTSpread: 3D variable present - this case not properly handled')
2156 end if
2157
2158 ! allocate statevector
2159 if (mmpi_myid < stateVector%numStep) then
2160 if(variableName == 'TM') then
2161 call gsv_allocate(stateVector_1step, 1, stateVector%hco, &
2162 stateVector%vco, mpi_local_opt = .false., &
2163 dataKind_opt = 8, varNames_opt = (/variableName/))
2164 else
2165 call utl_abort('gvt_SSTSpread: unrecognized variable name: '//trim(variableName))
2166 end if
2167 end if
2168
2169 allocate(isWaterValue(1:stateVector%hco%ni, latIndexBeg:latIndexEnd))
2170 allocate(updatedIsWaterValue(1:stateVector%hco%ni, latIndexBeg:latIndexEnd))
2171 allocate(updatedField(1:stateVector%hco%ni, latIndexBeg:latIndexEnd))
2172
2173 call gsv_transposeTilesToStep(stateVector_1step, stateVector, 1)
2174
2175 if (gsv_isAllocated(stateVector_1step)) then
2176
2177 call gsv_getField(stateVector_1step, field_ptr, variableName)
2178
2179 ! Initialisation
2180 do latIndex = latIndexBeg, latIndexEnd
2181 do lonIndex = 1, stateVector%hco%ni
2182 if (stateVector%oceanMask%mask(lonIndex, latIndex, 1)) then
2183 isWaterValue(lonIndex, latIndex) = .True.
2184 else
2185 isWaterValue(lonIndex, latIndex) = .False.
2186 end if
2187 updatedField(lonIndex,latIndex) = field_ptr(lonIndex, latIndex, 1, 1)
2188 end do
2189 end do
2190 updatedIsWaterValue(:,:) = isWaterValue(:,:)
2191
2192 boxSizeLoop: do boxSize = 1, maxBoxSize
2193
2194 do latIndex = latIndexBeg, latIndexEnd
2195 do lonIndex = 1, stateVector%hco%ni
2196
2197 if (isWaterValue(lonIndex,latIndex)) cycle
2198
2199 top = latIndex + boxSize
2200 bottom = latIndex - boxSize
2201 left = lonIndex - boxSize
2202 right = lonIndex + boxSize
2203
2204 np = 0
2205 l = bottom
2206 do k = left, right
2207 np = np + 1
2208 in(np) = k
2209 jn(np) = l
2210 end do
2211 k = right
2212 do l = bottom + 1, top
2213 np = np + 1
2214 in(np) = k
2215 jn(np) = l
2216 end do
2217 l = top
2218 do k = right - 1, left, -1
2219 np = np + 1
2220 in(np) = k
2221 jn(np) = l
2222 end do
2223 k = left
2224 do l = top - 1, bottom + 1, -1
2225 np = np + 1
2226 in(np) = k
2227 jn(np) = l
2228 end do
2229
2230 ngp = 0
2231 updatedValueSum = 0.0
2232 do m = 1, np
2233 if (in(m) >= 1 .and. in(m) <= stateVector%hco%ni .and. &
2234 jn(m) >= latIndexBeg .and. jn(m) <= latIndexEnd ) then
2235 if (isWaterValue(in(m), jn(m))) then
2236 updatedValueSum = updatedValueSum + field_ptr(in(m), jn(m), 1, 1)
2237 ngp = ngp + 1
2238 end if
2239 end if
2240 end do
2241
2242 ! mark the grid point as being filled by interpolation on the grid
2243 if (ngp /= 0) then
2244 updatedIsWaterValue(lonIndex, latIndex) = .True.
2245 updatedField(lonIndex, latIndex) = updatedValueSum / real(ngp)
2246 end if
2247
2248 end do
2249 end do
2250
2251 isWaterValue(:,:) = updatedIsWaterValue(:,:)
2252 field_ptr(:, latIndexBeg:latIndexEnd, 1, 1) = updatedField(:,:)
2253
2254 end do boxSizeLoop
2255
2256 write(*,*) 'gvt_SSTSpread: Field min value = ', minval(field_ptr(:, latIndexBeg:latIndexEnd, 1, 1))
2257 write(*,*) 'gvt_SSTSpread: Field max value = ', maxval(field_ptr(:, latIndexBeg:latIndexEnd, 1, 1))
2258
2259 end if
2260
2261 call gsv_transposeStepToTiles(stateVector_1step, stateVector, 1)
2262
2263 ! deallocate local arrays
2264 if (mmpi_myid < stateVector%numStep) then
2265 call gsv_deallocate(stateVector_1step)
2266 end if
2267 deallocate(isWaterValue)
2268 deallocate(updatedIsWaterValue)
2269 deallocate(updatedField)
2270
2271 end subroutine gvt_SSTSpread
2272
2273end module gridVariableTransforms_mod
