1module interpolation_mod
2 ! MODULE interpolation_mod (prefix='int' category='4. Data Object transformations')
3 !
4 !:Purpose: The grid-point state vector interpolation.
5 !
6 use midasMpi_mod
7 use gridstatevector_mod
8 use columnData_mod
9 use calcHeightAndPressure_mod
10 use varNameList_mod
11 use verticalCoord_mod
12 use horizontalCoord_mod
13 use mathPhysConstants_mod
14 use utilities_mod
15 use message_mod
16 use kdTree2_mod
17 implicit none
18 save
19 private
20
21 ! public subroutines and functions
22 public :: int_interp_gsv
23 public :: int_hInterp_gsv
24 public :: int_vInterp_gsv
25 public :: int_tInterp_gsv
26 public :: int_vInterp_col
27 public :: int_hInterpScalar, int_ezgdef, int_cxgaig
28
29 ! module interfaces
30 ! -----------------
31
32 interface int_hInterpScalar
33 module procedure int_hInterpScalar_gsv
34 module procedure int_hInterpScalar_r4_2d
35 module procedure int_hInterpScalar_r8_2d
36 end interface int_hInterpScalar
37
38 interface int_hInterpUV
39 module procedure int_hInterpUV_gsv
40 module procedure int_hInterpUV_r4_2d
41 module procedure int_hInterpUV_r8_2d
42 end interface int_hInterpUV
43
44 ! Namelist variables
45 ! ------------------
46 logical :: vInterpCopyLowestLevel ! overwrite values at lowest level to avoid extrapolation
47 logical :: checkCloudToGridUnassigned ! abort if unmasked points not assigned from cloudToGrid interp
48 integer :: maxBoxSize ! max size used to fill values for cloudToGrid interpolation
49
50contains
51
52 !--------------------------------------------------------------------------
53 ! int_readNml (private subroutine)
54 !--------------------------------------------------------------------------
55 subroutine int_readNml()
56 !
57 !:Purpose: Read the namelist block NAMINT.
58 !
59 !:Namelist parameters:
60 ! :vInterpCopyLowestLevel: if true, will overwrite values at the lowest
61 ! levels to avoid extrapolation
62 ! :maxBoxSize: maximum size in terms of grid points used for filling
63 ! in undefined values from neighbours when doing interpolation
64 ! from a cloud of points to a grid
65 ! :checkCloudToGridUnassigned: abort if any unmasked points are not assigned
66 ! after doing cloudToGrid interpolation
67 !
68 implicit none
69
70 ! Locals:
71 integer :: ierr, nulnam, fnom, fclos
72 logical, save :: alreadyRead = .false.
73
74 NAMELIST /NAMINT/ vInterpCopyLowestLevel, checkCloudToGridUnassigned, maxBoxSize
75
76 if (alreadyRead) then
77 return
78 else
79 alreadyRead = .true.
80 end if
81 call msg('int_readNml', 'START', verb_opt=2)
82
83 ! default values
84 vInterpCopyLowestLevel = .false.
85 checkCloudToGridUnassigned = .true.
86 maxBoxSize = 1
87
88 if ( .not. utl_isNamelistPresent('NAMINT','./flnml') ) then
89 call msg('int_readNml', 'namint is missing in the namelist.'&
90 //'The default values will be taken.', mpiAll_opt=.false.)
91 else
92 ! Read namelist NAMINT
93 nulnam = 0
94 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
95 read(nulnam,nml=namint,iostat=ierr)
96 if (ierr /= 0) call utl_abort('int_readlNml: Error reading namelist NAMINT')
97 if (mmpi_myid == 0) write(*,nml=namint)
98 ierr = fclos(nulnam)
99 end if
100
101 call msg('int_readNml', 'END', verb_opt=2)
102 end subroutine int_readNml
103
104 !--------------------------------------------------------------------------
105 ! int_interp_gsv
106 !--------------------------------------------------------------------------
107 subroutine int_interp_gsv(statevector_in, statevector_out, statevectorRef_opt, &
108 checkModelTop_opt)
109 !
110 !:Purpose: high-level interpolation subroutine that proceed with
111 ! horizontal and vertical interpolation
112 !
113 implicit none
114
115 ! Arguments:
116 type(struct_gsv), intent(in) :: statevector_in ! Statevector input
117 type(struct_gsv), intent(inout) :: statevector_out ! Statevector output (with target grids)
118 type(struct_gsv), optional, intent(in) :: statevectorRef_opt ! Reference statevector with fields P0, TT and HU
119 logical, optional, intent(in) :: checkModelTop_opt ! If true, model top consistency checked
120
121 ! Locals:
122 logical :: checkModelTop
123 character(len=4), pointer :: varNamesToInterpolate(:)
124 type(struct_gsv) :: statevector_in_varsLevs, statevector_in_varsLevs_hInterp
125 type(struct_gsv) :: statevector_in_hInterp
126
127 call msg('int_interp_gsv', 'START', verb_opt=2)
128
129 !
130 !- Error traps
131 !
132 if (.not.gsv_isAllocated(statevector_in)) then
133 call utl_abort('int_interp_gsv: gridStateVector_in not yet allocated')
134 end if
135 if (.not.gsv_isAllocated(statevector_out)) then
136 call utl_abort('int_interp_gsv: gridStateVector_out not yet allocated')
137 end if
138
139 !
140 !- Do the interpolation of statevector_in onto the grid of statevector_out
141 !
142 nullify(varNamesToInterpolate)
143 call gsv_varNamesList(varNamesToInterpolate, statevector_in)
144
145 !- Horizontal interpolation
146 call gsv_allocate(statevector_in_VarsLevs, statevector_in%numstep, &
147 statevector_in%hco, statevector_in%vco, &
148 mpi_local_opt=statevector_in%mpi_local, mpi_distribution_opt='VarsLevs', &
149 dataKind_opt=gsv_getDataKind(statevector_in), &
150 allocHeightSfc_opt=statevector_in%heightSfcPresent, &
151 varNames_opt=varNamesToInterpolate, &
152 hInterpolateDegree_opt=statevector_out%hInterpolateDegree, &
153 hExtrapolateDegree_opt=statevector_out%hExtrapolateDegree )
154
155 call gsv_transposeTilesToVarsLevs( statevector_in, statevector_in_VarsLevs )
156
157 call gsv_allocate(statevector_in_VarsLevs_hInterp, statevector_in%numstep, &
158 statevector_out%hco, statevector_in%vco, &
159 mpi_local_opt=statevector_out%mpi_local, mpi_distribution_opt='VarsLevs', &
160 dataKind_opt=gsv_getDataKind(statevector_out), &
161 allocHeightSfc_opt=statevector_out%heightSfcPresent, &
162 varNames_opt=varNamesToInterpolate, &
163 hInterpolateDegree_opt=statevector_out%hInterpolateDegree, &
164 hExtrapolateDegree_opt=statevector_out%hExtrapolateDegree )
165
166 call int_hInterp_gsv(statevector_in_VarsLevs, statevector_in_VarsLevs_hInterp)
167 call gsv_deallocate(statevector_in_VarsLevs)
168
169 call gsv_allocate(statevector_in_hInterp, statevector_in%numstep, &
170 statevector_out%hco, statevector_in%vco, &
171 mpi_local_opt=statevector_out%mpi_local, mpi_distribution_opt='Tiles', &
172 dataKind_opt=gsv_getDataKind(statevector_out), &
173 allocHeightSfc_opt=statevector_out%heightSfcPresent, &
174 varNames_opt=varNamesToInterpolate )
175
176 call gsv_transposeVarsLevsToTiles( statevector_in_varsLevs_hInterp, statevector_in_hInterp )
177 call gsv_deallocate(statevector_in_varsLevs_hInterp)
178
179 !- Vertical interpolation
180
181 ! the default is to ensure that the top of the output grid is ~equal or lower than the top of the input grid
182 if ( present(checkModelTop_opt) ) then
183 checkModelTop = checkModelTop_opt
184 else
185 checkModelTop = .true.
186 end if
187
188 call int_vInterp_gsv(statevector_in_hInterp,statevector_out,&
189 statevectorRef_opt=statevectorRef_opt, &
190 checkModelTop_opt=checkModelTop)
191
192 call gsv_deallocate(statevector_in_hInterp)
193 nullify(varNamesToInterpolate)
194
195 call msg('int_interp_gsv', 'END', verb_opt=2)
196 end subroutine int_interp_gsv
197
198 !--------------------------------------------------------------------------
199 ! int_hInterp_gsv
200 !--------------------------------------------------------------------------
201 subroutine int_hInterp_gsv(statevector_in,statevector_out)
202 !
203 !:Purpose: Horizontal interpolation
204 !
205 implicit none
206
207 ! Arguments:
208 type(struct_gsv), intent(inout) :: statevector_in ! Statevector input
209 type(struct_gsv), intent(inout) :: statevector_out ! Statevector with target horiz and vert grids and result
210
211 ! Locals:
212 integer :: varIndex, levIndex, nlev, stepIndex, ierr, kIndex
213 character(len=4) :: varName
214 character(len=12):: interpolationDegree, extrapolationDegree
215
216 call msg('int_hInterp_gsv', 'START', verb_opt=2)
217
218 if ( hco_equal(statevector_in%hco,statevector_out%hco) ) then
219 call msg('int_hInterp_gsv', 'The input and output statevectors are already on same horizontal grids')
220 call gsv_copy(statevector_in, statevector_out)
221 return
222 end if
223
224 if ( .not. vco_equal(statevector_in%vco, statevector_out%vco) ) then
225 call utl_abort('int_hInterp_gsv: The input and output statevectors are not on the same vertical levels.')
226 end if
227
228 ! set the interpolation degree
229 interpolationDegree = statevector_out%hInterpolateDegree
230 extrapolationDegree = statevector_out%hExtrapolateDegree
231
232 if ( .not.statevector_in%mpi_local .and. .not.statevector_out%mpi_local ) then
233
234 call msg('int_hInterp_gsv', 'before interpolation (no mpi)')
235
236 step_loop: do stepIndex = 1, statevector_out%numStep
237 ! copy over some time related parameters
238 statevector_out%deet = statevector_in%deet
239 statevector_out%dateOriginList(stepIndex) = statevector_in%dateOriginList(stepIndex)
240 statevector_out%npasList(stepIndex) = statevector_in%npasList(stepIndex)
241 statevector_out%ip2List(stepIndex) = statevector_in%ip2List(stepIndex)
242 statevector_out%etiket = statevector_in%etiket
243
244 ! Do horizontal interpolation for mpi global statevectors
245 var_loop: do varIndex = 1, vnl_numvarmax
246 varName = vnl_varNameList(varIndex)
247 if ( .not. gsv_varExist(statevector_in,varName) ) cycle var_loop
248 if ( trim(varName) == 'VV' ) cycle var_loop
249
250 nlev = statevector_out%varNumLev(varIndex)
251
252 ! horizontal interpolation
253
254 if ( trim(varName) == 'UU' ) then
255 ! interpolate both UV components and keep both UU and VV
256 do levIndex = 1, nlev
257 ierr = int_hInterpUV( statevector_out, statevector_in, 'BOTH', levIndex, stepIndex, &
258 interpDegree=trim(interpolationDegree), &
259 extrapDegree_opt=trim(extrapolationDegree) )
260 end do
261 else
262 ! interpolate scalar variable
263 do levIndex = 1, nlev
264 ierr = int_hInterpScalar( statevector_out, statevector_in, varName, levIndex, stepIndex, &
265 interpDegree=trim(interpolationDegree), &
266 extrapDegree_opt=trim(extrapolationDegree) )
267 end do
268 end if
269 end do var_loop
270
271 end do step_loop
272
273 else
274
275 call msg('int_hInterp_gsv', 'before interpolation (with mpi)')
276
277 if ( statevector_in%mpi_distribution /= 'VarsLevs' .or. &
278 statevector_out%mpi_distribution /= 'VarsLevs' ) then
279 call utl_abort('int_hInterp_gsv: The input or output statevector is not distributed by VarsLevs.')
280 end if
281
282 do stepIndex = 1, statevector_out%numStep
283 k_loop: do kIndex = statevector_in%mykBeg, statevector_in%mykEnd
284 varName = gsv_getVarNameFromK(statevector_in,kIndex)
285 if ( .not. gsv_varExist(statevector_in,varName) ) cycle k_loop
286
287 ! horizontal interpolation
288
289 if ( trim(varName) == 'UU' ) then
290 ! interpolate both UV components and keep UU in main vector
291 ierr = int_hInterpUV( statevector_out, statevector_in, 'UU', kIndex, stepIndex, &
292 interpDegree=trim(interpolationDegree), &
293 extrapDegree_opt=trim(extrapolationDegree) )
294 else if ( trim(varName) == 'VV' ) then
295 ! interpolate both UV components and keep VV in main vector
296 ierr = int_hInterpUV( statevector_out, statevector_in, 'VV', kIndex, stepIndex, &
297 interpDegree=trim(interpolationDegree), &
298 extrapDegree_opt=trim(extrapolationDegree) )
299 else
300 ! interpolate scalar variable
301 ierr = int_hInterpScalar( statevector_out, statevector_in, 'ALL', kIndex, stepIndex, &
302 interpDegree=trim(interpolationDegree), &
303 extrapDegree_opt=trim(extrapolationDegree) )
304 end if
305 end do k_loop
306
307 end do ! stepIndex
308
309 end if
310
311 if ( gsv_isAssocHeightSfc(statevector_in) .and. gsv_isAssocHeightSfc(statevector_out) ) then
312 call msg('int_hInterp_gsv','interpolating surface height')
313 ierr = int_hInterpScalar( statevector_out, statevector_in, 'ZSFC', 1, 1, &
314 interpDegree=trim(interpolationDegree), &
315 extrapDegree_opt=trim(extrapolationDegree) )
316 end if
317
318 call msg('int_hInterp_gsv', 'END', verb_opt=2)
319 end subroutine int_hInterp_gsv
320
321 !--------------------------------------------------------------------------
322 ! int_vInterp_gsv
323 !--------------------------------------------------------------------------
324 subroutine int_vInterp_gsv(statevector_in,statevector_out,statevectorRef_opt, &
325 Ps_in_hPa_opt,checkModelTop_opt)
326 !
327 !:Purpose: Vertical interpolation.
328 ! Interpolation coordinates are either height or log(P)
329 ! based on target vertical descriptor vcode.
330 ! * When target vertical coordinates are pressure (interpolating
331 ! to GEM-P), log(P) is used;
332 ! * When they are height (interpolating to GEM-H), height is used.
333 ! Call to ``calcHeightAndPressure_mod`` will return required
334 ! interpolation coordinates (log is then applied when required).
335 !
336 implicit none
337
338 ! Arguments:
339 type(struct_gsv), target, intent(in) :: statevector_in ! Statevector input
340 type(struct_gsv), intent(inout) :: statevector_out ! Statevector with target horiz/vert grids and result
341 type(struct_gsv), optional, target, intent(in) :: statevectorRef_opt ! Reference statevector with P0, TT and HU
342 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, surface pressure is in hPa, not Pa
343 logical, optional, intent(in) :: checkModelTop_opt ! Model top consistency will be checked
344
345 call msg('int_vInterp_gsv', 'START', verb_opt=2)
346
347 ! read the namelist
348 call int_readNml()
349
350 if ( gsv_getDataKind(statevector_in) == 8 &
351 .and. gsv_getDataKind(statevector_out) == 8 ) then
352 call vInterp_gsv_r8(statevector_in,statevector_out,statevectorRef_opt, &
353 Ps_in_hPa_opt,checkModelTop_opt)
354 else if ( gsv_getDataKind(statevector_in) == 4 &
355 .and. gsv_getDataKind(statevector_out) == 4 ) then
356 call vInterp_gsv_r4(statevector_in,statevector_out,statevectorRef_opt, &
357 Ps_in_hPa_opt,checkModelTop_opt)
358 else
359 call utl_abort('int_vInterp_gsv: input and output statevectors must be of same dataKind')
360 end if
361
362 call msg('int_vInterp_gsv', 'END', verb_opt=2)
363
364 end subroutine int_vInterp_gsv
365
366 !--------------------------------------------------------------------------
367 ! vInterp_gsv_r8
368 !--------------------------------------------------------------------------
369 subroutine vInterp_gsv_r8(statevector_in,statevector_out,statevectorRef_opt, &
370 Ps_in_hPa_opt,checkModelTop_opt)
371 !
372 !:Purpose: Vertical interpolation, ``real(8)`` version.
373 !
374 implicit none
375
376 ! Arguments:
377 type(struct_gsv), target, intent(in) :: statevector_in ! Statevector input
378 type(struct_gsv), intent(inout) :: statevector_out ! Statevector with target horiz/vert grids and result
379 type(struct_gsv), optional, target, intent(in) :: statevectorRef_opt ! Reference statevector with P0, TT and HU
380 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, surface pressure is in hPa, not Pa
381 logical, optional, intent(in) :: checkModelTop_opt ! Model top consistency will be checked
382
383 ! Locals:
384 logical :: checkModelTop, hLikeCalc
385 integer :: vcode_in, vcode_out
386 integer :: nlev_out, nlev_in
387 integer :: varIndex, stepIndex
388 type(struct_gsv), pointer :: statevectorRef
389 type(struct_gsv) :: statevectorRef_out
390 real(8), pointer :: hLikeT_in(:,:,:,:), hLikeM_in(:,:,:,:) ! abstract height dimensioned coordinate
391 real(8), pointer :: hLikeT_out(:,:,:,:), hLikeM_out(:,:,:,:) ! abstract height dimensioned coordinate
392 real(8), pointer :: field_in(:,:,:,:), field_out(:,:,:,:)
393 real(8), pointer :: heightSfcIn(:,:), heightSfcOut(:,:)
394 real(8), pointer :: tmpCoord_T(:,:,:,:), tmpCoord_M(:,:,:,:)
395 character(len=4) :: varName
396 type(struct_vco), pointer :: vco_in, vco_out
397
398 call msg('vInterp_gsv_r8', 'START', verb_opt=4)
399
400 vco_in => gsv_getVco(statevector_in)
401 vco_out => gsv_getVco(statevector_out)
402 vcode_in = vco_in%vcode
403 vcode_out = vco_out%vcode
404 nullify(hLikeT_in,hLikeM_in,hLikeT_out,hLikeM_out)
405
406 hLikeCalc = .true.
407
408 if (present(statevectorRef_opt)) then
409 if ( .not. vco_equal(gsv_getVco(statevectorRef_opt), gsv_getVco(statevector_in))) then
410 call utl_abort('vInterp_gsv_r8: reference must have input vertical structure')
411 end if
412 statevectorRef => statevectorRef_opt
413 else
414 statevectorRef => statevector_in
415 end if
416
417 if ( vco_equal(statevector_in%vco, statevector_out%vco) ) then
418 call msg('vInterp_gsv_r8', 'The input and output statevectors are already on same vertical levels')
419 call gsv_copy(statevector_in, statevector_out)
420 return
421 end if
422
423 if ( .not. hco_equal(statevector_in%hco, statevector_out%hco) ) then
424 call utl_abort('vInterp_gsv_r8: The input and output statevectors are not on the same horizontal grid.')
425 end if
426
427 if ( gsv_getDataKind(statevector_in) /= 8 .or. gsv_getDataKind(statevector_out) /= 8 ) then
428 call utl_abort('vInterp_gsv_r8: Incorrect value for dataKind. Only compatible with dataKind=8')
429 end if
430
431 if (gsv_isAssocHeightSfc(statevector_in) .and. gsv_isAssocHeightSfc(statevector_out) ) then
432 heightSfcIn => gsv_getHeightSfc(statevector_in)
433 heightSfcOut => gsv_getHeightSfc(statevector_out)
434 heightSfcOut(:,:) = heightSfcIn(:,:)
435 end if
436
437 ! the default is to ensure that the top of the output grid is ~equal or lower than the top of the input grid
438 if ( present(checkModelTop_opt) ) then
439 checkModelTop = checkModelTop_opt
440 else
441 checkModelTop = .true.
442 end if
443 if (checkModelTop) then
444 call msg('vInterp_gsv_r8', ' Checking that that the top of the destination grid is not higher than the top of the source grid.')
445 if ( vcode_in == 21001 .or. vcode_out == 21001 ) then
446 call msg('vInterp_gsv_r8', 'bypassing top check, '&
447 //'vcode_in='//str(vcode_in)//', vcode_out='//str(vcode_out))
448 ! Development notes (@mad001)
449 ! we should consider having a new criterion that works for GEM-H as well
450 else
451 call czp_ensureCompatibleTops(vco_in, vco_out)
452 end if
453 end if
454
455 ! create reference state on output vco with reference surface fields
456 call gsv_allocate(statevectorRef_out, statevectorRef%numstep, &
457 statevectorRef%hco, statevector_out%vco, &
458 mpi_local_opt=statevectorRef%mpi_local, mpi_distribution_opt='Tiles', &
459 dataKind_opt=gsv_getDataKind(statevectorRef), &
460 allocHeightSfc_opt=statevectorRef%heightSfcPresent, &
461 varNames_opt=(/'P0','P0LS'/) )
462 call gsv_copy(stateVectorRef, stateVectorRef_out, allowVcoMismatch_opt=.true., &
463 allowVarMismatch_opt=.true.)
464
465 var_loop: do varIndex = 1, vnl_numvarmax
466 varName = vnl_varNameList(varIndex)
467 if ( .not. gsv_varExist(statevector_in,varName) ) cycle var_loop
468
469 nlev_in = statevector_in%varNumLev(varIndex)
470 nlev_out = statevector_out%varNumLev(varIndex)
471
472 call gsv_getField(statevector_in ,field_in,varName)
473 call gsv_getField(statevector_out,field_out,varName)
474
475 ! for 2D fields and variables on "other" levels, just copy and cycle to next variable
476 if ( (nlev_in == 1 .and. nlev_out == 1) .or. &
477 (vnl_varLevelFromVarname(varName) == 'OT') ) then
478 field_out(:,:,:,:) = field_in(:,:,:,:)
479 cycle var_loop
480 end if
481
482 hLikeCalcIf: if (hLikeCalc) then
483 hLikeCalc = .false.
484
485 ! allocating for temporary pressure references
486 allocate(hLikeT_in( statevector_in%myLonBeg:statevector_in%myLonEnd, &
487 statevector_in%myLatBeg:statevector_in%myLatEnd, &
488 gsv_getNumLev(statevector_in,'TH'), statevector_in%numStep))
489 allocate(hLikeM_in( statevector_in%myLonBeg:statevector_in%myLonEnd, &
490 statevector_in%myLatBeg:statevector_in%myLatEnd, &
491 gsv_getNumLev(statevector_in,'MM'), statevector_in%numStep))
492 allocate(hLikeT_out(statevector_out%myLonBeg:statevector_out%myLonEnd, &
493 statevector_out%myLatBeg:statevector_out%myLatEnd, &
494 gsv_getNumLev(statevector_out,'TH'), statevector_out%numStep))
495 allocate(hLikeM_out(statevector_out%myLonBeg:statevector_out%myLonEnd, &
496 statevector_out%myLatBeg:statevector_out%myLatEnd, &
497 gsv_getNumLev(statevector_out,'MM'), statevector_out%numStep))
498 allocate(tmpCoord_T(statevector_in%myLonBeg:statevector_in%myLonEnd, &
499 statevector_in%myLatBeg:statevector_in%myLatEnd, &
500 gsv_getNumLev(statevector_in,'TH'), statevector_in%numStep))
501 allocate(tmpCoord_M(statevector_in%myLonBeg:statevector_in%myLonEnd, &
502 statevector_in%myLatBeg:statevector_in%myLatEnd, &
503 gsv_getNumLev(statevector_in,'MM'), statevector_in%numStep))
504
505 ! output grid GEM-P interpolation in log-pressure
506 if ( vcode_out==5002 .or. vcode_out==5005 .or. vcode_out==5100 ) then
507 call czp_calcReturnPressure_gsv_nl( statevectorRef_out, &
508 PTout_r8_opt=hLikeT_out, &
509 PMout_r8_opt=hLikeM_out, &
510 Ps_in_hPa_opt=Ps_in_hPa_opt)
511
512 if ( vcode_in==5002 .or. vcode_in==5005 .or. vcode_in==5100 ) then
513 call czp_calcReturnPressure_gsv_nl( statevectorRef, &
514 PTout_r8_opt=hLikeT_in, &
515 PMout_r8_opt=hLikeM_in, &
516 Ps_in_hPa_opt=Ps_in_hPa_opt)
517 else if ( vcode_in==21001 ) then
518 call czp_calcReturnHeight_gsv_nl( statevectorRef, &
519 ZTout_r8_opt=tmpCoord_T, &
520 ZMout_r8_opt=tmpCoord_M)
521 call czp_calcReturnPressure_gsv_nl( statevectorRef, &
522 ZTin_r8_opt=tmpCoord_T, &
523 ZMin_r8_opt=tmpCoord_M, &
524 PTout_r8_opt=hLikeT_in, &
525 PMout_r8_opt=hLikeM_in, &
526 Ps_in_hPa_opt=Ps_in_hPa_opt)
527 end if
528
529 call msg('vInterp_gsv_r8','converting pressure coordinates to height-like, '&
530 //'vcode_in='//str(vcode_in)//', vcode_out='//str(vcode_out))
531 call logP_r8(hLikeM_out)
532 call logP_r8(hLikeT_out)
533 call logP_r8(hLikeM_in)
534 call logP_r8(hLikeT_in)
535
536 ! output grid GEM-H interpolation in height
537 else if ( vcode_out==21001 ) then
538 call czp_calcReturnHeight_gsv_nl( statevectorRef_out, &
539 ZTout_r8_opt=hLikeT_out, &
540 ZMout_r8_opt=hLikeM_out)
541 if ( vcode_in==21001 ) then
542 call czp_calcReturnHeight_gsv_nl( statevectorRef, &
543 ZTout_r8_opt=hLikeT_in, &
544 ZMout_r8_opt=hLikeM_in)
545 else if ( vcode_in==5002 .or. vcode_in==5005 .or. vcode_in==5100 ) then
546 call czp_calcReturnPressure_gsv_nl( statevectorRef, &
547 PTout_r8_opt=tmpCoord_T, &
548 PMout_r8_opt=tmpCoord_M, &
549 Ps_in_hPa_opt=Ps_in_hPa_opt)
550 call czp_calcReturnHeight_gsv_nl( statevectorRef, &
551 PTin_r8_opt=tmpCoord_T, &
552 PMin_r8_opt=tmpCoord_M, &
553 ZTout_r8_opt=hLikeT_in, &
554 ZMout_r8_opt=hLikeM_in)
555 end if
556 end if
557 deallocate(tmpCoord_T)
558 deallocate(tmpCoord_M)
559
560 end if hLikeCalcIf
561
562 step_loop: do stepIndex = 1, statevector_out%numStep
563
564 ! copy over some time related and other parameters
565 statevector_out%deet = statevector_in%deet
566 statevector_out%dateOriginList(stepIndex) = statevector_in%dateOriginList(stepIndex)
567 statevector_out%npasList(stepIndex) = statevector_in%npasList(stepIndex)
568 statevector_out%ip2List(stepIndex) = statevector_in%ip2List(stepIndex)
569 statevector_out%etiket = statevector_in%etiket
570 statevector_out%onPhysicsGrid(:) = statevector_in%onPhysicsGrid(:)
571 statevector_out%hco_physics => statevector_in%hco_physics
572
573 ! do the vertical interpolation
574 field_out(:,:,:,stepIndex) = 0.0d0
575 if (vnl_varLevelFromVarname(varName) == 'TH') then
576 call hLike_interpolation_r8(hLikeT_in, hLikeT_out)
577 else
578 call hLike_interpolation_r8(hLikeM_in, hLikeM_out)
579 end if
580
581 ! overwrite values at the lowest levels to avoid extrapolation
582 if (vInterpCopyLowestLevel) then
583 field_out(:,:,nlev_out,stepIndex) = field_in(:,:,nlev_in,stepIndex)
584 end if
585
586 end do step_loop
587
588 end do var_loop
589
590 if (associated(hLikeT_in)) then
591 deallocate(hLikeT_in, hLikeM_in, hLikeT_out, hLikeM_out)
592 end if
593 call gsv_deallocate(statevectorRef_out)
594
595 call msg('vInterp_gsv_r8', 'END', verb_opt=4)
596
597 contains
598
599 subroutine hLike_interpolation_r8(hLike_in, hLike_out)
600 !
601 !:Purpose: Proceed to actual interpolation in H-logP representation
602 !
603 implicit none
604
605 ! Arguments:
606 real(8), pointer, intent(in) :: hLike_in(:,:,:,:) ! abstract height dimensioned input coordinate
607 real(8), pointer, intent(in) :: hLike_out(:,:,:,:) ! abstract height dimensioned target coordinate
608
609 ! Locals:
610 integer :: latIndex, lonIndex, levIndex_out, levIndex_in
611 real(8) :: zwb, zwt
612
613 !$OMP PARALLEL DO PRIVATE(latIndex,lonIndex,levIndex_in,levIndex_out,zwb,zwt)
614 do latIndex = statevector_out%myLatBeg, statevector_out%myLatEnd
615 do lonIndex = statevector_out%myLonBeg, statevector_out%myLonEnd
616 levIndex_in = 1
617 do levIndex_out = 1, nlev_out
618 levIndex_in = levIndex_in + 1
619 do while(hLike_out(lonIndex,latIndex,levIndex_out,stepIndex) &
620 .gt.hLike_in(lonIndex,latIndex,levIndex_in,stepIndex) &
621 .and.levIndex_in.lt.nlev_in)
622 levIndex_in = levIndex_in + 1
623 end do
624 levIndex_in = levIndex_in - 1
625 zwb = (hLike_out(lonIndex,latIndex,levIndex_out,stepIndex) &
626 - hLike_in(lonIndex,latIndex,levIndex_in,stepIndex)) &
627 /(hLike_in(lonIndex,latIndex,levIndex_in+1,stepIndex) &
628 - hLike_in(lonIndex,latIndex,levIndex_in,stepIndex))
629 zwt = 1.d0 - zwb
630 field_out(lonIndex,latIndex,levIndex_out,stepIndex) = &
631 zwb*field_in(lonIndex,latIndex,levIndex_in+1,stepIndex) &
632 + zwt*field_in(lonIndex,latIndex,levIndex_in,stepIndex)
633 end do
634 end do
635 end do
636 !$OMP END PARALLEL DO
637
638 end subroutine hLike_interpolation_r8
639
640 end subroutine vInterp_gsv_r8
641
642 !--------------------------------------------------------------------------
643 ! logP_r8
644 !--------------------------------------------------------------------------
645 subroutine logP_r8(presInLogOut)
646 !
647 !:Purpose: compute log of pressurce field, real(8) version
648 !
649 implicit none
650
651 ! Arguments:
652 real(8), intent(inout) :: presInLogOut(:,:,:,:)
653
654 ! Locals:
655 integer :: latIdx, lonIdx, levIdx, stepIdx
656
657 !$OMP PARALLEL DO PRIVATE(lonIdx,latIdx,levIdx,stepIdx)
658 do stepIdx = lbound(presInLogOut,4), ubound(presInLogOut,4)
659 do levIdx = lbound(presInLogOut,3), ubound(presInLogOut,3)
660 do latIdx = lbound(presInLogOut,2), ubound(presInLogOut,2)
661 do lonIdx = lbound(presInLogOut,1), ubound(presInLogOut,1)
662 presInLogOut(lonIdx,latIdx,levIdx,stepIdx) = &
663 log(presInLogOut(lonIdx,latIdx,levIdx,stepIdx))
664 end do
665 end do
666 end do
667 end do
668 !$OMP END PARALLEL DO
669
670 end subroutine logP_r8
671
672 !--------------------------------------------------------------------------
673 ! vInterp_gsv_r4
674 !--------------------------------------------------------------------------
675 subroutine vInterp_gsv_r4(statevector_in,statevector_out,statevectorRef_opt, &
676 Ps_in_hPa_opt,checkModelTop_opt)
677 !
678 !:Purpose: Vertical interpolation, ``real(4)`` version.
679 !
680 !
681 implicit none
682
683 ! Arguments:
684 type(struct_gsv), target, intent(in) :: statevector_in ! Statevector input
685 type(struct_gsv), intent(inout) :: statevector_out ! Statevector with the target horiz/vert grids and result
686 type(struct_gsv), optional, target, intent(in) :: statevectorRef_opt ! Reference statevector with P0, TT and HU
687 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, surface pressure in in hPa, not Pa
688 logical, optional, intent(in) :: checkModelTop_opt ! Model top consistency will be checked
689
690 ! Locals:
691 logical :: checkModelTop, hLikeCalc
692 integer :: vcode_in, vcode_out
693 integer :: nlev_out, nlev_in
694 integer :: varIndex, stepIndex
695 type(struct_gsv), pointer :: statevectorRef
696 type(struct_gsv) :: statevectorRef_out
697 real(4), pointer :: hLikeT_in(:,:,:,:), hLikeM_in(:,:,:,:) ! abstract height dimensioned coordinate
698 real(4), pointer :: hLikeT_out(:,:,:,:), hLikeM_out(:,:,:,:) ! abstract height dimensioned coordinate
699 real(4), pointer :: field_in(:,:,:,:), field_out(:,:,:,:)
700 real(8), pointer :: heightSfcIn(:,:), heightSfcOut(:,:)
701 real(4), pointer :: tmpCoord_T(:,:,:,:), tmpCoord_M(:,:,:,:)
702 character(len=4) :: varName
703 type(struct_vco), pointer :: vco_in, vco_out
704
705 call msg('vInterp_gsv_r4', 'START', verb_opt=4)
706
707 vco_in => gsv_getVco(statevector_in)
708 vco_out => gsv_getVco(statevector_out)
709 vcode_in = vco_in%vcode
710 vcode_out = vco_out%vcode
711 nullify(hLikeT_in,hLikeM_in,hLikeT_out,hLikeM_out)
712
713 hLikeCalc = .true.
714
715 if (present(statevectorRef_opt)) then
716 if ( .not. vco_equal(gsv_getVco(statevectorRef_opt), gsv_getVco(statevector_in))) then
717 call utl_abort('vInterp_gsv_r4: reference must have input vertical structure')
718 end if
719 statevectorRef => statevectorRef_opt
720 else
721 statevectorRef => statevector_in
722 end if
723
724 if ( vco_equal(statevector_in%vco, statevector_out%vco) ) then
725 call msg('vInterp_gsv_r4', 'The input and output statevectors are already on same vertical levels')
726 call gsv_copy(statevector_in, statevector_out)
727 return
728 end if
729
730 if ( .not. hco_equal(statevector_in%hco, statevector_out%hco) ) then
731 call utl_abort('vInterp_gsv_r4: The input and output statevectors are not on the same horizontal grid.')
732 end if
733
734 ! DBGmad remove?
735 if ( gsv_getDataKind(statevector_in) /= 4 .or. gsv_getDataKind(statevector_out) /= 4 ) then
736 call utl_abort('vInterp_gsv_r4: Incorrect value for dataKind. Only compatible with dataKind=4')
737 end if
738
739 if (gsv_isAssocHeightSfc(statevector_in) .and. gsv_isAssocHeightSfc(statevector_out) ) then
740 heightSfcIn => gsv_getHeightSfc(statevector_in)
741 heightSfcOut => gsv_getHeightSfc(statevector_out)
742 heightSfcOut(:,:) = heightSfcIn(:,:)
743 end if
744
745 ! DBGmad move to int_vInterp_gsv?
746 ! the default is to ensure that the top of the output grid is ~equal or lower than the top of the input grid
747 if ( present(checkModelTop_opt) ) then
748 checkModelTop = checkModelTop_opt
749 else
750 checkModelTop = .true.
751 end if
752 if (checkModelTop) then
753 call msg('vInterp_gsv_r4', ' Checking that that the top of the destination grid is not higher than the top of the source grid.')
754 if ( vcode_in == 21001 .or. vcode_out == 21001 ) then
755 call msg('vInterp_gsv_r4', 'bypassing top check, '&
756 //'vcode_in='//str(vcode_in)//', vcode_out='//str(vcode_out))
757 ! Development notes (@mad001)
758 ! we should consider having a new criterion that works for GEM-H as well
759 else
760 call czp_ensureCompatibleTops(vco_in, vco_out)
761 end if
762 end if
763
764 ! create reference state on output vco with reference surface fields
765 call gsv_allocate(statevectorRef_out, statevectorRef%numstep, &
766 statevectorRef%hco, statevector_out%vco, &
767 mpi_local_opt=statevectorRef%mpi_local, &
768 mpi_distribution_opt=statevectorRef%mpi_distribution, &
769 dataKind_opt=gsv_getDataKind(statevectorRef), &
770 allocHeightSfc_opt=statevectorRef%heightSfcPresent, &
771 varNames_opt=(/'P0','P0LS'/) )
772 call gsv_copy(stateVectorRef, stateVectorRef_out, allowVcoMismatch_opt=.true., &
773 allowVarMismatch_opt=.true.)
774
775 var_loop: do varIndex = 1, vnl_numvarmax
776 varName = vnl_varNameList(varIndex)
777 if ( .not. gsv_varExist(statevector_in,varName) ) cycle var_loop
778
779 nlev_in = statevector_in%varNumLev(varIndex)
780 nlev_out = statevector_out%varNumLev(varIndex)
781
782 call gsv_getField(statevector_in ,field_in,varName)
783 call gsv_getField(statevector_out,field_out,varName)
784
785 ! for 2D fields and variables on "other" levels, just copy and cycle to next variable
786 if ( (nlev_in == 1 .and. nlev_out == 1) .or. &
787 (vnl_varLevelFromVarname(varName) == 'OT') ) then
788 field_out(:,:,:,:) = field_in(:,:,:,:)
789 cycle var_loop
790 end if
791
792 hLikeCalcIf: if (hLikeCalc) then
793 hLikeCalc = .false.
794
795 ! allocating for temporary pressure references
796 allocate(hLikeT_in( statevector_in%myLonBeg:statevector_in%myLonEnd, &
797 statevector_in%myLatBeg:statevector_in%myLatEnd, &
798 gsv_getNumLev(statevector_in,'TH'), statevector_in%numStep))
799 allocate(hLikeM_in( statevector_in%myLonBeg:statevector_in%myLonEnd, &
800 statevector_in%myLatBeg:statevector_in%myLatEnd, &
801 gsv_getNumLev(statevector_in,'MM'), statevector_in%numStep))
802 allocate(hLikeT_out(statevector_out%myLonBeg:statevector_out%myLonEnd, &
803 statevector_out%myLatBeg:statevector_out%myLatEnd, &
804 gsv_getNumLev(statevector_out,'TH'), statevector_out%numStep))
805 allocate(hLikeM_out(statevector_out%myLonBeg:statevector_out%myLonEnd, &
806 statevector_out%myLatBeg:statevector_out%myLatEnd, &
807 gsv_getNumLev(statevector_out,'MM'), statevector_out%numStep))
808 allocate(tmpCoord_T(statevector_in%myLonBeg:statevector_in%myLonEnd, &
809 statevector_in%myLatBeg:statevector_in%myLatEnd, &
810 gsv_getNumLev(statevector_in,'TH'), statevector_in%numStep))
811 allocate(tmpCoord_M(statevector_in%myLonBeg:statevector_in%myLonEnd, &
812 statevector_in%myLatBeg:statevector_in%myLatEnd, &
813 gsv_getNumLev(statevector_in,'MM'), statevector_in%numStep))
814
815 ! output grid GEM-P interpolation in log-pressure
816 if ( vcode_out==5002 .or. vcode_out==5005 .or. vcode_out==5100 ) then
817 call czp_calcReturnPressure_gsv_nl( statevectorRef_out, &
818 PTout_r4_opt=hLikeT_out, &
819 PMout_r4_opt=hLikeM_out, &
820 Ps_in_hPa_opt=Ps_in_hPa_opt)
821
822 if ( vcode_in==5002 .or. vcode_in==5005 .or. vcode_in==5100 ) then
823 call czp_calcReturnPressure_gsv_nl( statevectorRef, &
824 PTout_r4_opt=hLikeT_in, &
825 PMout_r4_opt=hLikeM_in, &
826 Ps_in_hPa_opt=Ps_in_hPa_opt)
827 else if ( vcode_in==21001 ) then
828 call czp_calcReturnHeight_gsv_nl( statevectorRef, &
829 ZTout_r4_opt=tmpCoord_T, &
830 ZMout_r4_opt=tmpCoord_M)
831 call czp_calcReturnPressure_gsv_nl( statevectorRef, &
832 ZTin_r4_opt=tmpCoord_T, &
833 ZMin_r4_opt=tmpCoord_M, &
834 PTout_r4_opt=hLikeT_in, &
835 PMout_r4_opt=hLikeM_in, &
836 Ps_in_hPa_opt=Ps_in_hPa_opt)
837 end if
838
839 call msg('vInterp_gsv_r4','converting pressure coordinates to height-like, '&
840 //'vcode_in='//str(vcode_in)//', vcode_out='//str(vcode_out))
841 call logP_r4(hLikeM_out)
842 call logP_r4(hLikeT_out)
843 call logP_r4(hLikeM_in)
844 call logP_r4(hLikeT_in)
845
846 ! output grid GEM-H interpolation in height
847 else if ( vcode_out==21001 ) then
848 call czp_calcReturnHeight_gsv_nl( statevectorRef_out, &
849 ZTout_r4_opt=hLikeT_out, &
850 ZMout_r4_opt=hLikeM_out)
851 if ( vcode_in==21001 ) then
852 call czp_calcReturnHeight_gsv_nl( statevectorRef, &
853 ZTout_r4_opt=hLikeT_in, &
854 ZMout_r4_opt=hLikeM_in)
855 else if ( vcode_in==5002 .or. vcode_in==5005 .or. vcode_in==5100 ) then
856 call czp_calcReturnPressure_gsv_nl( statevectorRef, &
857 PTout_r4_opt=tmpCoord_T, &
858 PMout_r4_opt=tmpCoord_M, &
859 Ps_in_hPa_opt=Ps_in_hPa_opt)
860 call czp_calcReturnHeight_gsv_nl( statevectorRef, &
861 PTin_r4_opt=tmpCoord_T, &
862 PMin_r4_opt=tmpCoord_M, &
863 ZTout_r4_opt=hLikeT_in, &
864 ZMout_r4_opt=hLikeM_in)
865 end if
866 end if
867 deallocate(tmpCoord_T)
868 deallocate(tmpCoord_M)
869
870 end if hLikeCalcIf
871
872 step_loop: do stepIndex = 1, statevector_out%numStep
873
874 ! copy over some time related and other parameters
875 statevector_out%deet = statevector_in%deet
876 statevector_out%dateOriginList(stepIndex) = statevector_in%dateOriginList(stepIndex)
877 statevector_out%npasList(stepIndex) = statevector_in%npasList(stepIndex)
878 statevector_out%ip2List(stepIndex) = statevector_in%ip2List(stepIndex)
879 statevector_out%etiket = statevector_in%etiket
880 statevector_out%onPhysicsGrid(:) = statevector_in%onPhysicsGrid(:)
881 statevector_out%hco_physics => statevector_in%hco_physics
882
883 ! do the vertical interpolation
884 field_out(:,:,:,stepIndex) = 0.0d0
885 if (vnl_varLevelFromVarname(varName) == 'TH') then
886 call hLike_interpolation_r4(hLikeT_in, hLikeT_out)
887 else
888 call hLike_interpolation_r4(hLikeM_in, hLikeM_out)
889 end if
890
891 ! overwrite values at the lowest levels to avoid extrapolation
892 if (vInterpCopyLowestLevel) then
893 field_out(:,:,nlev_out,stepIndex) = field_in(:,:,nlev_in,stepIndex)
894 end if
895
896 end do step_loop
897
898 end do var_loop
899
900 if (associated(hLikeT_in)) then
901 deallocate(hLikeT_in, hLikeM_in, hLikeT_out, hLikeM_out)
902 end if
903 call gsv_deallocate(statevectorRef_out)
904
905 call msg('vInterp_gsv_r4', 'END', verb_opt=4)
906
907 contains
908
909 subroutine hLike_interpolation_r4(hLike_in, hLike_out)
910 !
911 !:Purpose: Proceed to actual interpolation in H-logP representation
912 !
913 implicit none
914
915 ! Arguments:
916 real(4), pointer, intent(in) :: hLike_in(:,:,:,:) ! abstract height dimensioned input coordinate
917 real(4), pointer, intent(in) :: hLike_out(:,:,:,:) ! abstract height dimensioned target coordinate
918
919 ! Locals:
920 integer :: latIndex, lonIndex, levIndex_out, levIndex_in
921 real(4) :: zwb, zwt
922
923 !$OMP PARALLEL DO PRIVATE(latIndex,lonIndex,levIndex_in,levIndex_out,zwb,zwt)
924 do latIndex = statevector_out%myLatBeg, statevector_out%myLatEnd
925 do lonIndex = statevector_out%myLonBeg, statevector_out%myLonEnd
926 levIndex_in = 1
927 do levIndex_out = 1, nlev_out
928 levIndex_in = levIndex_in + 1
929 do while(hLike_out(lonIndex,latIndex,levIndex_out,stepIndex) &
930 .gt.hLike_in(lonIndex,latIndex,levIndex_in,stepIndex) &
931 .and.levIndex_in.lt.nlev_in)
932 levIndex_in = levIndex_in + 1
933 end do
934 levIndex_in = levIndex_in - 1
935 zwb = (hLike_out(lonIndex,latIndex,levIndex_out,stepIndex) &
936 - hLike_in(lonIndex,latIndex,levIndex_in,stepIndex)) &
937 /(hLike_in(lonIndex,latIndex,levIndex_in+1,stepIndex) &
938 - hLike_in(lonIndex,latIndex,levIndex_in,stepIndex))
939 zwt = 1.d0 - zwb
940 field_out(lonIndex,latIndex,levIndex_out,stepIndex) = &
941 zwb*field_in(lonIndex,latIndex,levIndex_in+1,stepIndex) &
942 + zwt*field_in(lonIndex,latIndex,levIndex_in,stepIndex)
943 end do
944 end do
945 end do
946 !$OMP END PARALLEL DO
947
948 end subroutine hLike_interpolation_r4
949
950 end subroutine vInterp_gsv_r4
951
952 !--------------------------------------------------------------------------
953 ! logP_r4
954 !--------------------------------------------------------------------------
955 subroutine logP_r4(presInLogOut)
956 !
957 !:Purpose: compute log of pressurce field
958 !
959 implicit none
960
961 ! Arguments:
962 real(4), intent(inout) :: presInLogOut(:,:,:,:)
963
964 ! Locals:
965 integer :: latIdx, lonIdx, levIdx, stepIdx
966
967 !$OMP PARALLEL DO PRIVATE(lonIdx,latIdx,levIdx,stepIdx)
968 do stepIdx = lbound(presInLogOut,4), ubound(presInLogOut,4)
969 do levIdx = lbound(presInLogOut,3), ubound(presInLogOut,3)
970 do latIdx = lbound(presInLogOut,2), ubound(presInLogOut,2)
971 do lonIdx = lbound(presInLogOut,1), ubound(presInLogOut,1)
972 presInLogOut(lonIdx,latIdx,levIdx,stepIdx) = &
973 log(presInLogOut(lonIdx,latIdx,levIdx,stepIdx))
974 end do
975 end do
976 end do
977 end do
978 !$OMP END PARALLEL DO
979
980 end subroutine logP_r4
981
982 !--------------------------------------------------------------------------
983 ! int_tInterp_gsv
984 !--------------------------------------------------------------------------
985 subroutine int_tInterp_gsv(statevector_in,statevector_out)
986 !
987 !:Purpose: Time interpolation from statevector with low temporal resolution
988 ! to statevector with high temporal resolution.
989 !
990 implicit none
991
992 ! Arguments:
993 type(struct_gsv), intent(in) :: statevector_in ! Statevector input
994 type(struct_gsv), intent(inout) :: statevector_out ! Statevector with target temporal structure and results
995
996 ! Locals:
997 integer :: kIndex, latIndex, lonIndex
998 integer :: stepIndexIn1, stepIndexIn2, stepIndexOut, numStepIn, numStepOut
999 integer :: lon1, lon2, lat1, lat2, k1, k2
1000 integer :: dateStampIn, dateStampOut
1001 real(8) :: weight1, weight2
1002 real(8) :: deltaHour, deltaHourInOut
1003 real(4), pointer :: gdIn_r4(:,:,:,:), gdOut_r4(:,:,:,:)
1004 real(8), pointer :: gdIn_r8(:,:,:,:), gdOut_r8(:,:,:,:)
1005
1006 call msg('int_tInterp_gsv', 'START', verb_opt=2)
1007
1008 ! read the namelist
1009 call int_readNml()
1010
1011 if ( .not. vco_equal(statevector_in%vco, statevector_out%vco) ) then
1012 call utl_abort('int_tInterp_gsv: The input and output statevectors are not on the same vertical levels')
1013 end if
1014
1015 if ( .not. hco_equal(statevector_in%hco, statevector_out%hco) ) then
1016 call utl_abort('int_tInterp_gsv: The input and output statevectors are not on the same horizontal grid.')
1017 end if
1018
1019 if ( statevector_in%numStep > statevector_out%numStep ) then
1020 call msg('int_tInterp_gsv', 'numStep_out is less than numStep_in, calling gsv_copy.')
1021 call gsv_copy(statevector_in, statevector_out, allowTimeMismatch_opt=.true.)
1022 return
1023 else if ( statevector_in%numStep == statevector_out%numStep ) then
1024 call msg('int_tInterp_gsv', 'numStep_out is equal to numStep_in, calling gsv_copy.')
1025 call gsv_copy(statevector_in, statevector_out, allowVarMismatch_opt=.true.)
1026 return
1027 end if
1028
1029 lon1 = statevector_in%myLonBeg
1030 lon2 = statevector_in%myLonEnd
1031 lat1 = statevector_in%myLatBeg
1032 lat2 = statevector_in%myLatEnd
1033 k1 = statevector_in%mykBeg
1034 k2 = statevector_in%mykEnd
1035
1036 numStepIn = statevector_in%numStep
1037 numStepOut = statevector_out%numStep
1038 call msg('int_tInterp_gsv', 'numStepIn='//str(numStepIn)&
1039 //', numStepOut='//str(numStepOut), mpiAll_opt=.false.)
1040
1041 ! compute positive deltaHour between two first stepIndex of statevector_in (input temporal grid).
1042 ! If numStepIn == 1, no time interpolation needed (weights are set to zero).
1043 if ( numStepIn > 1 ) then
1044 call difdatr(statevector_in%dateStampList(1), statevector_in%dateStampList(2), deltaHour)
1045 deltaHour = abs(deltaHour)
1046 else
1047 deltaHour = 0.0d0
1048 end if
1049
1050 do stepIndexOut = 1, numStepOut
1051 if ( numStepIn == 1 ) then
1052 stepIndexIn1 = 1
1053 stepIndexIn2 = 1
1054 weight1 = 1.0d0
1055 weight2 = 0.0d0
1056 deltaHourInOut = 0.0d0
1057 else
1058 ! find staevector_in%dateStamp on the left and right
1059 if ( statevector_in%dateStampList(numStepIn) == statevector_out%dateStampList(stepIndexOut) ) then
1060 stepIndexIn2 = numStepIn
1061 else
1062 stepInLoop: do stepIndexIn2 = 1, numStepIn
1063 dateStampIn = statevector_in%dateStampList(stepIndexIn2)
1064 dateStampOut = statevector_out%dateStampList(stepIndexOut)
1065 call difdatr(dateStampIn, dateStampOut, deltaHourInOut)
1066 if ( deltaHourInOut > 0.0d0 ) exit stepInLoop
1067 end do stepInLoop
1068 end if
1069 stepIndexIn1 = stepIndexIn2 - 1
1070
1071 ! compute deltaHour between left stepIndex of statevector_in and statevector_out
1072 dateStampIn = statevector_in%dateStampList(stepIndexIn1)
1073 dateStampOut = statevector_out%dateStampList(stepIndexOut)
1074 call difdatr(dateStampOut, dateStampIn, deltaHourInOut)
1075 if ( deltaHourInOut < 0.0d0 ) then
1076 call utl_abort('int_tInterp_gsv: deltaHourInOut should be greater or equal to 0')
1077 end if
1078
1079 ! compute the interpolation weights for left stepIndex of statevector_in (weight1) and right stepIndex of statevector_in (weight2)
1080 weight1 = 1.0d0 - deltaHourInOut / deltaHour
1081 weight2 = deltaHourInOut / deltaHour
1082 end if
1083
1084 call msg('int_tInterp_gsv', 'for stepIndexOut='//str(stepIndexOut) &
1085 //', stepIndexIn1='//str(stepIndexIn1)//', stepIndexIn2='//str(stepIndexIn2) &
1086 //', weight1='//str(weight1)//', weight2='//str(weight2) &
1087 //', deltaHourInOut/deltaHour='//str(deltaHourInOut)//'/'//str(deltaHour), &
1088 mpiAll_opt=.false.)
1089
1090 if ( gsv_getDataKind(statevector_in) == 4 .and. gsv_getDataKind(statevector_out) == 4 ) then
1091 call gsv_getField(statevector_in, gdIn_r4)
1092 call gsv_getField(statevector_out, gdOut_r4)
1093 !$OMP PARALLEL DO PRIVATE (latIndex,kIndex,lonIndex)
1094 do kIndex = k1, k2
1095 do latIndex = lat1, lat2
1096 do lonIndex = lon1, lon2
1097 gdOut_r4(lonIndex,latIndex,kIndex,stepIndexOut) = &
1098 real(weight1,4) * gdIn_r4(lonIndex,latIndex,kIndex,stepIndexIn1) + &
1099 real(weight2,4) * gdIn_r4(lonIndex,latIndex,kIndex,stepIndexIn2)
1100 end do
1101 end do
1102 end do
1103 !$OMP END PARALLEL DO
1104 else if ( gsv_getDataKind(statevector_in) == 4 .and. gsv_getDataKind(statevector_out) == 8 ) then
1105 call gsv_getField(statevector_in, gdIn_r4)
1106 call gsv_getField(statevector_out, gdOut_r8)
1107 !$OMP PARALLEL DO PRIVATE (latIndex,kIndex,lonIndex)
1108 do kIndex = k1, k2
1109 do latIndex = lat1, lat2
1110 do lonIndex = lon1, lon2
1111 gdOut_r8(lonIndex,latIndex,kIndex,stepIndexOut) = &
1112 weight1 * real(gdIn_r4(lonIndex,latIndex,kIndex,stepIndexIn1),8) + &
1113 weight2 * real(gdIn_r4(lonIndex,latIndex,kIndex,stepIndexIn2),8)
1114 end do
1115 end do
1116 end do
1117 !$OMP END PARALLEL DO
1118 else if ( gsv_getDataKind(statevector_in) == 8 .and. gsv_getDataKind(statevector_out) == 4 ) then
1119 call gsv_getField(statevector_in, gdIn_r8)
1120 call gsv_getField(statevector_out, gdOut_r4)
1121 !$OMP PARALLEL DO PRIVATE (latIndex,kIndex,lonIndex)
1122 do kIndex = k1, k2
1123 do latIndex = lat1, lat2
1124 do lonIndex = lon1, lon2
1125 gdOut_r4(lonIndex,latIndex,kIndex,stepIndexOut) = &
1126 real(weight1 * gdIn_r8(lonIndex,latIndex,kIndex,stepIndexIn1),4) + &
1127 real(weight2 * gdIn_r8(lonIndex,latIndex,kIndex,stepIndexIn2),4)
1128 end do
1129 end do
1130 end do
1131 !$OMP END PARALLEL DO
1132 else if ( gsv_getDataKind(statevector_in) == 8 .and. gsv_getDataKind(statevector_out) == 8 ) then
1133 call gsv_getField(statevector_in, gdIn_r8)
1134 call gsv_getField(statevector_out, gdOut_r8)
1135 !$OMP PARALLEL DO PRIVATE (latIndex,kIndex,lonIndex)
1136 do kIndex = k1, k2
1137 do latIndex = lat1, lat2
1138 do lonIndex = lon1, lon2
1139 gdOut_r8(lonIndex,latIndex,kIndex,stepIndexOut) = &
1140 weight1 * gdIn_r8(lonIndex,latIndex,kIndex,stepIndexIn1) + &
1141 weight2 * gdIn_r8(lonIndex,latIndex,kIndex,stepIndexIn2)
1142 end do
1143 end do
1144 end do
1145 !$OMP END PARALLEL DO
1146 end if
1147
1148 end do
1149
1150 call msg('int_tInterp_gsv', 'END', verb_opt=2)
1151
1152 end subroutine int_tInterp_gsv
1153
1154 !--------------------------------------------------------------------------
1155 ! int_vInterp_col
1156 !--------------------------------------------------------------------------
1157 subroutine int_vInterp_col(column_in,column_out,varName,useColumnPressure_opt)
1158 !
1159 !:Purpose: Vertical interpolation of a columData object
1160 !
1161 implicit none
1162
1163 ! Arguments:
1164 type(struct_columnData), intent(in) :: column_in ! ColumnData input
1165 type(struct_columnData), intent(inout) :: column_out ! columnData with the vert structure and results
1166 character(len=*), intent(in) :: varName ! variable name to be interpolated
1167 logical, optional, intent(in) :: useColumnPressure_opt ! if .true. use P_* instead of the pressure provided by calcHeightAndPressure_mod
1168
1169 ! Locals:
1170 real(8), pointer :: varInterp_in(:), varInterp_out(:)
1171 real(8), pointer :: coordRef_in(:) , coordRef_out(:)
1172 real(8), pointer :: ZT_in(:,:), ZM_in(:,:), PT_in(:,:), PM_in(:,:)
1173 real(8), pointer :: PT_out(:,:), PM_out(:,:)
1174 character(len=4) :: varLevel
1175 real(8) :: zwb, zwt
1176 integer :: levIndex_out, levIndex_in, columnIndex
1177 integer :: vcode_in, vcode_out
1178 integer :: nLevIn_T, nLevIn_M, nLevOut_T, nLevOut_M
1179 logical :: vInterp, useColumnPressure
1180 integer, allocatable, target :: THlevelWanted(:), MMlevelWanted(:)
1181 integer, pointer :: levelWanted(:)
1182
1183 call msg('int_vInterp_col', 'START', verb_opt=2)
1184 varLevel = vnl_varLevelFromVarname(varName)
1185
1186 if ( present(useColumnPressure_opt) ) then
1187 useColumnPressure = useColumnPressure_opt
1188 else
1189 useColumnPressure = .true.
1190 end if
1191
1192 call msg('int_vInterp_col', varName//' ('//varLevel &
1193 //'), useColumnPressure='//str(useColumnPressure), verb_opt=3)
1194
1195 vInterp = .true.
1196 if ( .not. col_varExist(column_in,'P0' ) ) then
1197 call msg('int_vInterp_col', 'P0 is missing. Vertical interpolation WILL NOT BE PERFORMED')
1198 vInterp = .false.
1199 else if ( col_getNumLev(column_in ,'TH') <= 1 .or. &
1200 col_getNumLev(column_in ,'MM') <= 1 ) then
1201 vInterp = .false.
1202 call msg('int_vInterp_col', 'The input backgrounds are 2D. Vertical interpolation WILL NOT BE PERFORMED')
1203 end if
1204
1205 vcode_in = column_in%vco%vcode
1206 vcode_out = column_out%vco%vcode
1207 nLevIn_T = col_getNumLev(column_in, 'TH')
1208 nLevIn_M = col_getNumLev(column_in, 'MM')
1209 nLevOut_T = col_getNumLev(column_out, 'TH')
1210 nLevOut_M = col_getNumLev(column_out, 'MM')
1211
1212 if_vInterp: if (vInterp) then
1213 if ( .not. useColumnPressure ) then
1214
1215 call msg('int_vInterp_col', 'vcode_in='//str(vcode_in)//', '&
1216 //'vcode_out='//str(vcode_out), verb_opt=3)
1217
1218 ! Compute interpolation variables (pressures and or heights)
1219 allocate(PT_in( col_getNumCol(column_in), nLevIn_T))
1220 allocate(PM_in( col_getNumCol(column_in), nLevIn_M))
1221 allocate(ZT_in( col_getNumCol(column_in), nLevIn_T))
1222 allocate(ZM_in( col_getNumCol(column_in), nLevIn_M))
1223 allocate(PT_out(col_getNumCol(column_out), nLevOut_T))
1224 allocate(PM_out(col_getNumCol(column_out), nLevOut_M))
1225
1226 call czp_calcReturnPressure_col_nl(column_in, PT_in, PM_in)
1227 call czp_calcReturnHeight_col_nl(column_in, ZT_in, ZM_in)
1228 call czp_calcReturnPressure_col_nl(column_out, PT_out, PM_out)
1229
1230 end if ! useColumnPressure
1231
1232 do columnIndex = 1, col_getNumCol(column_out)
1233
1234 ! coordRef_{in,out}
1235 if ( varLevel == 'TH' ) then
1236 if ( .not. useColumnPressure ) then
1237 coordRef_in => PT_in(columnIndex,:)
1238 coordRef_out => PT_out(columnIndex,:)
1239 else
1240 coordRef_in => col_getColumn(column_in,columnIndex,'P_T')
1241 coordRef_out => col_getColumn(column_out,columnIndex,'P_T')
1242 end if
1243 else if ( varLevel == 'MM' ) then
1244 if ( .not. useColumnPressure ) then
1245 coordRef_in => PM_in(columnIndex,:)
1246 coordRef_out => PM_out(columnIndex,:)
1247 else
1248 coordRef_in => col_getColumn(column_in,columnIndex,'P_M')
1249 coordRef_out => col_getColumn(column_out,columnIndex,'P_M')
1250 end if
1251 else
1252 call utl_abort('int_vInterp_col: only varLevel TH/MM is allowed')
1253 end if
1254
1255 varInterp_in => col_getColumn(column_in ,columnIndex,varName)
1256 varInterp_out => col_getColumn(column_out,columnIndex,varName)
1257
1258 if ( columnIndex == 1 .and. (trim(varName) == 'P_T' ) ) then
1259 call msg('int_vInterp_col', 'useColumnPressure='//str(useColumnPressure) &
1260 //', '//trim(varName)//':' &
1261 //new_line('')//'COLUMN_IN(1):'//str(varInterp_in(:))&
1262 //new_line('')//'COLUMN_OUT(1):'//str(varInterp_out(:)), &
1263 mpiAll_opt=.false.)
1264 end if
1265
1266 ! actual interpolation
1267 ! Development notes (@mad001)
1268 ! Potential issue with GEM-H height based interpolation
1269 ! we should consider to convert to pure logP interpolation
1270 ! as we did for int_vInterp_gsv
1271 ! see also #466 https://gitlab.science.gc.ca/atmospheric-data-assimilation/midas/issues/466#note_497052
1272 levIndex_in = 1
1273 do levIndex_out = 1, col_getNumLev(column_out,varLevel)
1274 levIndex_in = levIndex_in + 1
1275 do while( coordRef_out(levIndex_out) .gt. coordRef_in(levIndex_in) .and. &
1276 levIndex_in .lt. col_getNumLev(column_in,varLevel) )
1277 levIndex_in = levIndex_in + 1
1278 end do
1279 levIndex_in = levIndex_in - 1
1280 zwb = log(coordRef_out(levIndex_out)/coordRef_in(levIndex_in))/ &
1281 log(coordRef_in(levIndex_in+1)/coordRef_in(levIndex_in))
1282 zwt = 1. - zwb
1283 if ( useColumnPressure .and. &
1284 (trim(varName) == 'P_T' .or. trim(varName) == 'P_M' ) ) then
1285 ! do nothing, i.e. use the pressures from column_in
1286 else if ( .not. useColumnPressure .and. &
1287 (trim(varName) == 'P_T' .or. trim(varName) == 'P_M' ) ) then
1288 varInterp_out(levIndex_out) = exp(zwb*log(varInterp_in(levIndex_in+1)) + zwt*log(varInterp_in(levIndex_in)))
1289 else
1290 varInterp_out(levIndex_out) = zwb*varInterp_in(levIndex_in+1) + zwt*varInterp_in(levIndex_in)
1291 end if
1292 end do
1293 end do
1294
1295 call msg('int_vInterp_col', trim(varName)//' (in): '//str(varInterp_in)&
1296 //new_line('')//trim(varName)//' (out): '//str(varInterp_out), verb_opt=3)
1297
1298 else if_vInterp
1299
1300 if (column_out%vco%nlev_T > 0 .and. column_out%vco%nlev_M > 0) then
1301
1302 ! Find which levels in column_in matches column_out
1303 allocate(THlevelWanted(column_out%vco%nlev_T))
1304 allocate(MMlevelWanted(column_out%vco%nlev_M))
1305
1306 call vco_levelMatchingList( THlevelWanted, MMlevelWanted, & ! OUT
1307 column_out%vco, column_in%vco ) ! IN
1308
1309 if ( any(THlevelWanted == -1) .or. any(MMlevelWanted == -1) ) then
1310 call utl_abort('int_vInterp_col: column_out is not a subsets of column_in!')
1311 end if
1312
1313 ! Transfer the corresponding data
1314 do columnIndex = 1, col_getNumCol(column_out)
1315 varInterp_in => col_getColumn(column_in ,columnIndex,varName)
1316 varInterp_out => col_getColumn(column_out,columnIndex,varName)
1317 if (vnl_varLevelFromVarname(varName) == 'TH') then
1318 levelWanted => THlevelWanted
1319 else
1320 levelWanted => MMlevelWanted
1321 end if
1322 do levIndex_out = 1, col_getNumLev(column_out,varLevel)
1323 varInterp_out(levIndex_out) = varInterp_in(levelWanted(levIndex_out))
1324 end do
1325 end do
1326
1327 deallocate(THlevelWanted)
1328 deallocate(MMlevelWanted)
1329
1330 else if (column_out%vco%nlev_depth > 0) then
1331 call msg('int_vInterp_col', 'vco_levelMatchingList: no MM and TH levels, but depth levels exist')
1332 if (any(column_out%vco%depths(:) /= column_in%vco%depths(:))) then
1333 call utl_abort('int_vInterp_col: some depth levels not equal')
1334 else
1335 ! copy over depth levels
1336 do columnIndex = 1, col_getNumCol(column_out)
1337 varInterp_in => col_getColumn(column_in ,columnIndex,varName)
1338 varInterp_out => col_getColumn(column_out,columnIndex,varName)
1339 do levIndex_out = 1, col_getNumLev(column_out,varLevel)
1340 varInterp_out(levIndex_out) = varInterp_in(levIndex_out)
1341 end do
1342 end do
1343 end if
1344
1345 end if
1346
1347 end if if_vInterp
1348
1349 call msg('int_vInterp_col', 'END', verb_opt=2)
1350 end subroutine int_vInterp_col
1351
1352 !--------------------------------------------------------------------------
1353 ! int_setezopt
1354 !--------------------------------------------------------------------------
1355 subroutine int_setezopt(interpDegree, extrapDegree_opt)
1356 !
1357 !:Purpose: Wrapper subroutine for rmnlib routine setezopt.
1358 !
1359 implicit none
1360
1361 ! Arguments:
1362 character(len=*), intent(in) :: interpDegree
1363 character(len=*), optional, intent(in) :: extrapDegree_opt
1364
1365 ! Locals:
1366 character(len=12) :: extrapDegree
1367 integer :: ierr, ezsetopt, ezsetval
1368
1369 call msg('int_setezopt', 'START', verb_opt=4)
1370 if ( trim(interpDegree) /= 'LINEAR' .and. &
1371 trim(interpDegree) /= 'CUBIC' .and. &
1372 trim(interpDegree) /= 'NEAREST' ) then
1373
1374 call utl_abort('int_setezopt: invalid interpolation degree = '//trim(interpDegree))
1375 end if
1376
1377 if ( present(extrapDegree_opt) ) then
1378 extrapDegree = extrapDegree_opt
1379 else
1380 extrapDegree = 'VALUE'
1381 end if
1382
1383 ierr = ezsetopt('INTERP_DEGREE', interpDegree)
1384 if ( trim(extrapDegree) == 'VALUE' ) then
1385 ierr = ezsetval('EXTRAP_VALUE', 0.0)
1386 end if
1387 ierr = ezsetopt('EXTRAP_DEGREE', extrapDegree)
1388
1389 call msg('int_setezopt', 'END', verb_opt=4)
1390 end subroutine int_setezopt
1391
1392 !--------------------------------------------------------------------------
1393 ! int_hInterpScalar_gsv
1394 !--------------------------------------------------------------------------
1395 function int_hInterpScalar_gsv(stateVectorOut, stateVectorIn, varName, levIndex, stepIndex, &
1396 interpDegree, extrapDegree_opt) result(ierr)
1397 !
1398 !:Purpose: Horizontal interpolation of 2D scalar field that use stateVector
1399 ! objects for input and output. Accessed through int_hInterpScalar.
1400 !
1401 implicit none
1402
1403 ! Arguments:
1404 type(struct_gsv), intent(inout) :: stateVectorOut
1405 type(struct_gsv), intent(inout) :: stateVectorIn
1406 character(len=*), intent(in) :: varName
1407 integer , intent(in) :: levIndex
1408 integer , intent(in) :: stepIndex
1409 character(len=*), intent(in) :: interpDegree
1410 character(len=*), optional, intent(in) :: extrapDegree_opt
1411 ! Result:
1412 integer :: ierr
1413
1414 ! Locals:
1415 real(4), pointer :: fieldOut_r4(:,:,:,:), fieldIn_r4(:,:,:,:)
1416 real(8), pointer :: fieldOut_r8(:,:,:,:), fieldIn_r8(:,:,:,:)
1417 real(8), pointer :: heightSfcOut(:,:), heightSfcIn(:,:)
1418 integer :: ezsint, ezdefset
1419
1420 call msg('int_hInterpScalar_gsv', 'START', verb_opt=4)
1421 ! read the namelist
1422 call int_readNml()
1423
1424 ! check if special interpolation is required
1425 if (stateVectorIn%hco%initialized .and. stateVectorOut%hco%initialized) then
1426 if (stateVectorIn%hco%grtyp == 'Y') then
1427 ! for now, only comptable for real(4)
1428 if( gsv_getDataKind(stateVectorOut) /= 4 .or. gsv_getDataKind(stateVectorIn) /= 4) then
1429 call utl_abort('int_hInterpScalar_gsv: cloudToGrid only implemented for real(4)')
1430 end if
1431 ierr = int_sintCloudToGrid_gsv(stateVectorOut, stateVectorIn, varName, levIndex, stepIndex)
1432 return
1433 end if
1434 end if
1435
1436 ! do the standard interpolation
1437
1438 ierr = ezdefset(stateVectorOut%hco%EZscintID, stateVectorIn%hco%EZscintID)
1439 call int_setezopt(interpDegree, extrapDegree_opt)
1440
1441 if (trim(varName) == 'ZSFC') then
1442
1443 heightSfcIn => gsv_getHeightSfc(stateVectorIn)
1444 heightSfcOut => gsv_getHeightSfc(stateVectorOut)
1445
1446 ! allocate real(4) buffers and copy to/from for interpolation
1447 allocate(fieldIn_r4(stateVectorIn%hco%ni,stateVectorIn%hco%nj,1,1))
1448 allocate(fieldOut_r4(stateVectorOut%hco%ni,stateVectorOut%hco%nj,1,1))
1449 fieldIn_r4(:,:,1,1) = heightSfcIn(:,:)
1450 ierr = ezsint(fieldOut_r4(:,:,1,1),fieldIn_r4(:,:,1,1))
1451 heightSfcOut(:,:) = fieldOut_r4(:,:,1,1)
1452 deallocate(fieldIn_r4,fieldOut_r4)
1453
1454 else if ( gsv_getDataKind(stateVectorOut) == 4 .and. gsv_getDataKind(stateVectorIn) == 4) then
1455
1456 if (trim(varName) == 'ALL') then
1457 call gsv_getField(stateVectorOut, fieldOut_r4)
1458 call gsv_getField(stateVectorIn, fieldIn_r4)
1459 else
1460 call gsv_getField(stateVectorOut, fieldOut_r4, varName)
1461 call gsv_getField(stateVectorIn, fieldIn_r4, varName)
1462 end if
1463
1464 ierr = ezsint(fieldOut_r4(:,:,levIndex,stepIndex),fieldIn_r4(:,:,levIndex,stepIndex))
1465
1466 else if ( gsv_getDataKind(stateVectorOut) == 8 .and. gsv_getDataKind(stateVectorIn) == 8) then
1467
1468 if (trim(varName) == 'ALL') then
1469 call gsv_getField(stateVectorOut, fieldOut_r8)
1470 call gsv_getField(stateVectorIn, fieldIn_r8)
1471 else
1472 call gsv_getField(stateVectorOut, fieldOut_r8, varName)
1473 call gsv_getField(stateVectorIn, fieldIn_r8, varName)
1474 end if
1475
1476 ! allocate real(4) buffers and copy to/from for interpolation
1477 allocate(fieldIn_r4(stateVectorIn%hco%ni,stateVectorIn%hco%nj,1,1))
1478 allocate(fieldOut_r4(stateVectorOut%hco%ni,stateVectorOut%hco%nj,1,1))
1479 fieldIn_r4(:,:,1,1) = fieldIn_r8(:,:,levIndex,stepIndex)
1480 ierr = ezsint(fieldOut_r4(:,:,1,1),fieldIn_r4(:,:,1,1))
1481 fieldOut_r8(:,:,levIndex,stepIndex) = fieldOut_r4(:,:,1,1)
1482 deallocate(fieldIn_r4,fieldOut_r4)
1483
1484 else
1485
1486 call utl_abort('int_hInterpScalar_gsv: not implemented for mixed dataKind')
1487
1488 end if
1489
1490 call msg('int_hInterpScalar_gsv', 'END', verb_opt=4)
1491 end function int_hInterpScalar_gsv
1492
1493 !--------------------------------------------------------------------------
1494 ! int_hInterpScalar_r4_2d
1495 !--------------------------------------------------------------------------
1496 function int_hInterpScalar_r4_2d(fieldOut_r4, fieldIn_r4, interpDegree, extrapDegree_opt) result(ierr)
1497 !
1498 !:Purpose: Horizontal interpolation of 2D scalar field that use real(4) arrays
1499 ! for input and output. Accessed through int_hInterpScalar.
1500 !
1501 implicit none
1502
1503 ! Arguments:
1504 real(4), intent(inout) :: fieldOut_r4(:,:)
1505 real(4), intent(in) :: fieldIn_r4(:,:)
1506 character(len=*), intent(in) :: interpDegree
1507 character(len=*), optional, intent(in) :: extrapDegree_opt
1508 ! Result:
1509 integer :: ierr
1510
1511 ! Locals:
1512 integer :: ezsint
1513
1514 call msg('int_hInterpScalar_r4_2d', 'START', verb_opt=4)
1515 ! read the namelist
1516 call int_readNml()
1517
1518 ! do the standard interpolation
1519 call int_setezopt(interpDegree, extrapDegree_opt)
1520 ierr = ezsint(fieldOut_r4,fieldIn_r4)
1521
1522 call msg('int_hInterpScalar_r4_2d', 'END', verb_opt=4)
1523 end function int_hInterpScalar_r4_2d
1524
1525 !--------------------------------------------------------------------------
1526 ! int_sintCloudToGrid_gsv
1527 !--------------------------------------------------------------------------
1528 function int_sintCloudToGrid_gsv(stateVectorGrid, stateVectorCloud, varName, levIndex, stepIndex) result(ierr)
1529 !
1530 !:Purpose: Perform horizontal interpolation for 1 level and time step (and variable)
1531 ! in the case where the input data is a cloud of points (i.e. a Y grid) and
1532 ! the output is on a regular grid. Accessed through int_hInterpScalar.
1533 !
1534 !:Note: When varName=='ALL', the argument levIndex is actually kIndex
1535 !
1536 implicit none
1537
1538 ! Arguments:
1539 type(struct_gsv), intent(inout) :: stateVectorGrid
1540 type(struct_gsv), intent(inout) :: stateVectorCloud
1541 character(len=*), intent(in) :: varName
1542 integer, intent(in) :: levIndex
1543 integer, intent(in) :: stepIndex
1544 ! Result:
1545 integer :: ierr
1546
1547 ! Locals:
1548 integer :: gdxyfll, omp_get_thread_num
1549 integer :: niCloud, njCloud, niGrid, njGrid, myThreadNum
1550 integer :: top, bottom, left, right, numBoxIndexes, lonIndexCloud, latIndexCloud
1551 integer :: boxSize, lonBoxIndex, latBoxIndex, boxIndex, lonIndexGrid, latIndexGrid
1552 integer :: lonBoxIndexes(100), latBoxIndexes(100), ngp
1553 integer :: nfill(mmpi_numThread), nhole(mmpi_numThread), nextrap0, nextrap1
1554 integer, allocatable :: numFilledByAvg(:,:), filledByInterp(:,:), maskGrid(:,:), maskCloud(:,:)
1555 real(4), pointer :: fieldCloud_4d(:,:,:,:), fieldGrid_4d(:,:,:,:)
1556 real(4), pointer :: fieldCloud(:,:), fieldGrid(:,:)
1557 real(4), allocatable :: fieldGrid_tmp(:,:)
1558 real(4), allocatable :: xCloud(:,:), yCloud(:,:)
1559 real(8) :: extrapValue
1560 integer, parameter :: maxNumLocalGridPointsSearch = 200000
1561 integer :: numLocalGridPointsFound, gridIndex
1562 type(kdtree2), pointer :: tree => null()
1563 real(kdkind), allocatable :: positionArray(:,:)
1564 type(kdtree2_result) :: searchResults(maxNumLocalGridPointsSearch)
1565 real(kdkind) :: searchRadiusSquared
1566 real(kdkind) :: refPosition(3)
1567
1568 call utl_tmg_start(176, 'low-level--int_sintCloudToGrid_gsv')
1569 call msg('int_sintCloudToGrid_gsv', 'START', verb_opt=2)
1570
1571 niCloud = stateVectorCloud%hco%ni
1572 njCloud = stateVectorCloud%hco%nj
1573 niGrid = stateVectorGrid%hco%ni
1574 njGrid = stateVectorGrid%hco%nj
1575
1576 if (trim(varName) == 'ALL') then
1577 call gsv_getField(stateVectorGrid, fieldGrid_4d)
1578 call gsv_getField(stateVectorCloud, fieldCloud_4d)
1579 else
1580 call gsv_getField(stateVectorGrid, fieldGrid_4d, varName)
1581 call gsv_getField(stateVectorCloud, fieldCloud_4d, varName)
1582 end if
1583 fieldGrid => fieldGrid_4d(:,:,levIndex,stepIndex)
1584 fieldCloud => fieldCloud_4d(:,:,levIndex,stepIndex)
1585
1586 allocate(xCloud(niCloud, njCloud))
1587 allocate(yCloud(niCloud, njCloud))
1588 allocate(numFilledByAvg(niGrid, njGrid))
1589 allocate(filledByInterp(niGrid, njGrid))
1590 allocate(maskCloud(niCloud, njCloud))
1591 allocate(maskGrid(niGrid, njGrid))
1592 allocate(fieldGrid_tmp(niGrid,njGrid))
1593
1594 ! set masks based on oceanMasks (if present)
1595 if (stateVectorCloud%oceanMask%maskPresent) then
1596 maskCloud(:,:) = 0
1597 if (trim(varName) == 'ALL') then
1598 ! when varName==ALL, the argument levIndex is actually kIndex
1599 where(stateVectorCloud%oceanMask%mask(:,:,gsv_getLevFromK(stateVectorCloud,levIndex))) maskCloud(:,:) = 1
1600 else
1601 where(stateVectorCloud%oceanMask%mask(:,:,levIndex)) maskCloud(:,:) = 1
1602 end if
1603 else
1604 maskCloud(:,:) = 1
1605 end if
1606 if (stateVectorGrid%oceanMask%maskPresent) then
1607 maskGrid(:,:) = 0
1608 if (trim(varName) == 'ALL') then
1609 ! when varName==ALL, the argument levIndex is actually kIndex
1610 where(stateVectorGrid%oceanMask%mask(:,:,gsv_getLevFromK(stateVectorGrid,levIndex))) maskGrid(:,:) = 1
1611 else
1612 where(stateVectorGrid%oceanMask%mask(:,:,levIndex)) maskGrid(:,:) = 1
1613 end if
1614 else
1615 maskGrid(:,:) = 1
1616 end if
1617
1618 ! Calcul des pos. x-y des eclairs sur la grille modele
1619 ierr = gdxyfll(stateVectorGrid%hco%EZscintID, xCloud, yCloud, &
1620 stateVectorCloud%hco%lat2d_4*MPC_DEGREES_PER_RADIAN_R4, &
1621 stateVectorCloud%hco%lon2d_4*MPC_DEGREES_PER_RADIAN_R4, &
1622 stateVectorCloud%hco%ni*stateVectorCloud%hco%nj)
1623
1624 ! Average values of cloud points neighbouring a grid location
1625 fieldGrid(:,:) = 0.0
1626 numFilledByAvg(:,:) = 0
1627 do latIndexCloud = 1, njCloud
1628 do lonIndexCloud = 1, niCloud
1629 lonIndexGrid = nint(xCloud(lonIndexCloud,latIndexCloud))
1630 latIndexGrid = nint(yCloud(lonIndexCloud,latIndexCloud))
1631 if ( lonIndexGrid >= 1 .and. latIndexGrid >= 1 .and. &
1632 lonIndexGrid <= niGrid .and. latIndexGrid <= njGrid ) then
1633 if ( maskCloud(lonIndexCloud,latIndexCloud) == 1 ) then
1634 fieldGrid(lonIndexGrid,latIndexGrid) = fieldGrid(lonIndexGrid,latIndexGrid) + &
1635 fieldCloud(lonIndexCloud,latIndexCloud)
1636 numFilledByAvg(lonIndexGrid,latIndexGrid) = numFilledByAvg(lonIndexGrid,latIndexGrid) + 1
1637 end if
1638 end if
1639 end do
1640 end do
1641
1642 do latIndexGrid = 1, njGrid
1643 do lonIndexGrid = 1, niGrid
1644 if(numFilledByAvg(lonIndexGrid,latIndexGrid) > 0) then
1645 fieldGrid(lonIndexGrid,latIndexGrid) = fieldGrid(lonIndexGrid,latIndexGrid)/ &
1646 real(numFilledByAvg(lonIndexGrid,latIndexGrid))
1647 end if
1648 end do
1649 end do
1650
1651 ! Now do something for grid points that don't have any value assigned
1652 fieldGrid_tmp(:,:) = fieldGrid(:,:)
1653 nfill(:) = 0
1654 nhole(:) = 0
1655 boxSizeLoop: do boxSize = 1, maxBoxSize
1656 filledByInterp(:,:) = 0
1657 !$OMP PARALLEL DO PRIVATE(latIndexGrid,lonIndexGrid,myThreadNum,top,bottom,left,right,numBoxIndexes,lonBoxIndex,latBoxIndex,lonBoxIndexes,latBoxIndexes,ngp,boxIndex)
1658 do latIndexGrid = 1, njGrid
1659 myThreadNum = 1 + omp_get_thread_num()
1660 do lonIndexGrid = 1, niGrid
1661 if (numFilledByAvg(lonIndexGrid,latIndexGrid) > 0) cycle
1662
1663 if (boxSize == 1) nhole(myThreadNum) = nhole(myThreadNum) + 1
1664
1665 top = latIndexGrid + boxSize
1666 bottom = latIndexGrid - boxSize
1667 left = lonIndexGrid - boxSize
1668 right = lonIndexGrid + boxSize
1669
1670 numBoxIndexes = 0
1671 latBoxIndex = bottom
1672 do lonBoxIndex = left, right
1673 numBoxIndexes = numBoxIndexes + 1
1674 lonBoxIndexes(numBoxIndexes) = lonBoxIndex
1675 latBoxIndexes(numBoxIndexes) = latBoxIndex
1676 end do
1677 lonBoxIndex = right
1678 do latBoxIndex = bottom + 1, top
1679 numBoxIndexes = numBoxIndexes + 1
1680 lonBoxIndexes(numBoxIndexes) = lonBoxIndex
1681 latBoxIndexes(numBoxIndexes) = latBoxIndex
1682 end do
1683 latBoxIndex = top
1684 do lonBoxIndex = right - 1, left, -1
1685 numBoxIndexes = numBoxIndexes + 1
1686 lonBoxIndexes(numBoxIndexes) = lonBoxIndex
1687 latBoxIndexes(numBoxIndexes) = latBoxIndex
1688 end do
1689 lonBoxIndex = left
1690 do latBoxIndex = top - 1, bottom + 1, -1
1691 numBoxIndexes = numBoxIndexes + 1
1692 lonBoxIndexes(numBoxIndexes) = lonBoxIndex
1693 latBoxIndexes(numBoxIndexes) = latBoxIndex
1694 end do
1695
1696 ngp = 0
1697 do boxIndex = 1, numBoxIndexes
1698
1699 if (lonBoxIndexes(boxIndex) >= 1 .and. lonBoxIndexes(boxIndex) <= niGrid .and. &
1700 latBoxIndexes(boxIndex) >= 1 .and. latBoxIndexes(boxIndex) <= njGrid) then
1701 if ( numFilledByAvg(lonBoxIndexes(boxIndex),latBoxIndexes(boxIndex)) > 0 ) then
1702
1703 fieldGrid_tmp(lonIndexGrid,latIndexGrid) = &
1704 fieldGrid_tmp(lonIndexGrid,latIndexGrid) + &
1705 fieldGrid(lonBoxIndexes(boxIndex),latBoxIndexes(boxIndex))
1706 ngp = ngp + 1
1707
1708 end if
1709 end if
1710
1711 end do
1712
1713 if (ngp /= 0) then
1714 ! mark the grid point as being filled by interpolation on the grid
1715 filledByInterp(lonIndexGrid,latIndexGrid) = 1
1716 fieldGrid_tmp(lonIndexGrid,latIndexGrid) = fieldGrid_tmp(lonIndexGrid,latIndexGrid)/real(ngp)
1717 nfill(myThreadNum) = nfill(myThreadNum) + 1
1718 end if
1719
1720 end do
1721 end do
1722 !$OMP END PARALLEL DO
1723
1724 fieldGrid(:,:) = fieldGrid_tmp(:,:)
1725 numFilledByAvg(:,:) = numFilledByAvg(:,:) + filledByInterp(:,:)
1726 end do boxSizeLoop
1727
1728 ! find any remaining grid points that are likely water and assign value
1729 if ( checkCloudToGridUnassigned .and. &
1730 stateVectorCloud%oceanMask%maskPresent .and. &
1731 .not. stateVectorGrid%oceanMask%maskPresent ) then
1732
1733 ! create a kdtree to index all cloud locations
1734 allocate(positionArray(3,niCloud*njCloud))
1735 gridIndex = 0
1736 do lonIndexCloud = 1, niCloud
1737 do latIndexCloud = 1, njCloud
1738 gridIndex = gridIndex + 1
1739 positionArray(:,gridIndex) = &
1740 kdtree2_3dPosition(real(stateVectorCloud%hco%lon2d_4(lonIndexCloud,latIndexCloud),8), &
1741 real(stateVectorCloud%hco%lat2d_4(lonIndexCloud,latIndexCloud),8))
1742 end do
1743 end do
1744 tree => kdtree2_create(positionArray, sort=.true., rearrange=.true.)
1745
1746 ! assign grid mask to value of mask at nearest cloud location
1747 do latIndexGrid = 1, njGrid
1748 do lonIndexGrid = 1, niGrid
1749 if (numFilledByAvg(lonIndexGrid,latIndexGrid) > 0) cycle
1750
1751 ! find nearest cloud location
1752 refPosition(:) = kdtree2_3dPosition(real(stateVectorGrid%hco%lon2d_4(lonIndexGrid,latIndexGrid),8), &
1753 real(stateVectorGrid%hco%lat2d_4(lonIndexGrid,latIndexGrid),8))
1754
1755 searchRadiusSquared = (50.0D0*1000.0D0)**2
1756 call kdtree2_r_nearest(tp=tree, qv=refPosition, r2=searchRadiusSquared, &
1757 nfound=numLocalGridPointsFound, &
1758 nalloc=maxNumLocalGridPointsSearch, results=searchResults)
1759 if (numLocalGridPointsFound > maxNumLocalGridPointsSearch) then
1760 call utl_abort('The parameter maxNumLocalGridPointsSearch must be increased')
1761 end if
1762 if (numLocalGridPointsFound == 0) then
1763 ! very far from a cloud location, assume it is land
1764 maskGrid(lonIndexGrid,latIndexGrid) = 0
1765 else
1766 gridIndex = searchResults(1)%idx
1767 lonIndexCloud = 1 + ((gridIndex-1)/njCloud)
1768 latIndexCloud = gridIndex - ((gridIndex-1)/njCloud)*njCloud
1769 if(lonIndexCloud < 1 .or. lonIndexCloud > niCloud) call utl_abort('lonIndexCloud wrong')
1770 if(latIndexCloud < 1 .or. latIndexCloud > njCloud) call utl_abort('latIndexCloud wrong')
1771 maskGrid(lonIndexGrid,latIndexGrid) = maskCloud(lonIndexCloud,latIndexCloud)
1772 end if
1773
1774 end do
1775 end do
1776
1777 deallocate(positionArray)
1778 call kdtree2_destroy(tree)
1779
1780 end if
1781
1782 ! fill in remaining grid points
1783 if (count(numFilledByAvg(:,:) > 0) > 0) then
1784
1785 ! compute spatial mean of assigned grid values
1786 extrapValue = 0.0d0
1787 do latIndexGrid = 1, njGrid
1788 do lonIndexGrid = 1, niGrid
1789 if (numFilledByAvg(lonIndexGrid,latIndexGrid) == 0) cycle
1790 extrapValue = extrapValue + real(fieldGrid(lonIndexGrid,latIndexGrid),8)
1791 end do
1792 end do
1793 extrapValue = extrapValue / real(count(numFilledByAvg(:,:) > 0),8)
1794
1795 ! set unassigned grid points and count them
1796 nextrap0 = 0
1797 nextrap1 = 0
1798 do latIndexGrid = 1, njGrid
1799 do lonIndexGrid = 1, niGrid
1800
1801 if (numFilledByAvg(lonIndexGrid,latIndexGrid) > 0) cycle
1802
1803 fieldGrid(lonIndexGrid,latIndexGrid) = extrapValue
1804
1805 if (maskGrid(lonIndexGrid,latIndexGrid) == 1) then
1806 nextrap1 = nextrap1 + 1
1807 else if (maskGrid(lonIndexGrid,latIndexGrid) == 0) then
1808 nextrap0 = nextrap0 + 1
1809 else
1810 call msg('int_sintCloudToGrid_gsv', &
1811 'Expecting 0 or 1 for the mask field at grid point: ' &
1812 //str(lonIndexGrid)//', '//str(latIndexGrid)//', ' &
1813 //str(maskGrid(lonIndexGrid,latIndexGrid)))
1814 call utl_abort('int_sintCloudToGrid_gsv')
1815 end if
1816
1817 end do
1818 end do
1819 end if
1820
1821 if ( checkCloudToGridUnassigned ) then
1822 call msg('int_sintCloudToGrid_gsv', &
1823 new_line('')//'Total number of grid points: '//str(niGrid*njGrid) &
1824 //new_line('')//'Number of grid points not covered by the cloud of points: '//str(sum(nhole(:))) &
1825 //new_line('')//'Number of grid points filled by neighbours: '//str(sum(nfill(:))) &
1826 //new_line('')//'Number of grid points with extrapolated value in masked area: '//str(nextrap0) &
1827 //new_line('')//'Number of grid points with extrapolated value in visible area: '//str(nextrap1))
1828
1829 if ( nextrap1 > 0 ) then
1830 call utl_abort('int_sintCloudToGrid_gsv: Values at some unmasked grid points were not assigned')
1831 end if
1832 end if
1833
1834 deallocate(fieldGrid_tmp)
1835 deallocate(xCloud)
1836 deallocate(yCloud)
1837 deallocate(numFilledByAvg)
1838 deallocate(filledByInterp)
1839 deallocate(maskCloud)
1840 deallocate(maskGrid)
1841
1842 ierr = 0
1843
1844 call msg('int_sintCloudToGrid_gsv', 'END', verb_opt=2)
1845 call utl_tmg_stop(176)
1846
1847 end function int_sintCloudToGrid_gsv
1848
1849 !--------------------------------------------------------------------------
1850 ! int_hInterpScalar_r8_2d
1851 !--------------------------------------------------------------------------
1852 function int_hInterpScalar_r8_2d(fieldOut_r8, fieldIn_r8, interpDegree, extrapDegree_opt) result(ierr)
1853 !
1854 !:Purpose: Horizontal interpolation of 2D scalar field that use real(8) arrays
1855 ! for input and output. Accessed through int_hInterpScalar.
1856 !
1857 implicit none
1858
1859 ! Arguments:
1860 real(8), intent(inout) :: fieldOut_r8(:,:)
1861 real(8), intent(in) :: fieldIn_r8(:,:)
1862 character(len=*), intent(in) :: interpDegree
1863 character(len=*), optional, intent(in) :: extrapDegree_opt
1864 ! Result:
1865 integer :: ierr
1866
1867 ! Locals:
1868 integer :: nii, nji, nio, njo
1869 integer :: jk1, jk2
1870 real(4), allocatable :: bufferi4(:,:), buffero4(:,:)
1871 integer :: ezsint
1872
1873 call msg('int_hInterpScalar_r8_2d', 'START', verb_opt=4)
1874 ! read the namelist
1875 call int_readNml()
1876
1877 ! do the standard interpolation
1878
1879 call int_setezopt(interpDegree, extrapDegree_opt)
1880
1881 nii = size(fieldIn_r8,1)
1882 nji = size(fieldIn_r8,2)
1883
1884 nio = size(fieldOut_r8,1)
1885 njo = size(fieldOut_r8,2)
1886
1887 allocate(bufferi4(nii,nji))
1888 allocate(buffero4(nio,njo))
1889
1890 do jk2 = 1,nji
1891 do jk1 = 1,nii
1892 bufferi4(jk1,jk2) = fieldIn_r8(jk1,jk2)
1893 end do
1894 end do
1895
1896 ierr = ezsint(buffero4,bufferi4)
1897
1898 do jk2 = 1,njo
1899 do jk1 = 1,nio
1900 fieldOut_r8(jk1,jk2) = buffero4(jk1,jk2)
1901 end do
1902 end do
1903
1904 deallocate(bufferi4)
1905 deallocate(buffero4)
1906
1907 call msg('int_hInterpScalar_r8_2d', 'END', verb_opt=4)
1908 end function int_hInterpScalar_r8_2d
1909
1910 !--------------------------------------------------------------------------
1911 ! int_hInterpUV_gsv
1912 !--------------------------------------------------------------------------
1913 function int_hInterpUV_gsv(stateVectorOut, stateVectorIn, varName, levIndex, stepIndex, &
1914 interpDegree, extrapDegree_opt) result(ierr)
1915 !
1916 !:Purpose: Horizontal interpolation of 2D vector field that use stateVector objects
1917 ! for input and output. Accessed through int_hInterpUV.
1918 !
1919 implicit none
1920
1921 ! Arguments:
1922 type(struct_gsv), intent(inout) :: stateVectorOut
1923 type(struct_gsv), intent(inout) :: stateVectorIn
1924 character(len=*), intent(in) :: varName
1925 integer , intent(in) :: levIndex
1926 integer , intent(in) :: stepIndex
1927 character(len=*), intent(in) :: interpDegree
1928 character(len=*), optional, intent(in) :: extrapDegree_opt
1929 ! Result:
1930 integer :: ierr
1931
1932 ! Locals:
1933 real(4), pointer :: UUout4(:,:,:,:), VVout4(:,:,:,:), UUin4(:,:,:,:), VVin4(:,:,:,:)
1934 real(8), pointer :: UUout8(:,:,:,:), VVout8(:,:,:,:), UUin8(:,:,:,:), VVin8(:,:,:,:)
1935 real(4), pointer :: UVout4(:,:,:), UVin4(:,:,:)
1936 real(8), pointer :: UVout8(:,:,:), UVin8(:,:,:)
1937 integer :: ezuvint, ezdefset
1938
1939 call msg('int_hInterpUV_gsv', 'START', verb_opt=4)
1940
1941 ! read the namelist
1942 call int_readNml()
1943
1944 ! check if special interpolation is required
1945 if (stateVectorIn%hco%initialized .and. stateVectorOut%hco%initialized) then
1946 if (stateVectorIn%hco%grtyp == 'Y') then
1947 call utl_abort('int_hInterpUV_gsv: cloudToGrid not implemented')
1948 end if
1949 end if
1950
1951 ! do the standard interpolation
1952
1953 ierr = ezdefset(stateVectorOut%hco%EZscintID, stateVectorIn%hco%EZscintID)
1954 call int_setezopt(interpDegree, extrapDegree_opt)
1955
1956 if ( gsv_getDataKind(stateVectorOut) == 4 .and. gsv_getDataKind(stateVectorIn) == 4) then
1957
1958 if (trim(varName) == 'BOTH') then
1959 call gsv_getField(stateVectorOut, UUout4, 'UU')
1960 call gsv_getField(stateVectorOut, VVout4, 'VV')
1961 call gsv_getField(stateVectorIn, UUin4, 'UU')
1962 call gsv_getField(stateVectorIn, VVin4, 'VV')
1963 ierr = ezuvint(UUout4(:,:,levIndex,stepIndex),VVout4(:,:,levIndex,stepIndex), &
1964 UUin4(:,:,levIndex,stepIndex), VVin4(:,:,levIndex,stepIndex))
1965 else if (trim(varName) == 'UU') then
1966 call gsv_getField (stateVectorIn, UUin4)
1967 call gsv_getField (stateVectorOut, UUout4)
1968 call gsv_getFieldUV(stateVectorIn, UVin4, levIndex)
1969 call gsv_getFieldUV(stateVectorOut, UVout4, levIndex)
1970 ierr = ezuvint(UUout4(:,:,levIndex,stepIndex),UVout4(:,:,stepIndex), &
1971 UUin4(:,:,levIndex,stepIndex), UVin4(:,:,stepIndex))
1972 else if (trim(varName) == 'VV') then
1973 call gsv_getField (stateVectorIn, VVin4)
1974 call gsv_getField (stateVectorOut, VVout4)
1975 call gsv_getFieldUV(stateVectorIn, UVin4, levIndex)
1976 call gsv_getFieldUV(stateVectorOut, UVout4, levIndex)
1977 ierr = ezuvint(UVout4(:,:,stepIndex),VVout4(:,:,levIndex,stepIndex), &
1978 UVin4(:,:,stepIndex), VVin4(:,:,levIndex,stepIndex))
1979 else
1980 call utl_abort('int_hInterpUV_gsv: unexpected varName: '//trim(varName))
1981 end if
1982
1983 else if ( gsv_getDataKind(stateVectorOut) == 8 .and. gsv_getDataKind(stateVectorIn) == 8) then
1984
1985 ! allocate real(4) buffers for copying to/from for interpolation
1986 allocate(UUin4(stateVectorIn%hco%ni,stateVectorIn%hco%nj,1,1))
1987 allocate(VVin4(stateVectorIn%hco%ni,stateVectorIn%hco%nj,1,1))
1988 allocate(UUout4(stateVectorOut%hco%ni,stateVectorOut%hco%nj,1,1))
1989 allocate(VVout4(stateVectorOut%hco%ni,stateVectorOut%hco%nj,1,1))
1990
1991 if (trim(varName) == 'BOTH') then
1992 call gsv_getField(stateVectorOut, UUout8, 'UU')
1993 call gsv_getField(stateVectorOut, VVout8, 'VV')
1994 call gsv_getField(stateVectorIn, UUin8, 'UU')
1995 call gsv_getField(stateVectorIn, VVin8, 'VV')
1996 UUin4(:,:,1,1) = UUin8(:,:,levIndex,stepIndex)
1997 VVin4(:,:,1,1) = VVin8(:,:,levIndex,stepIndex)
1998 ierr = ezuvint(UUout4(:,:,1,1),VVout4(:,:,1,1),UUin4(:,:,1,1),VVin4(:,:,1,1))
1999 UUout8(:,:,levIndex,stepIndex) = UUout4(:,:,1,1)
2000 VVout8(:,:,levIndex,stepIndex) = VVout4(:,:,1,1)
2001 else if (trim(varName) == 'UU') then
2002 call gsv_getField (stateVectorIn, UUin8)
2003 call gsv_getField (stateVectorOut, UUout8)
2004 call gsv_getFieldUV(stateVectorIn, UVin8, levIndex)
2005 call gsv_getFieldUV(stateVectorOut, UVout8, levIndex)
2006 UUin4(:,:,1,1) = UUin8(:,:,levIndex,stepIndex)
2007 VVin4(:,:,1,1) = UVin8(:,:,stepIndex)
2008 ierr = ezuvint(UUout4(:,:,1,1),VVout4(:,:,1,1),UUin4(:,:,1,1),VVin4(:,:,1,1))
2009 UUout8(:,:,levIndex,stepIndex) = UUout4(:,:,1,1)
2010 UVout8(:,:,stepIndex) = VVout4(:,:,1,1)
2011 else if (trim(varName) == 'VV') then
2012 call gsv_getField (stateVectorIn, VVin8)
2013 call gsv_getField (stateVectorOut, VVout8)
2014 call gsv_getFieldUV(stateVectorIn, UVin8, levIndex)
2015 call gsv_getFieldUV(stateVectorOut, UVout8, levIndex)
2016 UUin4(:,:,1,1) = UVin8(:,:,stepIndex)
2017 VVin4(:,:,1,1) = VVin8(:,:,levIndex,stepIndex)
2018 ierr = ezuvint(UUout4(:,:,1,1),VVout4(:,:,1,1),UUin4(:,:,1,1),VVin4(:,:,1,1))
2019 UVout8(:,:,stepIndex) = UUout4(:,:,1,1)
2020 VVout8(:,:,levIndex,stepIndex) = VVout4(:,:,1,1)
2021 else
2022 call utl_abort('int_hInterpUV_gsv: unexpected varName: '//trim(varName))
2023 end if
2024
2025 deallocate(UUin4,VVin4,UUout4,VVout4)
2026
2027 else
2028
2029 call utl_abort('int_hInterpUV_gsv: not implemented for mixed dataKind')
2030
2031 end if
2032
2033 call msg('int_hInterpUV_gsv', 'END', verb_opt=4)
2034 end function int_hInterpUV_gsv
2035
2036 !--------------------------------------------------------------------------
2037 ! int_hInterpUV_r4_2d
2038 !--------------------------------------------------------------------------
2039 function int_hInterpUV_r4_2d(uuout, vvout, uuin, vvin, interpDegree, extrapDegree_opt) result(ierr)
2040 !
2041 !:Purpose: Horizontal interpolation of 2D vector field that use real(4) arrays
2042 ! for input and output. Accessed through int_hInterpUV.
2043 !
2044 implicit none
2045
2046 ! Arguments:
2047 real(4), intent(inout) :: uuout(:,:)
2048 real(4), intent(inout) :: vvout(:,:)
2049 real(4), intent(in) :: uuin(:,:)
2050 real(4), intent(in) :: vvin(:,:)
2051 character(len=*), intent(in) :: interpDegree
2052 character(len=*), optional, intent(in) :: extrapDegree_opt
2053 ! Result:
2054 integer :: ierr
2055
2056 ! Locals:
2057 integer :: ezuvint
2058
2059 call msg('int_hInterpUV_r4_2d', 'START', verb_opt=4)
2060 ! read the namelist
2061 call int_readNml()
2062
2063 ! do the standard interpolation
2064 call int_setezopt(interpDegree, extrapDegree_opt)
2065 ierr = ezuvint(uuout, vvout, uuin, vvin)
2066
2067 call msg('int_hInterpUV_r4_2d', 'END', verb_opt=4)
2068 end function int_hInterpUV_r4_2d
2069
2070 !--------------------------------------------------------------------------
2071 ! int_hInterpUV_r8_2d
2072 !--------------------------------------------------------------------------
2073 function int_hInterpUV_r8_2d(uuout, vvout, uuin, vvin, interpDegree, extrapDegree_opt) result(ierr)
2074 !
2075 !:Purpose: Horizontal interpolation of 2D vector field that use real(8) arrays
2076 ! for input and output. Accessed through int_hInterpUV.
2077 !
2078 implicit none
2079
2080 ! Arguments:
2081 real(8), intent(inout) :: uuout(:,:)
2082 real(8), intent(inout) :: vvout(:,:)
2083 real(8), intent(in) :: uuin(:,:)
2084 real(8), intent(in) :: vvin(:,:)
2085 character(len=*), intent(in) :: interpDegree
2086 character(len=*), optional, intent(in) :: extrapDegree_opt
2087 ! Result:
2088 integer :: ierr
2089
2090 ! Locals:
2091 integer :: nio, njo, nii, nji
2092 integer :: jk1, jk2
2093 real, allocatable :: bufuuout4(:,:), bufvvout4(:,:)
2094 real, allocatable :: bufuuin4(:,:), bufvvin4(:,:)
2095 integer :: ezuvint
2096
2097 call msg('int_hInterpUV_r8_2d', 'START', verb_opt=4)
2098 ! read the namelist
2099 call int_readNml()
2100
2101 ! do the standard interpolation
2102
2103 call int_setezopt(interpDegree, extrapDegree_opt)
2104
2105 nii = size(uuin,1)
2106 nji = size(uuin,2)
2107
2108 nio = size(uuout,1)
2109 njo = size(uuout,2)
2110
2111 allocate(bufuuout4(nio,njo))
2112 allocate(bufvvout4(nio,njo))
2113 allocate(bufuuin4(nii,nji))
2114 allocate(bufvvin4(nii,nji))
2115
2116 do jk2 = 1,nji
2117 do jk1 = 1,nii
2118 bufuuin4(jk1,jk2) = uuin(jk1,jk2)
2119 bufvvin4(jk1,jk2) = vvin(jk1,jk2)
2120 end do
2121 end do
2122
2123 ierr = ezuvint(bufuuout4, bufvvout4, bufuuin4, bufvvin4)
2124
2125 do jk2 = 1,njo
2126 do jk1 = 1,nio
2127 uuout(jk1,jk2) = bufuuout4(jk1,jk2)
2128 vvout(jk1,jk2) = bufvvout4(jk1,jk2)
2129 end do
2130 end do
2131
2132 deallocate(bufuuin4)
2133 deallocate(bufvvin4)
2134 deallocate(bufuuout4)
2135 deallocate(bufvvout4)
2136
2137 call msg('int_hInterpUV_r8_2d', 'END', verb_opt=4)
2138 end function int_hInterpUV_r8_2d
2139
2140 !--------------------------------------------------------------------------
2141 ! int_ezgdef
2142 !--------------------------------------------------------------------------
2143 function int_ezgdef(ni, nj, grtyp, grtypref, ig1, ig2, ig3, ig4, ax, ay) result(vezgdef)
2144 !
2145 !:Purpose: Subroutine wrapper for rmnlib procedure ezgdef.
2146 !
2147 implicit none
2148
2149 ! Arguments:
2150 integer, intent(in) :: ni
2151 integer, intent(in) :: nj
2152 integer, intent(in) :: ig1
2153 integer, intent(in) :: ig2
2154 integer, intent(in) :: ig3
2155 integer, intent(in) :: ig4
2156 real(8), intent(in) :: ax(:)
2157 real(8), intent(in) :: ay(:)
2158 character(len=*), intent(in) :: grtyp
2159 character(len=*), intent(in) :: grtypref
2160 ! Result:
2161 integer :: vezgdef
2162
2163 ! Locals:
2164 integer :: ier2, jk, ilenx, ileny
2165 real(4), allocatable :: bufax4(:), bufay4(:)
2166 integer :: ezgdef
2167
2168 if (grtyp .eq. 'Y') then
2169 ilenx = max(1,ni*nj)
2170 ileny = ilenx
2171 else if (grtyp .eq. 'Z') then
2172 ilenx = max(1,ni)
2173 ileny = max(1,nj)
2174 else
2175 call utl_abort('VEZGDEF: Grid type not supported')
2176 end if
2177
2178 allocate(bufax4(ilenx))
2179 allocate(bufay4(ileny))
2180
2181 do jk = 1,ilenx
2182 bufax4(jk) = ax(jk)
2183 end do
2184 do jk = 1,ileny
2185 bufay4(jk) = ay(jk)
2186 end do
2187
2188 ier2 = ezgdef(ni, nj, grtyp, grtypref, ig1, ig2, ig3, ig4, &
2189 bufax4, bufay4)
2190
2191 deallocate(bufax4)
2192 deallocate(bufay4)
2193
2194 vezgdef = ier2
2195
2196 end function int_ezgdef
2197
2198 !--------------------------------------------------------------------------
2199 ! int_cxgaig
2200 !--------------------------------------------------------------------------
2201 subroutine int_cxgaig(grtyp, ig1, ig2, ig3, ig4, xlat0, xlon0, dlat, dlon)
2202 !
2203 !:Purpose: Subroutine wrapper for rmnlib procedure cxgaig.
2204 !
2205 implicit none
2206
2207 ! Arguments:
2208 integer, intent(in) :: ig1
2209 integer, intent(in) :: ig2
2210 integer, intent(in) :: ig3
2211 integer, intent(in) :: ig4
2212 real(8), intent(in) :: xlat0
2213 real(8), intent(in) :: xlon0
2214 real(8), intent(in) :: dlat
2215 real(8), intent(in) :: dlon
2216 character(len=*), intent(in) :: grtyp
2217
2218 ! Locals:
2219 real(4) :: xlat04, xlon04, dlat4, dlon4
2220
2221 xlat04 = xlat0
2222 xlon04 = xlon0
2223 dlat4 = dlat
2224 dlon4 = dlon
2225
2226 call cxgaig(grtyp, ig1, ig2, ig3, ig4, xlat04, xlon04, dlat4, dlon4)
2227
2228 end subroutine int_cxgaig
2229
2230end module interpolation_mod
