1module calcHeightAndPressure_mod
2 ! MODULE calcHeightAndPressure_mod (prefix='czp' category='4. Data Object transformations')
3 !
4 !:Purpose: Subroutines for computing height and/or pressure on statevectors
5 ! and columns depending on the vgrid kind.
6 ! Nonlinear, tangent-linear and adjoint versions of these
7 ! transformations are included in separate subroutines.
8 ! Depending on the vertical representation of the state or column,
9 ! pressure or height values are either computed or retrieved using
10 ! the vgrid (https://gitlab.science.gc.ca/RPN-SI/vgrid) library.
11 ! When computation is required (for instance to compute height on a
12 ! GEM-P, represented on pressure coordinates), thermodynamical
13 ! variables are required, typically `P0`, `TT` and `HU`.
14 ! Height and pressure values are obtained for both thermodynamical
15 ! and momentum levels and labeled `Z_T` (`P_T`) and `Z_M` (`P_M`).
16 !
17 use codePrecision_mod
18 use mathPhysConstants_mod
19 use earthConstants_mod
20 use physicsFunctions_mod
21 use verticalCoord_mod
22 use gridstatevector_mod
23 use columnData_mod
24 use utilities_mod
25 use message_mod
26 use gps_mod
27 use HorizontalCoord_mod
28 use Vgrid_Descriptors
29 implicit none
30 save
31 private
32
33 ! public procedures
34 public :: czp_calcZandP_nl, czp_calcZandP_tl, czp_calcZandP_ad
35 public :: czp_calcHeight_nl, czp_calcHeight_tl, czp_calcHeight_ad
36 public :: czp_calcPressure_nl, czp_calcPressure_tl, czp_calcPressure_ad
37 public :: czp_calcReturnHeight_gsv_nl, czp_calcReturnPressure_gsv_nl
38 public :: czp_calcReturnHeight_col_nl, czp_calcReturnPressure_col_nl
39 public :: czp_ensureCompatibleTops
40 public :: czp_fetch3DLevels, czp_fetch1DLevels, czp_fetch1DdPdPs
41
42 interface czp_fetch3DLevels
43 module procedure fetch3DLevels_r8
44 module procedure fetch3DLevels_r4
45 end interface czp_fetch3DLevels
46 interface czp_fetch1DLevels
47 module procedure fetch1DLevels_r8
48 end interface czp_fetch1DLevels
49 interface czp_fetch1DdPdPs
50 module procedure fetch1DdPdPs_r8
51 end interface czp_fetch1DdPdPs
52 interface czp_calcZandP_nl
53 module procedure calcZandP_gsv_nl
54 module procedure calcZandP_col_nl
55 end interface czp_calcZandP_nl
56 interface czp_calcZandP_tl
57 module procedure calcZandP_gsv_tl
58 module procedure calcZandP_col_tl
59 end interface czp_calcZandP_tl
60 interface czp_calcZandP_ad
61 module procedure calcZandP_gsv_ad
62 module procedure calcZandP_col_ad
63 end interface czp_calcZandP_ad
64
65 interface czp_calcHeight_nl
66 module procedure calcHeight_gsv_nl
67 module procedure calcHeight_col_nl
68 end interface czp_calcHeight_nl
69 interface czp_calcHeight_tl
70 module procedure calcHeight_gsv_tl
71 module procedure calcHeight_col_tl
72 end interface czp_calcHeight_tl
73 interface czp_calcHeight_ad
74 module procedure calcHeight_gsv_ad
75 module procedure calcHeight_col_ad
76 end interface czp_calcHeight_ad
77
78 interface czp_calcPressure_nl
79 module procedure calcPressure_gsv_nl
80 module procedure calcPressure_col_nl
81 end interface czp_calcPressure_nl
82 interface czp_calcPressure_tl
83 module procedure calcPressure_gsv_tl
84 module procedure calcPressure_col_tl
85 end interface czp_calcPressure_tl
86 interface czp_calcPressure_ad
87 module procedure calcPressure_gsv_ad
88 module procedure calcPressure_col_ad
89 end interface czp_calcPressure_ad
90
91 ! private module variables
92 real(8), allocatable :: coeff_M_TT_gsv(:,:,:,:), coeff_M_HU_gsv(:,:,:,:)
93 real(8), allocatable :: coeff_T_TT_gsv(:,:,:), coeff_T_HU_gsv(:,:,:)
94 real(8), allocatable :: coeff_M_P0_delPM_gsv(:,:,:,:)
95 real(8), allocatable :: coeff_M_P0_dP_delPT_gsv(:,:,:,:)
96 real(8), allocatable :: coeff_M_P0_dP_delP0_gsv(:,:,:,:)
97 real(8), allocatable :: coeff_T_P0_delP1_gsv(:,:,:)
98 real(8), allocatable :: coeff_T_P0_dP_delPT_gsv(:,:,:)
99 real(8), allocatable :: coeff_T_P0_dP_delP0_gsv(:,:,:)
100
101 real(8), allocatable :: coeff_M_TT_col(:,:), coeff_M_HU_col(:,:)
102 real(8), allocatable :: coeff_T_TT_col(:), coeff_T_HU_col(:)
103 real(8), allocatable :: coeff_M_P0_delPM_col(:,:)
104 real(8), allocatable :: coeff_M_P0_dP_delPT_col(:,:)
105 real(8), allocatable :: coeff_M_P0_dP_delP0_col(:,:)
106 real(8), allocatable :: coeff_T_P0_delP1_col(:), coeff_T_P0_dP_delPT_col(:)
107 real(8), allocatable :: coeff_T_P0_dP_delP0_col(:)
108
109contains
110 !---------------------------------------------------------------------
111 ! subroutines operating on struct_gsv
112 !---------------------------------------------------------------------
113
114 !---------------------------------------------------------
115 ! calcZandP_gsv_nl
116 !---------------------------------------------------------
117 subroutine calcZandP_gsv_nl(statevector)
118 !
119 !:Purpose: Pressure and height computation on the grid in proper order
120 ! depending on the vgrid kind.
121 ! Depending on the vcode, the routine will check the existence of
122 ! P_* (vcode=5xxx) or Z_* (vcode=2100x) first and proceed with
123 ! pressure (height) computation. Then, if the other variables are
124 ! also present, it will secondly compute height (pressure).
125 ! Hence if only P_* (Z_*) is present, only these are computed.
126 ! If the first variable P_* (Z_*) is not present, nothing is done.
127 !
128 implicit none
129
130 ! Arguments:
131 type(struct_gsv), intent(inout) :: statevector ! statevector that will contain the Z_*/P_* fields
132
133 ! Locals:
134 integer :: Vcode
135
136 call msg('calcZandP_gsv_nl (czp)', 'START', verb_opt=2)
137
138 Vcode = gsv_getVco(statevector)%vcode
139
140 if (Vcode == 5002 .or. Vcode == 5005 .or. Vcode == 5100) then
141 ! if P_T, P_M not allocated : do nothing
142 if (gsv_varExist(statevector, 'P_*')) then
143 call calcPressure_gsv_nl(statevector)
144 if (gsv_varExist(statevector, 'Z_*')) then
145 call calcHeight_gsv_nl(statevector)
146 end if
147 end if
148 else if (Vcode == 21001) then
149 ! if Z_T, Z_M not allocated : do nothing
150 if (gsv_varExist(statevector, 'Z_*')) then
151 call calcHeight_gsv_nl(statevector)
152 if (gsv_varExist(statevector, 'P_*')) then
153 call calcPressure_gsv_nl(statevector)
154 end if
155 end if
156 end if
157
158 call msg('calcZandP_gsv_nl (czp)', 'END', verb_opt=2)
159 end subroutine calcZandP_gsv_nl
160
161 !---------------------------------------------------------
162 ! calcZandP_gsv_tl
163 !---------------------------------------------------------
164 subroutine calcZandP_gsv_tl(statevector, statevectorRef)
165 !
166 !:Purpose: Pressure and height increment computation on the grid in proper
167 ! order depending on the vgrid kind.
168 ! Depending on the vcode, the routine will check the existence of
169 ! P_* (vcode=5xxx) or Z_* (vcode=2100x) first and proceed with
170 ! pressure (height) computation. Then, if the other variables are
171 ! also present, it will secondly compute height (pressure).
172 ! Hence if only P_* (Z_*) is present, only these are computed.
173 ! If the first variable P_* (Z_*) is not present, nothing is done.
174 !
175 implicit none
176
177 ! Arguments:
178 type(struct_gsv), intent(inout) :: statevector ! statevector that will contain the Z_*/P_* increments
179 type(struct_gsv), intent(in) :: statevectorRef ! statevector containing needed reference fields
180
181 ! Locals:
182 type(struct_vco), pointer :: vco
183 integer :: Vcode
184
185 call msg('calcZandP_gsv_tl (czp)', 'START', verb_opt=2)
186
187 vco => gsv_getVco(statevector)
188 Vcode = vco%vcode
189
190 if (Vcode == 0) return
191
192 if (Vcode == 5002 .or. Vcode == 5005) then
193 ! if P_T, P_M not allocated : do nothing
194 if (gsv_varExist(statevector, 'P_*')) then
195
196 if ( .not. gsv_containsNonZeroValues(stateVectorRef) ) then
197 call utl_abort('calcZandP_gsv_tl: stateVectorRef not initialized')
198 end if
199 call calcPressure_gsv_tl(statevector, statevectorRef)
200
201 if (gsv_varExist(statevector, 'Z_*')) then
202 call calcHeight_gsv_tl(statevector, statevectorRef)
203 end if
204
205 end if
206 else if (Vcode == 21001) then
207 ! if Z_T, Z_M not allocated : do nothing
208 if (gsv_varExist(statevector, 'Z_*')) then
209
210 if ( .not. gsv_containsNonZeroValues(stateVectorRef) ) then
211 call utl_abort('calcZandP_gsv_tl: stateVectorRef not initialized')
212 end if
213 call calcHeight_gsv_tl(statevector, statevectorRef)
214
215 if (gsv_varExist(statevector, 'P_*')) then
216 call calcPressure_gsv_tl(statevector, statevectorRef)
217 end if
218
219 end if
220 else
221 call utl_abort('calcZandP_gsv_tl (czp): not implemented')
222 end if
223
224 call msg('calcZandP_gsv_tl (czp)', 'END', verb_opt=2)
225 end subroutine calcZandP_gsv_tl
226
227 !---------------------------------------------------------
228 ! calcZandP_gsv_ad
229 !---------------------------------------------------------
230 subroutine calcZandP_gsv_ad(statevector, statevectorRef)
231 !
232 !:Purpose: Pressure and height increment adjoint computation on the grid
233 ! in proper order depending on the vgrid kind
234 ! Depending on the vcode, the routine will check the existence of
235 ! Z_* (vcode=5xxx) or P_* (vcode=2100x) first and proceed with
236 ! height (pressure) adjoint computation. Then, if the other
237 ! variables are also present, it will secondly proceed with
238 ! adjoint computation of pressure (height).
239 ! Hence if only Z_* (P_*) is present, only these are computed.
240 ! If the first variable Z_* (P_*) is not present, nothing is done.
241 !
242 implicit none
243
244 ! Arguments:
245 type(struct_gsv), intent(inout) :: statevector ! statevector that will contain the Z_*/P_* increments
246 type(struct_gsv), intent(in) :: statevectorRef ! statevector containing needed reference fields
247
248 ! Locals:
249 type(struct_vco), pointer :: vco
250 integer :: Vcode
251
252 call msg('calcZandP_gsv_ad (czp)', 'START', verb_opt=2)
253
254 vco => gsv_getVco(statevector)
255 Vcode = vco%vcode
256
257 if (Vcode == 0) return
258
259 if (Vcode == 5002 .or. Vcode == 5005) then
260 ! if Z_T, Z_M not allocated : do nothing
261 if (gsv_varExist(statevector, 'Z_*')) then
262
263 if ( .not. gsv_containsNonZeroValues(stateVectorRef) ) then
264 call utl_abort('calcZandP_gsv_ad: stateVectorRef not initialized')
265 end if
266 call calcHeight_gsv_ad(statevector, statevectorRef)
267
268 if (gsv_varExist(statevector, 'P_*')) then
269 call calcPressure_gsv_ad(statevector, statevectorRef)
270 end if
271
272 end if
273 else if (Vcode == 21001) then
274 ! if P_T, P_M not allocated : do nothing
275 if (gsv_varExist(statevector, 'P_*')) then
276
277 if ( .not. gsv_containsNonZeroValues(stateVectorRef) ) then
278 call utl_abort('calcZandP_gsv_ad: stateVectorRef not initialized')
279 end if
280 call calcPressure_gsv_ad(statevector, statevectorRef)
281
282 if (gsv_varExist(statevector, 'Z_*')) then
283 call calcHeight_gsv_ad(statevector, statevectorRef)
284 end if
285
286 end if
287 else
288 call utl_abort('calcZandP_gsv_ad (czp): not implemented')
289 end if
290
291 call msg('calcZandP_gsv_ad (czp)', 'END', verb_opt=2)
292 end subroutine calcZandP_gsv_ad
293
294 !---------------------------------------------------------
295 ! calcHeight_gsv_nl
296 !---------------------------------------------------------
297 subroutine calcHeight_gsv_nl(statevector)
298 !
299 !:Purpose: Compute or retrieve heights and store values in statevector.
300 !
301 implicit none
302
303 ! Arguments:
304 type(struct_gsv), intent(inout) :: statevector
305
306 ! Locals:
307 integer :: Vcode
308 real(4), pointer :: ptr_PT_r4(:,:,:,:), ptr_PM_r4(:,:,:,:)
309 real(8), pointer :: ptr_PT_r8(:,:,:,:), ptr_PM_r8(:,:,:,:)
310 real(4), pointer :: ptr_ZT_r4(:,:,:,:), ptr_ZM_r4(:,:,:,:)
311 real(8), pointer :: ptr_ZT_r8(:,:,:,:), ptr_ZM_r8(:,:,:,:)
312
313 call utl_tmg_start(172,'low-level--czp_calcHeight_nl')
314 call msg('calcHeight_gsv_nl (czp)', 'START', verb_opt=2)
315
316 Vcode = gsv_getVco(statevector)%vcode
317 if (Vcode == 5005 .or. Vcode == 5002 .or. Vcode == 5100) then
318 if ( gsv_getDataKind(statevector) == 4 ) then
319 call gsv_getField(statevector, ptr_PT_r4, 'P_T')
320 call gsv_getField(statevector, ptr_PM_r4, 'P_M')
321 call gsv_getField(statevector, ptr_ZT_r4, 'Z_T')
322 call gsv_getField(statevector, ptr_ZM_r4, 'Z_M')
323 call calcHeight_gsv_nl_vcode5xxx( statevector, &
324 PTin_r4_opt=ptr_PT_r4, &
325 PMin_r4_opt=ptr_PM_r4, &
326 ZTout_r4_opt=ptr_ZT_r4, &
327 ZMout_r4_opt=ptr_ZM_r4)
328 else
329 call gsv_getField(statevector, ptr_PT_r8, 'P_T')
330 call gsv_getField(statevector, ptr_PM_r8, 'P_M')
331 call gsv_getField(statevector, ptr_ZT_r8, 'Z_T')
332 call gsv_getField(statevector, ptr_ZM_r8, 'Z_M')
333 call calcHeight_gsv_nl_vcode5xxx( statevector, &
334 PTin_r8_opt=ptr_PT_r8, &
335 PMin_r8_opt=ptr_PM_r8, &
336 ZTout_r8_opt=ptr_ZT_r8, &
337 ZMout_r8_opt=ptr_ZM_r8)
338 end if
339
340 else if (Vcode == 21001) then
341 ! Development notes (@mad001)
342 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
343 if ( gsv_getDataKind(statevector) == 4 ) then
344 call gsv_getField(statevector, ptr_ZT_r4, 'Z_T')
345 call gsv_getField(statevector, ptr_ZM_r4, 'Z_M')
346 call calcHeight_gsv_nl_vcode2100x_r4(statevector, ptr_ZT_r4, ptr_ZM_r4)
347 else
348 call gsv_getField(statevector, ptr_ZT_r8, 'Z_T')
349 call gsv_getField(statevector, ptr_ZM_r8, 'Z_M')
350 call calcHeight_gsv_nl_vcode2100x_r8(statevector, ptr_ZT_r8, ptr_ZM_r8)
351 end if
352 end if
353
354 if ( gsv_getDataKind(statevector) == 4 ) then
355 call msg('calcHeight_gsv_nl (czp)', &
356 new_line('')//'Z_M = '&
357 //str(ptr_ZM_r4(statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.) &
358 //new_line('')//'Z_T = '&
359 //str(ptr_ZT_r4( statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.), &
360 verb_opt=2)
361 else
362 call msg('calcHeight_gsv_nl (czp)', &
363 new_line('')//'Z_M = '&
364 //str(ptr_ZM_r8(statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.) &
365 //new_line('')//'Z_T = '&
366 //str(ptr_ZT_r8( statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.), &
367 verb_opt=2)
368 end if
369
370 call msg('calcHeight_gsv_nl (czp)', 'END', verb_opt=2)
371 call utl_tmg_stop(172)
372 end subroutine calcHeight_gsv_nl
373
374 !---------------------------------------------------------
375 ! czp_calcReturnHeight_gsv_nl
376 !---------------------------------------------------------
377 subroutine czp_calcReturnHeight_gsv_nl( statevector, &
378 PTin_r4_opt, PMin_r4_opt, &
379 PTin_r8_opt, PMin_r8_opt, &
380 ZTout_r4_opt, ZMout_r4_opt, &
381 ZTout_r8_opt, ZMout_r8_opt)
382 !
383 !:Purpose: Compute or retrieve heights and return values in pointer arguments.
384 ! Proceeds to vcode dispatching.
385 !
386 implicit none
387
388 ! Arguments:
389 type(struct_gsv), intent(in) :: statevector
390 real(4), optional, pointer, intent(in) :: PTin_r4_opt(:,:,:,:)
391 real(4), optional, pointer, intent(in) :: PMin_r4_opt(:,:,:,:)
392 real(8), optional, pointer, intent(in) :: PTin_r8_opt(:,:,:,:)
393 real(8), optional, pointer, intent(in) :: PMin_r8_opt(:,:,:,:)
394 real(4), optional, pointer, intent(inout) :: ZTout_r4_opt(:,:,:,:)
395 real(4), optional, pointer, intent(inout) :: ZMout_r4_opt(:,:,:,:)
396 real(8), optional, pointer, intent(inout) :: ZTout_r8_opt(:,:,:,:)
397 real(8), optional, pointer, intent(inout) :: ZMout_r8_opt(:,:,:,:)
398
399 ! Locals:
400 integer :: Vcode
401
402 call utl_tmg_start(172,'low-level--czp_calcHeight_nl')
403 call msg('czp_calcReturnHeight_gsv_nl', 'START', verb_opt=2)
404
405 Vcode = gsv_getVco(statevector)%vcode
406 if (Vcode == 5005 .or. Vcode == 5002 .or. Vcode == 5100) then
407 if ( gsv_getDataKind(statevector) == 4 ) then
408 if ( .not. (present(PTin_r4_opt) .and. present(PMin_r4_opt))) then
409 call utl_abort('czp_calcReturnHeight_gsv_nl: dataKind=4: P{T,M}out_r4_opt expected')
410 end if
411 if ( .not. (present(ZTout_r4_opt) .and. present(ZMout_r4_opt))) then
412 call utl_abort('czp_calcReturnHeight_gsv_nl: dataKind=4: Z{T,M}out_r4_opt expected')
413 end if
414 call calcHeight_gsv_nl_vcode5xxx( statevector, &
415 PTin_r4_opt=PTin_r4_opt, &
416 PMin_r4_opt=PMin_r4_opt, &
417 ZTout_r4_opt=ZTout_r4_opt, &
418 ZMout_r4_opt=ZMout_r4_opt)
419 else ! datakind = 8
420 if ( .not. (present(PTin_r8_opt) .and. present(PMin_r8_opt))) then
421 call utl_abort('czp_calcReturnHeight_gsv_nl: dataKind=8: P{T,M}out_r8_opt expected')
422 end if
423 if ( .not. (present(ZTout_r8_opt) .and. present(ZMout_r8_opt))) then
424 call utl_abort('czp_calcReturnHeight_gsv_nl: dataKind=8: Z{T,M}out_r8_opt expected')
425 end if
426 call calcHeight_gsv_nl_vcode5xxx( statevector, &
427 PTin_r8_opt=PTin_r8_opt, &
428 PMin_r8_opt=PMin_r8_opt, &
429 ZTout_r8_opt=ZTout_r8_opt, &
430 ZMout_r8_opt=ZMout_r8_opt)
431 end if
432
433 else if (Vcode == 21001) then
434 if ( gsv_getDataKind(statevector) == 4 ) then
435 if ( .not. (present(ZTout_r4_opt) .and. present(ZMout_r4_opt))) then
436 call utl_abort('czp_calcReturnHeight_gsv_nl: dataKind=4: Z{T,M}_r4 expected')
437 end if
438 call calcHeight_gsv_nl_vcode2100x_r4(statevector, ZTout_r4_opt, ZMout_r4_opt)
439 else
440 if ( .not. (present(ZTout_r8_opt) .and. present(ZMout_r8_opt))) then
441 call utl_abort('czp_calcReturnHeight_gsv_nl: dataKind=4: Z{T,M}_r4 expected')
442 end if
443 call calcHeight_gsv_nl_vcode2100x_r8(statevector, ZTout_r8_opt, ZMout_r8_opt)
444 end if
445 end if
446
447 call msg('czp_calcReturnHeight_gsv_nl', 'END', verb_opt=2)
448 call utl_tmg_stop(172)
449 end subroutine czp_calcReturnHeight_gsv_nl
450
451 !---------------------------------------------------------
452 ! calcHeight_gsv_nl_vcode2100x_r4
453 !---------------------------------------------------------
454 subroutine calcHeight_gsv_nl_vcode2100x_r4(statevector, Z_T, Z_M)
455 !
456 !:Purpose: Retrieve heights for GEM-H statevector, return height values
457 ! in pointer arguments.
458 ! real(4) version
459 !
460 implicit none
461
462 ! Arguments:
463 type(struct_gsv), intent(in) :: statevector
464 real(4), pointer, intent(inout) :: Z_T(:,:,:,:)
465 real(4), pointer, intent(inout) :: Z_M(:,:,:,:)
466
467 ! Locals:
468 integer :: numStep, stepIndex
469 real(kind=8), pointer :: Hsfc(:,:)
470 real(kind=4), allocatable :: Hsfc4(:,:)
471 real(kind=4), pointer :: GZHeightM_out(:,:,:), GZHeightT_out(:,:,:)
472
473 call msg('calcHeight_gsv_nl_vcode2100x_r4 (czp)', 'START', verb_opt=4)
474
475 if ( .not. gsv_varExist(statevector,'Z_*')) then
476 call utl_abort('calcHeight_gsv_nl_vcode2100x_r4 (czp): Z_T/Z_M do not exist in statevector!')
477 end if
478
479 allocate(Hsfc4( statevector%myLonBeg:statevector%myLonEnd, &
480 statevector%myLatBeg:statevector%myLatEnd))
481 Hsfc => gsv_getHeightSfc(statevector)
482 Hsfc4 = real(Hsfc,4)
483
484 numStep = statevector%numStep
485
486 do stepIndex = 1, numStep
487
488 call fetch3DLevels_r4(statevector%vco, Hsfc4, &
489 fldM_opt=GZHeightM_out, fldT_opt=GZHeightT_out)
490 Z_M(:,:,:,stepIndex) = gz2alt_r4(statevector, GZHeightM_out)
491 Z_T(:,:,:,stepIndex) = gz2alt_r4(statevector, GZHeightT_out)
492 deallocate(GZHeightM_out, GZHeightT_out)
493
494 end do
495 deallocate(Hsfc4)
496
497 call msg('calcHeight_gsv_nl_vcode2100x_r4 (czp)', 'END', verb_opt=4)
498 end subroutine calcHeight_gsv_nl_vcode2100x_r4
499
500 !---------------------------------------------------------
501 ! gz2alt_r4
502 !---------------------------------------------------------
503 function gz2alt_r4(statevector, gzHeight) result(alt)
504 !
505 !:Purpose: Iterative conversion of geopotential height to geometric
506 ! altitude. (solution proposed by J. Aparicio)
507 ! real(4) version.
508 !
509 implicit none
510
511 ! Arguments:
512 type(struct_gsv), intent(in) :: statevector
513 real(kind=4), pointer, intent(in) :: gzHeight(:,:,:)
514 ! Result:
515 real(kind=4), allocatable :: alt(:,:,:)
516
517 ! Locals:
518 integer :: nLon, nLat, nLev
519 type(struct_hco), pointer :: hco
520 real(kind=4) :: latitude
521 real(kind=4) :: gzH, b1, b2, A2, A3
522 integer :: lonIndex, latIndex, lvlIndex
523
524 ! gzHeight comes from external `vgd_levels` which does not know the
525 ! mpi shifted indexes
526 nLon = ubound(gzHeight, 1)
527 nLat = ubound(gzHeight, 2)
528 nLev = ubound(gzHeight, 3)
529 allocate(alt(nLon, nLat, nLev))
530
531 hco => gsv_getHco(statevector)
532
533 do lonIndex = 1, nLon
534 do latIndex = 1, nLat
535 do lvlIndex = 1, nLev
536 ! explicit shift of indexes
537 latitude = hco%lat2d_4( lonIndex+statevector%myLonBeg-1,&
538 latIndex+statevector%myLatBeg-1)
539 gzH = gzHeight(lonIndex, latIndex, lvlIndex)
540 ! gzH(alt) = g0 * (1 + b1*alt + b2*alt**2)
541 b1 = -2.0/ec_wgs_a*(1.0+ec_wgs_f+ec_wgs_m-2*ec_wgs_f*latitude**2)
542 b2 = 3.0/ec_wgs_a**2
543 ! reversed series coefficients (Abramowitz and Stegun 3.6.25)
544 A2 = -b1/2.0
545 A3 = b1**2/2.0 - b2/3.0
546 alt(lonIndex, latIndex, lvlIndex) = gzH + A2*gzH**2 + A3*gzH**3
547 end do
548 end do
549 end do
550
551 end function gz2alt_r4
552
553 !---------------------------------------------------------
554 ! calcHeight_gsv_nl_vcode2100x_r8
555 !---------------------------------------------------------
556 subroutine calcHeight_gsv_nl_vcode2100x_r8(statevector, Z_T, Z_M)
557 !
558 !:Purpose: Retrieve heights for GEM-H statevector, return height values
559 ! in pointer arguments.
560 ! real(8) version
561 !
562 implicit none
563
564 ! Arguments:
565 type(struct_gsv), intent(in) :: statevector
566 real(8), pointer, intent(inout) :: Z_T(:,:,:,:)
567 real(8), pointer, intent(inout) :: Z_M(:,:,:,:)
568
569 ! Locals:
570 integer :: numStep, stepIndex
571 real(kind=8), pointer :: Hsfc(:,:), GZHeightM_out(:,:,:), GZHeightT_out(:,:,:)
572
573 call msg('calcHeight_gsv_nl_vcode2100x_r8 (czp)', 'START', verb_opt=4)
574
575 Hsfc => gsv_getHeightSfc(statevector)
576 numStep = statevector%numStep
577
578 do stepIndex = 1, numStep
579
580 call fetch3DLevels_r8(statevector%vco, Hsfc, &
581 fldM_opt=GZHeightM_out, fldT_opt=GZHeightT_out)
582 Z_M(:,:,:,stepIndex) = gz2alt_r8(statevector, GZHeightM_out)
583 Z_T(:,:,:,stepIndex) = gz2alt_r8(statevector, GZHeightT_out)
584 deallocate(GZHeightM_out, GZHeightT_out)
585 end do
586
587 call msg('calcHeight_gsv_nl_vcode2100x_r8 (czp)', 'END', verb_opt=4)
588 end subroutine calcHeight_gsv_nl_vcode2100x_r8
589
590 !---------------------------------------------------------
591 ! gz2alt_r8
592 !---------------------------------------------------------
593 function gz2alt_r8(statevector, gzHeight) result(alt)
594 !
595 !:Purpose: Iterative conversion of geopotential height to geometric
596 ! altitude. (solution proposed by J. Aparicio)
597 ! real(8) version.
598 !
599 implicit none
600
601 ! Arguments:
602 type(struct_gsv), intent(in) :: statevector
603 real(kind=8), pointer, intent(in) :: gzHeight(:,:,:)
604 ! Result:
605 real(kind=8), allocatable :: alt(:,:,:)
606
607 ! Locals:
608 integer :: nLon, nLat, nLev
609 type(struct_hco), pointer :: hco
610 real(kind=8) :: latitude
611 real(kind=8) :: gzH, b1, b2, A2, A3
612 integer :: lonIndex, latIndex, lvlIndex
613
614 ! gzHeight comes from external `vgd_levels` which does not know the
615 ! mpi shifted indexes
616 nLon = ubound(gzHeight, 1)
617 nLat = ubound(gzHeight, 2)
618 nLev = ubound(gzHeight, 3)
619 allocate(alt(nLon, nLat, nLev))
620
621 hco => gsv_getHco(statevector)
622
623 do lonIndex = 1, nLon
624 do latIndex = 1, nLat
625 do lvlIndex = 1, nLev
626 ! explicit shift of indexes
627 latitude = hco%lat2d_4( lonIndex+statevector%myLonBeg-1,&
628 latIndex+statevector%myLatBeg-1)
629 gzH = gzHeight(lonIndex, latIndex, lvlIndex)
630 ! gzH(alt) = g0 * (1 + b1*alt + b2*alt**2)
631 b1 = -2.0D0/ec_wgs_a*(1.0D0+ec_wgs_f+ec_wgs_m-2*ec_wgs_f*latitude**2)
632 b2 = 3.0D0/ec_wgs_a**2
633 ! reversed series coefficients (Abramowitz and Stegun 3.6.25)
634 A2 = -b1/2.0D0
635 A3 = b1**2/2.0D0 - b2/3.0D0
636 alt(lonIndex, latIndex, lvlIndex) = gzH + A2*gzH**2 + A3*gzH**3
637 end do
638 end do
639 end do
640
641 end function gz2alt_r8
642
643 !---------------------------------------------------------
644 ! calcHeight_gsv_nl_vcode5xxx
645 !---------------------------------------------------------
646 subroutine calcHeight_gsv_nl_vcode5xxx( statevector, &
647 PTin_r4_opt, PMin_r4_opt, &
648 PTin_r8_opt, PMin_r8_opt, &
649 ZTout_r4_opt, ZMout_r4_opt, &
650 ZTout_r8_opt, ZMout_r8_opt)
651 !
652 !:Purpose: Compute heights for GEM-P statevector, return height values
653 ! in pointer arguments.
654 ! Assumptions:
655 ! 1) nlev_T = nlev_M+1 (for vcode=5002)
656 ! 2) alt_T(nlev_T) = alt_M(nlev_M), both at the surface
657 ! 3) a thermo level exists at the top, higher than the highest
658 ! momentum level
659 ! 4) the placement of the thermo levels means that alt_T is the
660 ! average of 2 nearest alt_M (according to Ron and Claude)
661 !
662 implicit none
663
664 ! Arguments:
665 type(struct_gsv), intent(in) :: statevector
666 real(4), pointer, optional, intent(in) :: PTin_r4_opt(:,:,:,:)
667 real(4), pointer, optional, intent(in) :: PMin_r4_opt(:,:,:,:)
668 real(8), pointer, optional, intent(in) :: PTin_r8_opt(:,:,:,:)
669 real(8), pointer, optional, intent(in) :: PMin_r8_opt(:,:,:,:)
670 real(4), pointer, optional, intent(inout) :: ZTout_r4_opt(:,:,:,:)
671 real(4), pointer, optional, intent(inout) :: ZMout_r4_opt(:,:,:,:)
672 real(8), pointer, optional, intent(inout) :: ZTout_r8_opt(:,:,:,:)
673 real(8), pointer, optional, intent(inout) :: ZMout_r8_opt(:,:,:,:)
674
675 ! Locals:
676 integer :: lev_M,lev_T,nlev_M,nlev_T,status,Vcode
677 integer :: numStep, stepIndex, latIndex,lonIndex
678 real(4) :: lat_4, heightSfcOffset_T_r4, heightSfcOffset_M_r4
679 real(8) :: delThick, ratioP
680 real(8) :: ScaleFactorBottom, ScaleFactorTop
681 real(8) :: P_M, P_M1, P_Mm1, P_T
682 real(8) :: hu, tt, Pr, cmp, h0, Rgh, P0, dh, rMt
683 real(8) :: sLat, cLat, lat_8
684 real(8), allocatable :: tv(:), height_T(:), height_M(:)
685 real(4), pointer :: height_T_ptr_r4(:,:,:,:)
686 real(4), pointer :: height_M_ptr_r4(:,:,:,:)
687 real(4), pointer :: hu_ptr_r4(:,:,:,:),tt_ptr_r4(:,:,:,:)
688 real(4), pointer :: P_T_ptr_r4(:,:,:,:),P_M_ptr_r4(:,:,:,:)
689 real(4), pointer :: P0_ptr_r4(:,:,:,:)
690 real(8), pointer :: height_T_ptr_r8(:,:,:,:)
691 real(8), pointer :: height_M_ptr_r8(:,:,:,:)
692 real(8), pointer :: P_T_ptr_r8(:,:,:,:),P_M_ptr_r8(:,:,:,:)
693 real(8), pointer :: P0_ptr_r8(:,:,:,:)
694 real(8), pointer :: hu_ptr_r8(:,:,:,:),tt_ptr_r8(:,:,:,:)
695 real(8), pointer :: HeightSfc_ptr_r8(:,:)
696
697 call msg('calcHeight_gsv_nl_vcode5xxx (czp)', 'START', verb_opt=4)
698
699 nlev_T = gsv_getNumLev(statevector,'TH')
700 nlev_M = gsv_getNumLev(statevector,'MM')
701 Vcode = gsv_getVco(statevector)%vcode
702 numStep = statevector%numStep
703
704 allocate(tv(nlev_T))
705
706 if (Vcode == 5002 .and. nlev_T /= nlev_M+1) then
707 call utl_abort('calcHeight_gsv_nl_vcode5xxx (czp): nlev_T is not equal to nlev_M+1!')
708 end if
709 if ((Vcode == 5005 .or. Vcode == 5100) .and. nlev_T /= nlev_M) then
710 call utl_abort('calcHeight_gsv_nl_vcode5xxx (czp): nlev_T is not equal to nlev_M!')
711 end if
712
713 if (Vcode == 5005 .or. Vcode == 5100) then
714 status = vgd_get( statevector%vco%vgrid, &
715 key='DHM - height of the diagnostic level (m)', &
716 value=heightSfcOffset_M_r4)
717 status = vgd_get( statevector%vco%vgrid, &
718 key='DHT - height of the diagnostic level (t)', &
719 value=heightSfcOffset_T_r4)
720 call msg('calcHeight_gsv_nl_vcode5xxx (czp)', &
721 'height offset for near-sfc momentum level is:'//str(heightSfcOffset_M_r4)//' meters'&
722 //new_line('')//'height offset for near-sfc thermo level is:'//str(heightSfcOffset_T_r4)//' meters', &
723 verb_opt=2, mpiAll_opt=.false.)
724 if ( .not.statevector%addHeightSfcOffset ) then
725 call msg('calcHeight_gsv_nl_vcode5xxx (czp)', new_line('') &
726 //'--------------------------------------------------------------------------'//new_line('')&
727 //'BUT HEIGHT OFFSET REMOVED FOR DIAGNOSTIC LEVELS FOR BACKWARD COMPATIBILITY'//new_line('')&
728 //'--------------------------------------------------------------------------', &
729 verb_opt=2, mpiAll_opt=.false.)
730 end if
731 end if
732
733 allocate(height_T(nlev_T))
734 allocate(height_M(nlev_M))
735
736 if ( statevector%dataKind == 4 ) then
737 if ( .not. (present(PTin_r4_opt) .and. present(PMin_r4_opt))) then
738 call utl_abort('calcHeight_gsv_nl_vcode5xxx (czp): dataKind=4: P{T,M}in_r4_opt expected')
739 end if
740 P_T_ptr_r4 => PTin_r4_opt
741 P_M_ptr_r4 => PMin_r4_opt
742
743 if ( .not. (present(ZTout_r4_opt) .and. present(ZMout_r4_opt))) then
744 call utl_abort('calcHeight_gsv_nl_vcode5xxx (czp): dataKind=4: Z{T,M}out_r4_opt expected')
745 end if
746 height_M_ptr_r4 => ZMout_r4_opt
747 height_T_ptr_r4 => ZTout_r4_opt
748
749 ! initialize the height pointer to zero
750 height_M_ptr_r4(:,:,:,:) = 0.0
751 height_T_ptr_r4(:,:,:,:) = 0.0
752
753 call gsv_getField(statevector,hu_ptr_r4,'HU')
754 call gsv_getField(statevector,tt_ptr_r4,'TT')
755 call gsv_getField(statevector,P0_ptr_r4,'P0')
756
757 else ! datakind = 8
758 if ( .not. (present(PTin_r8_opt) .and. present(PMin_r8_opt))) then
759 call utl_abort('calcHeight_gsv_nl_vcode5xxx (czp): dataKind=8: P{T,M}in_r8_opt expected')
760 end if
761 P_T_ptr_r8 => PTin_r8_opt
762 P_M_ptr_r8 => PMin_r8_opt
763
764 if ( .not. (present(ZTout_r8_opt) .and. present(ZMout_r8_opt))) then
765 call utl_abort('calcHeight_gsv_nl_vcode5xxx (czp): dataKind=8: Z{T,M}out_r8_opt expected')
766 end if
767 height_M_ptr_r8 => ZMout_r8_opt
768 height_T_ptr_r8 => ZTout_r8_opt
769
770 ! initialize the height pointer to zero
771 height_M_ptr_r8(:,:,:,:) = 0.0d0
772 height_T_ptr_r8(:,:,:,:) = 0.0d0
773
774 call gsv_getField(statevector,hu_ptr_r8,'HU')
775 call gsv_getField(statevector,tt_ptr_r8,'TT')
776 call gsv_getField(statevector,P0_ptr_r8,'P0')
777 end if
778
779 HeightSfc_ptr_r8 => gsv_getHeightSfc(statevector)
780
781 ! compute virtual temperature on thermo levels (corrected of compressibility)
782 do_computeHeight_gsv_nl : do stepIndex = 1, numStep
783 !$OMP PARALLEL DO PRIVATE(latIndex, lonIndex, height_T, height_M, lat_4, lat_8, sLat, cLat, lev_T, &
784 !$OMP hu, tt, Pr, cmp, tv, rMT, P_T, P_M, P0, ratioP, h0, Rgh, dh, delThick, lev_M, P_M1, &
785 !$OMP scaleFactorBottom, scaleFactorTop, P_Mm1)
786 do latIndex = statevector%myLatBeg, statevector%myLatEnd
787 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
788
789 height_T(:) = 0.0D0
790 height_M(:) = 0.0D0
791
792 ! latitude
793 lat_4 = statevector%hco%lat2d_4(lonIndex,latIndex)
794 lat_8 = real(lat_4,8)
795 sLat = sin(lat_8)
796 cLat = cos(lat_8)
797
798 do lev_T = 1, nlev_T
799 if ( statevector%dataKind == 4 ) then
800 hu = real(hu_ptr_r4(lonIndex,latIndex,lev_T,stepIndex),8)
801 tt = real(tt_ptr_r4(lonIndex,latIndex,lev_T,stepIndex),8)
802 Pr = real(P_T_ptr_r4(lonIndex,latIndex,lev_T,stepIndex),8)
803 else
804 hu = hu_ptr_r8(lonIndex,latIndex,lev_T,stepIndex)
805 tt = tt_ptr_r8(lonIndex,latIndex,lev_T,stepIndex)
806 Pr = P_T_ptr_r8(lonIndex,latIndex,lev_T,stepIndex)
807 end if
808 cmp = gpscompressibility(Pr,tt,hu)
809
810 tv(lev_T) = phf_fotvt8(tt,hu) * cmp
811 end do
812
813 rMT = HeightSfc_ptr_r8(lonIndex,latIndex)
814
815 ! compute altitude on bottom momentum level
816 if (Vcode == 5002) then
817 height_M(nlev_M) = rMT
818 else if (Vcode == 5005 .or. Vcode == 5100) then
819 height_M(nlev_M) = rMT + heightSfcOffset_M_r4
820 end if
821
822 ! compute altitude on 2nd momentum level
823 if (nlev_M > 1) then
824 if ( statevector%dataKind == 4 ) then
825 P_M = real(P_M_ptr_r4(lonIndex,latIndex,nlev_M-1,stepIndex), 8)
826 P0 = real(P0_ptr_r4(lonIndex,latIndex,1,stepIndex), 8)
827 else
828 P_M = P_M_ptr_r8(lonIndex,latIndex,nlev_M-1,stepIndex)
829 P0 = P0_ptr_r8(lonIndex,latIndex,1,stepIndex)
830 end if
831
832 ratioP = log( P_M / P0 )
833
834 ! Gravity acceleration
835 h0 = rMT
836 Rgh = phf_gravityalt(sLat,h0)
837 dh = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(nlev_T-1) * ratioP
838 Rgh = phf_gravityalt(sLat, h0+0.5D0*dh)
839
840 delThick = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(nlev_T-1) * ratioP
841 height_M(nlev_M-1) = rMT + delThick
842 end if
843
844 ! compute altitude on rest of momentum levels
845 do lev_M = nlev_M-2, 1, -1
846 if ( statevector%dataKind == 4 ) then
847 P_M = real(P_M_ptr_r4(lonIndex,latIndex,lev_M,stepIndex),8)
848 P_M1 = real(P_M_ptr_r4(lonIndex,latIndex,lev_M+1,stepIndex),8)
849 else
850 P_M = P_M_ptr_r8(lonIndex,latIndex,lev_M,stepIndex)
851 P_M1 = P_M_ptr_r8(lonIndex,latIndex,lev_M+1,stepIndex)
852 end if
853
854 ratioP = log( P_M / P_M1 )
855
856 if (Vcode == 5002) then
857 lev_T = lev_M + 1
858 else if (Vcode == 5005 .or. Vcode == 5100) then
859 lev_T = lev_M
860 end if
861
862 ! Gravity acceleration
863 h0 = height_M(lev_M+1)
864 Rgh = phf_gravityalt(sLat,h0)
865 dh = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(lev_T) * ratioP
866 Rgh = phf_gravityalt(sLat, h0+0.5D0*dh)
867
868 delThick = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(lev_T) * ratioP
869 height_M(lev_M) = height_M(lev_M+1) + delThick
870 end do
871
872 ! compute Altitude on thermo levels
873 if_computeHeight_gsv_nl_vcodes : if (Vcode == 5002) then
874 height_T(nlev_T) = height_M(nlev_M)
875
876 do lev_T = 2, nlev_T-1
877 lev_M = lev_T ! momentum level just below thermo level being computed
878
879 if ( statevector%dataKind == 4 ) then
880 P_T = real(P_T_ptr_r4(&
881 lonIndex,latIndex,lev_T,stepIndex), 8)
882 P_M = real(P_M_ptr_r4(&
883 lonIndex,latIndex,lev_M,stepIndex), 8)
884 P_Mm1 = real(P_M_ptr_r4(&
885 lonIndex,latIndex,lev_M-1,stepIndex), 8)
886 else
887 P_T = P_T_ptr_r8(lonIndex,latIndex,lev_T ,stepIndex)
888 P_M = P_M_ptr_r8(lonIndex,latIndex,lev_M ,stepIndex)
889 P_Mm1 = P_M_ptr_r8(lonIndex,latIndex,lev_M-1,stepIndex)
890 end if
891
892 ScaleFactorBottom = log( P_T / P_Mm1 ) / log( P_M / P_Mm1 )
893 ScaleFactorTop = 1 - ScaleFactorBottom
894 height_T(lev_T) = ScaleFactorBottom * height_M(lev_M) &
895 + ScaleFactorTop * height_M(lev_M-1)
896 end do
897
898 ! compute altitude on top thermo level
899 if ( statevector%dataKind == 4 ) then
900 P_T = real(P_T_ptr_r4(lonIndex,latIndex,1,stepIndex),8)
901 P_M = real(P_M_ptr_r4(lonIndex,latIndex,1,stepIndex),8)
902 else
903 P_T = P_T_ptr_r8(lonIndex,latIndex,1,stepIndex)
904 P_M = P_M_ptr_r8(lonIndex,latIndex,1,stepIndex)
905 end if
906
907 ratioP = log( P_T / P_M )
908
909 ! Gravity acceleration
910 h0 = height_M(1)
911 Rgh = phf_gravityalt(sLat, h0)
912 dh = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(1) * ratioP
913 Rgh = phf_gravityalt(sLat, h0+0.5D0*dh)
914
915 delThick = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(1) * ratioP
916 height_T(1) = height_M(1) + delThick
917
918 else if (Vcode == 5005 .or. Vcode == 5100) then if_computeHeight_gsv_nl_vcodes
919 height_T(nlev_T) = rMT + heightSfcOffset_T_r4
920
921 do lev_T = 1, nlev_T-2
922 lev_M = lev_T + 1 ! momentum level just below thermo level being computed
923 if ( statevector%dataKind == 4 ) then
924 P_T = real(P_T_ptr_r4(&
925 lonIndex,latIndex,lev_T,stepIndex), 8)
926 P_M = real(P_M_ptr_r4(&
927 lonIndex,latIndex,lev_M,stepIndex), 8)
928 P_Mm1 = real(P_M_ptr_r4(&
929 lonIndex,latIndex,lev_M-1,stepIndex), 8)
930 else
931 P_T = P_T_ptr_r8(lonIndex,latIndex,lev_T ,stepIndex)
932 P_M = P_M_ptr_r8(lonIndex,latIndex,lev_M ,stepIndex)
933 P_Mm1 = P_M_ptr_r8(lonIndex,latIndex,lev_M-1,stepIndex)
934 end if
935
936 ScaleFactorBottom = log( P_T / P_Mm1 ) / log( P_M / P_Mm1 )
937 ScaleFactorTop = 1 - ScaleFactorBottom
938 height_T(lev_T) = ScaleFactorBottom * height_M(lev_M) &
939 + ScaleFactorTop * height_M(lev_M-1)
940 end do
941
942 ! compute altitude on next to bottom thermo level
943 if (nlev_T > 1) then
944 if ( statevector%dataKind == 4 ) then
945 P_T = real(P_T_ptr_r4(&
946 lonIndex,latIndex,nlev_T-1,stepIndex), 8)
947 P0 = real(P0_ptr_r4(&
948 lonIndex,latIndex,1,stepIndex), 8)
949 else
950 P_T = P_T_ptr_r8(lonIndex,latIndex,nlev_T-1,stepIndex)
951 P0 = P0_ptr_r8(lonIndex,latIndex,1,stepIndex)
952 end if
953
954 ratioP = log( P_T / P0 )
955
956 h0 = rMT
957 Rgh = phf_gravityalt(sLat,h0)
958 dh = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(nlev_T-1) * ratioP
959 Rgh = phf_gravityalt(sLat, h0+0.5D0*dh)
960
961 delThick = (-MPC_RGAS_DRY_AIR_R8 / Rgh) * tv(nlev_T-1) * &
962 ratioP
963 height_T(nlev_T-1) = rMT + delThick
964 end if
965 end if if_computeHeight_gsv_nl_vcodes
966
967 ! fill the height array
968 if ( statevector%dataKind == 4 ) then
969 do lev_T = 1, nlev_T
970 height_T_ptr_r4(lonIndex,latIndex,lev_T,stepIndex) &
971 = real(height_T(lev_T),4)
972 end do
973 do lev_M = 1, nlev_M
974 height_M_ptr_r4(lonIndex,latIndex,lev_M,stepIndex) &
975 = real(height_M(lev_M),4)
976 end do
977 else
978 height_T_ptr_r8(lonIndex,latIndex,1:nlev_T,stepIndex) &
979 = height_T(1:nlev_T)
980 height_M_ptr_r8(lonIndex,latIndex,1:nlev_M,stepIndex) &
981 = height_M(1:nlev_M)
982 end if
983
984 ! remove the height offset for the diagnostic levels for backward compatibility only
985 if ( .not. statevector%addHeightSfcOffset ) then
986 if ( statevector%dataKind == 4 ) then
987 height_T_ptr_r4(lonIndex,latIndex,nlev_T,stepIndex) = &
988 real(rMT,4)
989 height_M_ptr_r4(lonIndex,latIndex,nlev_M,stepIndex) = &
990 real(rMT,4)
991 else
992 height_T_ptr_r8(lonIndex,latIndex,nlev_T,stepIndex) = rMT
993 height_M_ptr_r8(lonIndex,latIndex,nlev_M,stepIndex) = rMT
994 end if
995 end if
996
997 end do
998 end do
999 !$OMP end parallel do
1000 end do do_computeHeight_gsv_nl
1001
1002 deallocate(height_M)
1003 deallocate(height_T)
1004 deallocate(tv)
1005
1006 call msg('calcHeight_gsv_nl_vcode5xxx (czp)', 'statevector%addHeightSfcOffset='&
1007 //str(statevector%addHeightSfcOffset), verb_opt=2)
1008 call msg('calcHeight_gsv_nl_vcode5xxx (czp)', 'END', verb_opt=4)
1009 end subroutine calcHeight_gsv_nl_vcode5xxx
1010
1011 !---------------------------------------------------------
1012 ! calcHeight_gsv_tl
1013 !---------------------------------------------------------
1014 subroutine calcHeight_gsv_tl(statevector,statevectorRef)
1015 !
1016 !:Purpose: Tangent height computation on statevector.
1017 !
1018 implicit none
1019
1020 ! Arguments:
1021 type(struct_gsv), intent(inout) :: statevector
1022 type(struct_gsv), intent(in) :: statevectorRef
1023
1024 ! Locals:
1025 integer :: Vcode
1026
1027 call utl_tmg_start(173,'low-level--czp_calcHeight_tl')
1028 call msg('calcHeight_gsv_tl (czp)', 'START', verb_opt=2)
1029
1030 Vcode = gsv_getVco(statevectorRef)%vcode
1031 if (Vcode == 5005 .or. Vcode == 5002) then
1032 if ( .not. gsv_varExist(statevector,'P_*') ) then
1033 call utl_abort('calcHeight_gsv_tl (czp): for vcode 5xxx, variables P_T and P_M must be allocated in gridstatevector')
1034 end if
1035 if ( .not. gsv_varExist(statevector,'Z_*') ) then
1036 call utl_abort('calcHeight_gsv_tl (czp): for vcode 5xxx, variables Z_T and Z_M must be allocated in gridstatevector')
1037 end if
1038 if ( .not. gsv_varExist(statevector,'TT') ) then
1039 call utl_abort('calcHeight_gsv_tl (czp): for vcode 5xxx, variable TT must be allocated in gridstatevector')
1040 end if
1041 if ( .not. gsv_varExist(statevector,'HU') ) then
1042 call utl_abort('calcHeight_gsv_tl (czp): for vcode 5xxx, variable HU must be allocated in gridstatevector')
1043 end if
1044 if ( .not. gsv_varExist(statevector,'P0') ) then
1045 call utl_abort('calcHeight_gsv_tl (czp): for vcode 5xxx, variable P0 must be allocated in gridstatevector')
1046 end if
1047 call calcHeight_gsv_tl_vcode5xxx
1048 else if (Vcode == 21001) then
1049 ! Development notes (@mad001)
1050 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
1051 call calcHeight_gsv_tl_vcode2100x
1052 else
1053 call utl_abort('calcHeight_gsv_tl (czp): not implemented')
1054 end if
1055
1056 call msg('calcHeight_gsv_tl (czp)', 'END', verb_opt=2)
1057 call utl_tmg_stop(173)
1058
1059 contains
1060 !---------------------------------------------------------
1061 ! calcHeight_gsv_tl_vcode2100x
1062 !---------------------------------------------------------
1063 subroutine calcHeight_gsv_tl_vcode2100x
1064 implicit none
1065
1066 call utl_abort('calcHeight_gsv_tl (czp): vcode 21001 not implemented yet')
1067
1068 end subroutine calcHeight_gsv_tl_vcode2100x
1069
1070 !---------------------------------------------------------
1071 ! calcHeight_gsv_tl_vcode5xxx
1072 !---------------------------------------------------------
1073 subroutine calcHeight_gsv_tl_vcode5xxx
1074 implicit none
1075
1076 ! Locals:
1077 integer :: lev_M,lev_T,nlev_M,nlev_T,Vcode_anl
1078 integer :: numStep,stepIndex, latIndex,lonIndex
1079 real(8) :: ScaleFactorBottom, ScaleFactorTop
1080 real(8), allocatable :: delThick(:,:,:,:)
1081 real(8), pointer :: height_T_ptr(:,:,:,:),height_M_ptr(:,:,:,:)
1082 real(8), pointer :: P_T(:,:,:,:), P_M(:,:,:,:)
1083 real(pre_incrReal), pointer :: delHeight_M_ptr_r48(:,:,:,:)
1084 real(pre_incrReal), pointer :: delHeight_T_ptr_r48(:,:,:,:)
1085 real(pre_incrReal), pointer :: delTT_r48(:,:,:,:), delHU_r48(:,:,:,:)
1086 real(pre_incrReal), pointer :: delP0_r48(:,:,:,:)
1087 real(pre_incrReal), pointer :: delP_T_r48(:,:,:,:), delP_M_r48(:,:,:,:)
1088
1089 call msg('calcHeight_gsv_tl_vcode5xxx (czp)', 'START', verb_opt=4)
1090 Vcode_anl = gsv_getVco(statevectorRef)%vcode
1091
1092 nlev_T = gsv_getNumLev(statevectorRef,'TH')
1093 nlev_M = gsv_getNumLev(statevectorRef,'MM')
1094 numStep = statevectorRef%numstep
1095
1096 allocate(delThick(statevectorRef%myLonBeg:statevectorRef%myLonEnd, &
1097 statevectorRef%myLatBeg:statevectorRef%myLatEnd, &
1098 nlev_T,numStep))
1099
1100 ! generate the height coefficients on the grid
1101 call calcHeightCoeff_gsv(statevectorRef)
1102
1103 ! loop over all lat/lon/step
1104
1105 call gsv_getField(statevectorRef,height_M_ptr,'Z_M')
1106 call gsv_getField(statevectorRef,height_T_ptr,'Z_T')
1107 call gsv_getField(statevectorRef,P_T,'P_T')
1108 call gsv_getField(statevectorRef,P_M,'P_M')
1109 call gsv_getField(statevector,delHeight_M_ptr_r48,'Z_M')
1110 call gsv_getField(statevector,delHeight_T_ptr_r48,'Z_T')
1111 call gsv_getField(statevector,delTT_r48,'TT')
1112 call gsv_getField(statevector,delHU_r48,'HU')
1113 call gsv_getField(statevector,delP0_r48,'P0')
1114 call gsv_getField(statevector,delP_T_r48,'P_T')
1115 call gsv_getField(statevector,delP_M_r48,'P_M')
1116 ! ensure increment at sfc is zero (fixed height level)
1117 delHeight_M_ptr_r48(:,:,nlev_M,:) = 0.0d0
1118 delHeight_T_ptr_r48(:,:,nlev_T,:) = 0.0d0
1119
1120 if_computeHeight_gsv_tl_vcodes : if(Vcode_anl == 5002) then
1121
1122 ! compute increment to thickness for each layer between the two momentum levels
1123 do stepIndex = 1, numStep
1124 do lev_T = 2, (nlev_T-1)
1125 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1126 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1127 delThick(lonIndex,latIndex,lev_T,stepIndex) = &
1128 coeff_M_TT_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1129 delTT_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1130 coeff_M_HU_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1131 delHU_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1132 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) *&
1133 ( delP_M_r48(lonIndex,latIndex,lev_T ,stepIndex) / &
1134 P_M(lonIndex,latIndex,lev_T ,stepIndex) - &
1135 delP_M_r48(lonIndex,latIndex,lev_T-1,stepIndex) / &
1136 P_M(lonIndex,latIndex,lev_T-1,stepIndex) ) + &
1137 coeff_M_P0_dP_delPT_gsv(&
1138 lonIndex,latIndex,lev_T,stepIndex) * &
1139 delP_T_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1140 coeff_M_P0_dP_delP0_gsv(&
1141 lonIndex,latIndex,lev_T,stepIndex) * &
1142 delP0_r48(lonIndex,latIndex,1,stepIndex)
1143 end do
1144 end do
1145 end do
1146 end do
1147
1148 ! compute height increment on momentum levels above the surface
1149 do stepIndex = 1, numStep
1150 do lev_M = (nlev_M-1), 1, -1
1151 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1152 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1153 lev_T = lev_M + 1 ! thermo level just below momentum level being computed
1154 delHeight_M_ptr_r48(lonIndex,latIndex,lev_M,stepIndex) = &
1155 delHeight_M_ptr_r48(&
1156 lonIndex,latIndex,lev_M+1,stepIndex) + &
1157 delThick(lonIndex,latIndex,lev_T,stepIndex)
1158 end do
1159 end do
1160 end do
1161 end do
1162
1163 ! compute height increment on thermo levels using weighted average of height increment of momentum levels
1164 do stepIndex = 1, numStep
1165 do lev_T = 1, (nlev_T-1)
1166 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1167 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1168 if ( lev_T == 1) then
1169 ! compute height increment for top thermo level (from top momentum level)
1170 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex) = &
1171 delHeight_M_ptr_r48(lonIndex,latIndex,1,stepIndex) + &
1172 coeff_T_TT_gsv(lonIndex,latIndex,stepIndex) * &
1173 delTT_r48(lonIndex,latIndex,1,stepIndex) + &
1174 coeff_T_HU_gsv(lonIndex,latIndex,stepIndex) * &
1175 delHU_r48(lonIndex,latIndex,1,stepIndex) + &
1176 coeff_T_P0_delP1_gsv(lonIndex,latIndex,stepIndex) * &
1177 ( delP_M_r48(lonIndex,latIndex,1,stepIndex) / &
1178 P_M(lonIndex,latIndex,1,stepIndex) - &
1179 delP_T_r48(lonIndex,latIndex,1,stepIndex) / &
1180 P_T(lonIndex,latIndex,1,stepIndex) ) + &
1181 coeff_T_P0_dP_delPT_gsv(lonIndex,latIndex,stepIndex) * &
1182 delP_T_r48(lonIndex,latIndex,1,stepIndex) + &
1183 coeff_T_P0_dP_delP0_gsv(lonIndex,latIndex,stepIndex) * &
1184 delP0_r48(lonIndex,latIndex,1,stepIndex)
1185 else
1186 lev_M = lev_T ! momentum level just below thermo level being computed
1187 ScaleFactorBottom = &
1188 (height_T_ptr(lonIndex,latIndex,lev_T,stepIndex) - &
1189 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex)) / &
1190 (height_M_ptr(lonIndex,latIndex,lev_M,stepIndex) - &
1191 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex))
1192 ScaleFactorTop = 1 - ScaleFactorBottom
1193 delHeight_T_ptr_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1194 ScaleFactorBottom * &
1195 delHeight_M_ptr_r48(&
1196 lonIndex,latIndex,lev_M ,stepIndex) + &
1197 ScaleFactorTop * &
1198 delHeight_M_ptr_r48(&
1199 lonIndex,latIndex,lev_M-1,stepIndex)
1200 end if
1201 end do
1202 end do
1203 end do
1204 end do
1205
1206 else if(Vcode_anl == 5005) then if_computeHeight_gsv_tl_vcodes
1207
1208 ! compute increment to thickness for each layer between the two momentum levels
1209 do stepIndex = 1, numStep
1210 do lev_T = 1, (nlev_T-1)
1211 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1212 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1213 delThick(lonIndex,latIndex,lev_T,stepIndex) = &
1214 coeff_M_TT_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1215 delTT_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1216 coeff_M_HU_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1217 delHU_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1218 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) *&
1219 ( delP_M_r48(lonIndex,latIndex,lev_T+1,stepIndex) / &
1220 P_M(lonIndex,latIndex,lev_T+1,stepIndex) - &
1221 delP_M_r48(lonIndex,latIndex,lev_T ,stepIndex) / &
1222 P_M(lonIndex,latIndex,lev_T ,stepIndex) ) + &
1223 coeff_M_P0_dP_delPT_gsv(&
1224 lonIndex,latIndex,lev_T,stepIndex) * &
1225 delP_T_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1226 coeff_M_P0_dP_delP0_gsv(&
1227 lonIndex,latIndex,lev_T,stepIndex) * &
1228 delP0_r48(lonIndex,latIndex,1,stepIndex)
1229 end do
1230 end do
1231 end do
1232 end do
1233
1234 ! compute height increment on momentum levels above the surface
1235 do stepIndex = 1, numStep
1236 do lev_M = (nlev_M-1), 1, -1
1237 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1238 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1239 lev_T = lev_M ! thermo level just below momentum level being computed
1240 delHeight_M_ptr_r48(lonIndex,latIndex,lev_M,stepIndex) = &
1241 delHeight_M_ptr_r48(lonIndex,latIndex,lev_M+1,stepIndex) + &
1242 delThick(lonIndex,latIndex,lev_T,stepIndex)
1243 end do
1244 end do
1245 end do
1246 end do
1247
1248 ! compute height increment on thermo levels using weighted average of height increment of momentum levels
1249 do stepIndex = 1, numStep
1250 do lev_T = 1, (nlev_T-1)
1251 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1252 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1253 lev_M = lev_T + 1 ! momentum level just below thermo level being computed
1254 ScaleFactorBottom = &
1255 (height_T_ptr(lonIndex,latIndex,lev_T,stepIndex) - &
1256 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex)) / &
1257 (height_M_ptr(lonIndex,latIndex,lev_M,stepIndex) - &
1258 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex))
1259 ScaleFactorTop = 1 - ScaleFactorBottom
1260 delHeight_T_ptr_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1261 ScaleFactorBottom * &
1262 delHeight_M_ptr_r48(lonIndex,latIndex,lev_M ,stepIndex) + &
1263 ScaleFactorTop * &
1264 delHeight_M_ptr_r48(lonIndex,latIndex,lev_M-1,stepIndex)
1265 end do
1266 end do
1267 end do
1268 end do
1269
1270 else
1271
1272 call utl_abort('calcHeight_gsv_tl_vcode5xxx (czp): not implemented')
1273
1274 end if if_computeHeight_gsv_tl_vcodes
1275
1276 deallocate(delThick)
1277
1278 call msg('calcHeight_gsv_tl_vcode5xxx (czp)', 'END', verb_opt=4)
1279 end subroutine calcHeight_gsv_tl_vcode5xxx
1280
1281 end subroutine calcHeight_gsv_tl
1282
1283 !---------------------------------------------------------
1284 ! calcHeight_gsv_ad
1285 !---------------------------------------------------------
1286 subroutine calcHeight_gsv_ad(statevector,statevectorRef)
1287 !
1288 !:Purpose: Adjoint of height computation on statevector.
1289 !
1290 implicit none
1291
1292 ! Arguments:
1293 type(struct_gsv), intent(inout) :: statevector
1294 type(struct_gsv), intent(in) :: statevectorRef
1295
1296 ! Locals:
1297 integer :: Vcode
1298
1299 call utl_tmg_start(174,'low-level--czp_calcHeight_ad')
1300 call msg('calcHeight_gsv_ad', 'START', verb_opt=2)
1301
1302 Vcode = gsv_getVco(statevectorRef)%vcode
1303 if (Vcode == 5005 .or. Vcode == 5002) then
1304 if ( .not. gsv_varExist(statevector,'P_*') ) then
1305 call utl_abort('calcHeight_gsv_ad (czp): for vcode 5xxx, variables P_M and P_T must be allocated in gridstatevector')
1306 end if
1307 if ( .not. gsv_varExist(statevector,'Z_*') ) then
1308 call utl_abort('calcHeight_gsv_ad (czp): for vcode 5xxx, variables Z_M and Z_T must be allocated in gridstatevector')
1309 end if
1310 if ( .not. gsv_varExist(statevector,'TT') ) then
1311 call utl_abort('calcHeight_gsv_ad (czp): for vcode 5xxx, variable TT must be allocated in gridstatevector')
1312 end if
1313 if ( .not. gsv_varExist(statevector,'HU') ) then
1314 call utl_abort('calcHeight_gsv_ad (czp): for vcode 5xxx, variable HU must be allocated in gridstatevector')
1315 end if
1316 if ( .not. gsv_varExist(statevector,'P0') ) then
1317 call utl_abort('calcHeight_gsv_ad (czp): for vcode 5xxx, variable P0 must be allocated in gridstatevector')
1318 end if
1319 call calcHeight_gsv_ad_vcode5xxx
1320 else if (Vcode == 21001) then
1321 ! Development notes (@mad001)
1322 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
1323 call calcHeight_gsv_ad_vcode2100x
1324 else
1325 call utl_abort('calcHeight_gsv_ad (czp): not implemented')
1326 end if
1327
1328 call msg('calcHeight_gsv_ad (czp)', 'END', verb_opt=2)
1329 call utl_tmg_stop(174)
1330
1331 contains
1332 !---------------------------------------------------------
1333 ! calcHeight_gsv_ad_vcode2100x
1334 !---------------------------------------------------------
1335 subroutine calcHeight_gsv_ad_vcode2100x
1336 implicit none
1337
1338 call utl_abort('calcHeight_gsv_ad (czp): vcode 21001 not implemented yet')
1339
1340 end subroutine calcHeight_gsv_ad_vcode2100x
1341
1342 !---------------------------------------------------------
1343 ! calcHeight_gsv_ad
1344 !---------------------------------------------------------
1345 subroutine calcHeight_gsv_ad_vcode5xxx
1346 implicit none
1347
1348 ! Locals:
1349 integer :: lev_M,lev_T,nlev_M,nlev_T
1350 integer :: numStep,stepIndex,latIndex,lonIndex
1351 real(8) :: ScaleFactorBottom, ScaleFactorTop
1352 real(8), allocatable :: delThick(:,:,:,:)
1353 real(8), pointer :: height_M_ptr(:,:,:,:),height_T_ptr(:,:,:,:)
1354 real(8), allocatable :: delHeight_M(:,:,:,:)
1355 real(8), pointer :: P_M(:,:,:,:),P_T(:,:,:,:)
1356 real(pre_incrReal), pointer :: delHeight_M_ptr_r48(:,:,:,:)
1357 real(pre_incrReal), pointer :: delHeight_T_ptr_r48(:,:,:,:)
1358 real(pre_incrReal), pointer :: delTT_r48(:,:,:,:),delHU_r48(:,:,:,:)
1359 real(pre_incrReal), pointer :: delP0_r48(:,:,:,:)
1360 real(pre_incrReal), pointer :: delP_M_r48(:,:,:,:),delP_T_r48(:,:,:,:)
1361
1362 call msg('calcHeight_gsv_ad_vcode5xxx (czp)', 'START', verb_opt=4)
1363
1364 nlev_T = gsv_getNumLev(statevectorRef,'TH')
1365 nlev_M = gsv_getNumLev(statevectorRef,'MM')
1366 numStep = statevectorRef%numstep
1367
1368 allocate(delHeight_M(statevectorRef%myLonBeg:statevectorRef%myLonEnd, &
1369 statevectorRef%myLatBeg:statevectorRef%myLatEnd, &
1370 nlev_M,numStep))
1371 allocate(delThick(statevectorRef%myLonBeg:statevectorRef%myLonEnd, &
1372 statevectorRef%myLatBeg:statevectorRef%myLatEnd, &
1373 nlev_T,numStep))
1374
1375 ! generate the height coefficients on the grid
1376 call calcHeightCoeff_gsv(statevectorRef)
1377
1378 ! loop over all lat/lon/step
1379 call gsv_getField(statevectorRef,height_M_ptr,'Z_M')
1380 call gsv_getField(statevectorRef,height_T_ptr,'Z_T')
1381 call gsv_getField(statevectorRef,P_T,'P_T')
1382 call gsv_getField(statevectorRef,P_M,'P_M')
1383
1384 call gsv_getField(statevector,delHeight_M_ptr_r48,'Z_M')
1385 call gsv_getField(statevector,delHeight_T_ptr_r48,'Z_T')
1386 call gsv_getField(statevector,delTT_r48,'TT')
1387 call gsv_getField(statevector,delHU_r48,'HU')
1388 call gsv_getField(statevector,delP0_r48,'P0')
1389 call gsv_getField(statevector,delP_T_r48,'P_T')
1390 call gsv_getField(statevector,delP_M_r48,'P_M')
1391 delHeight_M(:,:,:,:) = delHeight_M_ptr_r48(:,:,:,:)
1392
1393 if_computeHeight_gsv_ad_vcodes : if(Vcode == 5002) then
1394
1395 ! adjoint of compute height increment on thermo levels by simple averaging
1396 do stepIndex = 1, numStep
1397 do lev_T = 1, (nlev_T-1)
1398 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1399 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1400 lev_M = lev_T ! momentum level just below thermo level being computed
1401
1402 ! adjoint of compute height increment on top thermo level
1403 ! (from top momentum level)
1404 if (lev_T == 1) then
1405 delHeight_M(lonIndex,latIndex,1,stepIndex) = &
1406 delHeight_M(lonIndex,latIndex,1,stepIndex) + &
1407 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1408
1409 delTT_r48(lonIndex,latIndex,1,stepIndex) = &
1410 delTT_r48(lonIndex,latIndex,1,stepIndex) + &
1411 coeff_T_TT_gsv(lonIndex,latIndex,stepIndex) * &
1412 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1413
1414 delHU_r48(lonIndex,latIndex,1,stepIndex) = &
1415 delHU_r48(lonIndex,latIndex,1,stepIndex) + &
1416 coeff_T_HU_gsv (lonIndex,latIndex,stepIndex) * &
1417 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1418
1419 delP_M_r48(lonIndex,latIndex,1,stepIndex) = &
1420 delP_M_r48(lonIndex,latIndex,1,stepIndex) + &
1421 coeff_T_P0_delP1_gsv(lonIndex,latIndex,stepIndex) / &
1422 P_M(lonIndex,latIndex,1,stepIndex) * &
1423 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1424
1425 delP_T_r48(lonIndex,latIndex,1,stepIndex) = &
1426 delP_T_r48(lonIndex,latIndex,1,stepIndex) - &
1427 coeff_T_P0_delP1_gsv(lonIndex,latIndex,stepIndex) / &
1428 P_T(lonIndex,latIndex,1,stepIndex) * &
1429 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1430
1431 delP_T_r48(lonIndex,latIndex,1,stepIndex) = &
1432 delP_T_r48(lonIndex,latIndex,1,stepIndex) + &
1433 coeff_T_P0_dP_delPT_gsv(lonIndex,latIndex,stepIndex) * &
1434 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1435
1436 delP0_r48(lonIndex,latIndex,1,stepIndex) = &
1437 delP0_r48(lonIndex,latIndex,1,stepIndex) + &
1438 coeff_T_P0_dp_delP0_gsv(lonIndex,latIndex,stepIndex) * &
1439 delHeight_T_ptr_r48(lonIndex,latIndex,1,stepIndex)
1440 else
1441 ScaleFactorBottom = &
1442 (height_T_ptr(lonIndex,latIndex,lev_T,stepIndex) - &
1443 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex)) / &
1444 (height_M_ptr(lonIndex,latIndex,lev_M,stepIndex) - &
1445 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex))
1446 ScaleFactorTop = 1 - ScaleFactorBottom
1447
1448 delHeight_M(lonIndex,latIndex,lev_M-1,stepIndex) = &
1449 delHeight_M(lonIndex,latIndex,lev_M-1,stepIndex) + &
1450 ScaleFactorTop * &
1451 delHeight_T_ptr_r48(lonIndex,latIndex,lev_T,stepIndex)
1452
1453 delHeight_M(lonIndex,latIndex,lev_M,stepIndex) = &
1454 delHeight_M(lonIndex,latIndex,lev_M ,stepIndex) + &
1455 ScaleFactorBottom * &
1456 delHeight_T_ptr_r48(lonIndex,latIndex,lev_T,stepIndex)
1457 end if
1458 end do
1459 end do
1460 end do
1461 end do
1462
1463 ! adjoint of compute height increment on momentum levels above the surface
1464 delThick(:,:,1,:) = 0.0d0
1465 do stepIndex = 1, numStep
1466 do lev_M = 1, (nlev_M-1)
1467 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1468 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1469 lev_T = lev_M + 1 ! thermo level just below momentum level being computed
1470 delThick(lonIndex,latIndex,lev_T,stepIndex) = &
1471 delThick(lonIndex,latIndex,lev_T-1,stepIndex) + &
1472 delHeight_M(lonIndex,latIndex,lev_M ,stepIndex)
1473 end do
1474 end do
1475 end do
1476 end do
1477
1478 ! adjoint of compute increment to thickness for each layer between the two momentum levels
1479 do stepIndex = 1, numStep
1480 do lev_T = 2, nlev_T-1
1481 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1482 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1483 delTT_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1484 delTT_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1485 coeff_M_TT_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1486 delThick(lonIndex,latIndex,lev_T,stepIndex)
1487
1488 delHU_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1489 delHU_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1490 coeff_M_HU_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1491 delThick(lonIndex,latIndex,lev_T,stepIndex)
1492
1493 delP_M_r48(lonIndex,latIndex,lev_T,stepIndex)= &
1494 delP_M_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1495 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) /&
1496 P_M(lonIndex,latIndex,lev_T,stepIndex) * &
1497 delThick(lonIndex,latIndex,lev_T,stepIndex)
1498
1499 delP_M_r48(lonIndex,latIndex,lev_T-1,stepIndex) = &
1500 delP_M_r48(lonIndex,latIndex,lev_T-1,stepIndex) - &
1501 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) /&
1502 P_M(lonIndex,latIndex,lev_T-1,stepIndex) * &
1503 delThick(lonIndex,latIndex,lev_T,stepIndex)
1504
1505 delP_T_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1506 delP_T_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1507 coeff_M_P0_dP_delPT_gsv(&
1508 lonIndex,latIndex,lev_T,stepIndex) * &
1509 delThick(lonIndex,latIndex,lev_T,stepIndex)
1510
1511 delP0_r48(lonIndex,latIndex,1,stepIndex) = &
1512 delP0_r48(lonIndex,latIndex,1,stepIndex) + &
1513 coeff_M_P0_dP_delP0_gsv(&
1514 lonIndex,latIndex,lev_T,stepIndex) * &
1515 delThick(lonIndex,latIndex,lev_T,stepIndex)
1516 end do
1517 end do
1518 end do
1519 end do
1520
1521 else if(Vcode == 5005) then if_computeHeight_gsv_ad_vcodes
1522
1523 ! adjoint of compute height increment on thermo levels by simple averaging
1524 do stepIndex = 1, numStep
1525 do lev_T = 1, (nlev_T-1)
1526 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1527 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1528 lev_M = lev_T+1 ! momentum level just below thermo level being computed
1529 ScaleFactorBottom = &
1530 (height_T_ptr(lonIndex,latIndex,lev_T,stepIndex) - &
1531 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex)) / &
1532 (height_M_ptr(lonIndex,latIndex,lev_M,stepIndex) -&
1533 height_M_ptr(lonIndex,latIndex,lev_M-1,stepIndex))
1534 ScaleFactorTop = 1 - ScaleFactorBottom
1535 delHeight_M(lonIndex,latIndex,lev_M-1,stepIndex) = &
1536 delHeight_M(lonIndex,latIndex,lev_M-1,stepIndex) + &
1537 ScaleFactorTop * &
1538 delHeight_T_ptr_r48(lonIndex,latIndex,lev_T,stepIndex)
1539 delHeight_M(lonIndex,latIndex,lev_M,stepIndex) = &
1540 delHeight_M(lonIndex,latIndex,lev_M ,stepIndex) + &
1541 ScaleFactorBottom * &
1542 delHeight_T_ptr_r48(lonIndex,latIndex,lev_T,stepIndex)
1543 end do
1544 end do
1545 end do
1546 end do
1547
1548 ! adjoint of compute height increment on momentum levels
1549 do stepIndex = 1, numStep
1550 do lev_M = 1, (nlev_M-1)
1551 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1552 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1553 lev_T = lev_M ! thermo level just below momentum level being computed
1554 if (lev_T == 1) then
1555 delThick(lonIndex,latIndex,lev_T,stepIndex) = &
1556 delHeight_M (lonIndex,latIndex,lev_M,stepIndex)
1557 else
1558 delThick(lonIndex,latIndex,lev_T,stepIndex) = &
1559 delThick(lonIndex,latIndex,lev_T-1,stepIndex) + &
1560 delHeight_M (lonIndex,latIndex,lev_M,stepIndex)
1561 end if
1562 end do
1563 end do
1564 end do
1565 end do
1566
1567 do stepIndex = 1, numStep
1568 do lev_T = 1, nlev_T-1
1569 do latIndex = statevectorRef%myLatBeg, statevectorRef%myLatEnd
1570 do lonIndex = statevectorRef%myLonBeg, statevectorRef%myLonEnd
1571 delTT_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1572 delTT_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1573 coeff_M_TT_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1574 delThick(lonIndex,latIndex,lev_T,stepIndex)
1575
1576 delHU_r48(lonIndex,latIndex,lev_T,stepIndex) = &
1577 delHU_r48(lonIndex,latIndex,lev_T,stepIndex) + &
1578 coeff_M_HU_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1579 delThick(lonIndex,latIndex,lev_T,stepIndex)
1580
1581 delP_M_r48(lonIndex,latIndex,lev_T+1,stepIndex) = &
1582 delP_M_r48(lonIndex,latIndex,lev_T+1,stepIndex) + &
1583 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) /&
1584 P_M(lonIndex,latIndex,lev_T+1,stepIndex) * &
1585 delThick(lonIndex,latIndex,lev_T,stepIndex)
1586
1587 delP_M_r48(lonIndex,latIndex,lev_T ,stepIndex) = &
1588 delP_M_r48(lonIndex,latIndex,lev_T ,stepIndex) - &
1589 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) /&
1590 P_M(lonIndex,latIndex,lev_T ,stepIndex) * &
1591 delThick(lonIndex,latIndex,lev_T,stepIndex)
1592
1593 delP_T_r48(lonIndex,latIndex,lev_T ,stepIndex) = &
1594 delP_T_r48(lonIndex,latIndex,lev_T ,stepIndex) + &
1595 coeff_M_P0_dP_delPT_gsv( lonIndex,latIndex,lev_T,stepIndex) * &
1596 delThick(lonIndex,latIndex,lev_T,stepIndex)
1597
1598 delP0_r48(lonIndex,latIndex,1,stepIndex) = &
1599 delP0_r48(lonIndex,latIndex,1,stepIndex) + &
1600 coeff_M_P0_dP_delP0_gsv(lonIndex,latIndex,lev_T,stepIndex) * &
1601 delThick(lonIndex,latIndex,lev_T,stepIndex)
1602 end do
1603 end do
1604 end do
1605 end do
1606
1607 else
1608
1609 call utl_abort('calcHeight_gsv_ad_vcode5xxx (czp): not implemented')
1610
1611 end if if_computeHeight_gsv_ad_vcodes
1612
1613 deallocate(delThick)
1614 deallocate(delHeight_M)
1615
1616 call msg('calcHeight_gsv_ad_vcode5xxx (czp)', 'END', verb_opt=4)
1617 end subroutine calcHeight_gsv_ad_vcode5xxx
1618
1619 end subroutine calcHeight_gsv_ad
1620
1621 !---------------------------------------------------------
1622 ! calcPressure_gsv_nl
1623 !---------------------------------------------------------
1624 subroutine calcPressure_gsv_nl(statevector, Ps_in_hPa_opt)
1625 !
1626 !:Purpose: Pressure computation, values stored in statevector.
1627 !
1628 implicit none
1629
1630 ! Arguments:
1631 type(struct_gsv), intent(inout) :: statevector
1632 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, conversion from hPa to mbar done for surface pressure
1633
1634 ! Locals:
1635 integer :: Vcode
1636 real(4), pointer :: ptr_ZT_r4(:,:,:,:), ptr_ZM_r4(:,:,:,:)
1637 real(8), pointer :: ptr_ZT_r8(:,:,:,:), ptr_ZM_r8(:,:,:,:)
1638 real(4), pointer :: ptr_PT_r4(:,:,:,:), ptr_PM_r4(:,:,:,:)
1639 real(8), pointer :: ptr_PT_r8(:,:,:,:), ptr_PM_r8(:,:,:,:)
1640
1641 call utl_tmg_start(177,'low-level--czp_calcPressure_nl')
1642 call msg('calcPressure_gsv_nl (czp)', 'START', verb_opt=2)
1643
1644 Vcode = gsv_getVco(statevector)%vcode
1645 if (Vcode == 5005 .or. Vcode == 5002 .or. Vcode == 5100) then
1646 if ( gsv_getDataKind(statevector) == 4 ) then
1647 call gsv_getField(statevector, ptr_PT_r4, 'P_T')
1648 call gsv_getField(statevector, ptr_PM_r4, 'P_M')
1649 call calcPressure_gsv_nl_vcode5xxx_r4(statevector, &
1650 ptr_PT_r4, ptr_PM_r4, &
1651 Ps_in_hPa_opt=Ps_in_hPa_opt)
1652 else
1653 call gsv_getField(statevector, ptr_PT_r8, 'P_T')
1654 call gsv_getField(statevector, ptr_PM_r8, 'P_M')
1655 call calcPressure_gsv_nl_vcode5xxx_r8(statevector, &
1656 ptr_PT_r8, ptr_PM_r8, &
1657 Ps_in_hPa_opt=Ps_in_hPa_opt)
1658 end if
1659 else if (Vcode == 21001) then
1660 ! Development notes (@mad001)
1661 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
1662 if ( gsv_getDataKind(statevector) == 4 ) then
1663 call gsv_getField(statevector, ptr_ZT_r4, 'Z_T')
1664 call gsv_getField(statevector, ptr_ZM_r4, 'Z_M')
1665 call gsv_getField(statevector, ptr_PT_r4, 'P_T')
1666 call gsv_getField(statevector, ptr_PM_r4, 'P_M')
1667 call calcPressure_gsv_nl_vcode2100x(statevector, &
1668 ZTin_r4_opt=ptr_ZT_r4, &
1669 ZMin_r4_opt=ptr_ZM_r4, &
1670 PTout_r4_opt=ptr_PT_r4, &
1671 PMout_r4_opt=ptr_PM_r4)
1672 else
1673 call gsv_getField(statevector, ptr_ZT_r8, 'Z_T')
1674 call gsv_getField(statevector, ptr_ZM_r8, 'Z_M')
1675 call gsv_getField(statevector, ptr_PT_r8, 'P_T')
1676 call gsv_getField(statevector, ptr_PM_r8, 'P_M')
1677 call calcPressure_gsv_nl_vcode2100x(statevector, &
1678 ZTin_r8_opt=ptr_ZT_r8, &
1679 ZMin_r8_opt=ptr_ZM_r8, &
1680 PTout_r8_opt=ptr_PT_r8, &
1681 PMout_r8_opt=ptr_PM_r8)
1682 end if
1683 end if
1684
1685 if ( gsv_getDataKind(statevector) == 4 ) then
1686 call msg('calcPressure_gsv_nl (czp)', &
1687 new_line('')//'P_M = '&
1688 //str(ptr_PM_r4( statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.) &
1689 //new_line('')//'P_T = '&
1690 //str(ptr_PT_r4( statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.), &
1691 verb_opt=2)
1692 else
1693 call msg('calcPressure_gsv_nl (czp)', &
1694 new_line('')//'P_M = '&
1695 //str(ptr_PM_r8( statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.) &
1696 //new_line('')//'P_T = '&
1697 //str(ptr_PT_r8( statevector%myLonBeg,statevector%myLatBeg,:,1), vertical_opt=.false.), &
1698 verb_opt=2)
1699 end if
1700
1701 call msg('calcPressure_gsv_nl (czp)', 'END', verb_opt=2)
1702 call utl_tmg_stop(177)
1703 end subroutine calcPressure_gsv_nl
1704
1705 !---------------------------------------------------------
1706 ! czp_calcReturnPressure_gsv_nl
1707 !---------------------------------------------------------
1708 subroutine czp_calcReturnPressure_gsv_nl( statevector, &
1709 ZTin_r4_opt, ZMin_r4_opt, &
1710 ZTin_r8_opt, ZMin_r8_opt, &
1711 PTout_r4_opt, PMout_r4_opt, &
1712 PTout_r8_opt, PMout_r8_opt, &
1713 Ps_in_hPa_opt)
1714 !
1715 !:Purpose: Compute or retrieve pressures and return values in pointer arguments.
1716 ! Proceeds to vcode dispatching.
1717 !
1718 implicit none
1719
1720 ! Arguments:
1721 type(struct_gsv), intent(in) :: statevector
1722 real(4), optional, pointer, intent(in) :: ZTin_r4_opt(:,:,:,:)
1723 real(4), optional, pointer, intent(in) :: ZMin_r4_opt(:,:,:,:)
1724 real(8), optional, pointer, intent(in) :: ZTin_r8_opt(:,:,:,:)
1725 real(8), optional, pointer, intent(in) :: ZMin_r8_opt(:,:,:,:)
1726 real(4), optional, pointer, intent(inout) :: PTout_r4_opt(:,:,:,:)
1727 real(4), optional, pointer, intent(inout) :: PMout_r4_opt(:,:,:,:)
1728 real(8), optional, pointer, intent(inout) :: PTout_r8_opt(:,:,:,:)
1729 real(8), optional, pointer, intent(inout) :: PMout_r8_opt(:,:,:,:)
1730 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, conversion from hPa to mbar done for surface pressure
1731
1732 ! Locals:
1733 integer :: Vcode
1734
1735 call utl_tmg_start(177,'low-level--czp_calcPressure_nl')
1736 call msg('czp_calcReturnPressure_gsv_nl', 'START', verb_opt=2)
1737
1738 Vcode = gsv_getVco(statevector)%vcode
1739 if (Vcode == 5005 .or. Vcode == 5002 .or. Vcode == 5100) then
1740 if ( gsv_getDataKind(statevector) == 4 ) then
1741 if ( .not. (present(PTout_r4_opt) .and. present(PMout_r4_opt))) then
1742 call utl_abort('czp_calcReturnPressure_gsv_nl: dataKind=4: P{T,M}out_r4_opt expected')
1743 end if
1744 call calcPressure_gsv_nl_vcode5xxx_r4(statevector, PTout_r4_opt, PMout_r4_opt, &
1745 Ps_in_hPa_opt)
1746 else
1747 if ( .not. (present(PTout_r8_opt) .and. present(PMout_r8_opt))) then
1748 call utl_abort('czp_calcReturnPressure_gsv_nl: dataKind=8: P{T,M}out_r8_opt expected')
1749 end if
1750 call calcPressure_gsv_nl_vcode5xxx_r8(statevector, PTout_r8_opt, PMout_r8_opt, &
1751 Ps_in_hPa_opt)
1752 end if
1753 else if (Vcode == 21001) then
1754 if ( gsv_getDataKind(statevector) == 4 ) then
1755 if ( .not. (present(ZTin_r4_opt) .and. present(ZMin_r4_opt))) then
1756 call utl_abort('czp_calcReturnPressure_gsv_nl: dataKind=4: Z{T,M}out_r4_opt expected')
1757 end if
1758 if ( .not. (present(PTout_r4_opt) .and. present(PMout_r4_opt))) then
1759 call utl_abort('czp_calcReturnPressure_gsv_nl: dataKind=4: P{T,M}out_r4_opt expected')
1760 end if
1761 call calcPressure_gsv_nl_vcode2100x(statevector, &
1762 ZTin_r4_opt=ZTin_r4_opt, &
1763 ZMin_r4_opt=ZMin_r4_opt, &
1764 PTout_r4_opt=PTout_r4_opt, &
1765 PMout_r4_opt=PMout_r4_opt)
1766 else ! datakind = 8
1767 if ( .not. (present(ZTin_r8_opt) .and. present(ZMin_r8_opt))) then
1768 call utl_abort('czp_calcReturnPressure_gsv_nl: dataKind=8: Z{T,M}out_r8_opt expected')
1769 end if
1770 if ( .not. (present(PTout_r8_opt) .and. present(PMout_r8_opt))) then
1771 call utl_abort('czp_calcReturnPressure_gsv_nl: dataKind=8: P{T,M}out_r8_opt expected')
1772 end if
1773 call calcPressure_gsv_nl_vcode2100x(statevector, &
1774 ZTin_r8_opt=ZTin_r8_opt, &
1775 ZMin_r8_opt=ZMin_r8_opt, &
1776 PTout_r8_opt=PTout_r8_opt, &
1777 PMout_r8_opt=PMout_r8_opt)
1778 end if
1779 end if
1780
1781 call msg('czp_calcReturnPressure_gsv_nl', 'END', verb_opt=2)
1782 call utl_tmg_stop(177)
1783 end subroutine czp_calcReturnPressure_gsv_nl
1784
1785 !---------------------------------------------------------
1786 ! calcPressure_gsv_nl_vcode2100x
1787 !---------------------------------------------------------
1788 subroutine calcPressure_gsv_nl_vcode2100x(statevector, &
1789 ZTin_r4_opt, ZMin_r4_opt, &
1790 ZTin_r8_opt, ZMin_r8_opt, &
1791 PTout_r4_opt, PMout_r4_opt, &
1792 PTout_r8_opt, PMout_r8_opt)
1793 !
1794 !:Purpose: Compute pressure and return values in pointer arguments.
1795 ! GEM-H statevector input.
1796 !
1797 implicit none
1798
1799 ! Arguments:
1800 type(struct_gsv), intent(in) :: statevector
1801 real(4), pointer, optional, intent(in) :: ZTin_r4_opt(:,:,:,:)
1802 real(4), pointer, optional, intent(in) :: ZMin_r4_opt(:,:,:,:)
1803 real(8), pointer, optional, intent(in) :: ZTin_r8_opt(:,:,:,:)
1804 real(8), pointer, optional, intent(in) :: ZMin_r8_opt(:,:,:,:)
1805 real(4), pointer, optional, intent(inout) :: PTout_r4_opt(:,:,:,:)
1806 real(4), pointer, optional, intent(inout) :: PMout_r4_opt(:,:,:,:)
1807 real(8), pointer, optional, intent(inout) :: PTout_r8_opt(:,:,:,:)
1808 real(8), pointer, optional, intent(inout) :: PMout_r8_opt(:,:,:,:)
1809
1810 ! Locals:
1811 integer :: stepIndex, latIndex, lonIndex, numStep
1812 integer :: lev_M,lev_T,nlev_M,nlev_T,status
1813 real(4) :: heightSfcOffset_T_r4, heightSfcOffset_M_r4
1814 real(4) :: lat_4
1815 real(8) :: hu, tt, cmp, Rgh, P0, dh, tv0, rMt, Z_T, Z_M, Z_M1, logP
1816 real(8) :: sLat, cLat, lat_8
1817 real(8) :: ScaleFactorBottom
1818 real(8), allocatable :: tv(:), pressure_T(:), pressure_M(:)
1819 real(4), pointer :: height_T_ptr_r4(:,:,:,:)
1820 real(4), pointer :: height_M_ptr_r4(:,:,:,:)
1821 real(4), pointer :: hu_ptr_r4(:,:,:,:),tt_ptr_r4(:,:,:,:)
1822 real(4), pointer :: P_T_ptr_r4(:,:,:,:),P_M_ptr_r4(:,:,:,:)
1823 real(4), pointer :: P0_ptr_r4(:,:,:,:)
1824 real(8), pointer :: height_T_ptr_r8(:,:,:,:)
1825 real(8), pointer :: height_M_ptr_r8(:,:,:,:)
1826 real(8), pointer :: P_T_ptr_r8(:,:,:,:),P_M_ptr_r8(:,:,:,:)
1827 real(8), pointer :: P0_ptr_r8(:,:,:,:)
1828 real(8), pointer :: hu_ptr_r8(:,:,:,:),tt_ptr_r8(:,:,:,:)
1829 real(8), pointer :: HeightSfc_ptr_r8(:,:)
1830
1831 call msg('calcPressure_gsv_nl_vcode2100x (czp)', 'START', verb_opt=4)
1832
1833 nlev_T = gsv_getNumLev(statevector,'TH')
1834 nlev_M = gsv_getNumLev(statevector,'MM')
1835 numStep = statevector%numStep
1836
1837 allocate(tv(nlev_T))
1838
1839 if (nlev_T /= nlev_M) then
1840 call utl_abort('calcPressure_gsv_nl_vcode2100x: nlev_T is not equal to nlev_M!')
1841 end if
1842
1843 status = vgd_get( statevector%vco%vgrid, &
1844 key='DHM - height of the diagnostic level (m)', &
1845 value=heightSfcOffset_M_r4)
1846 status = vgd_get( statevector%vco%vgrid, &
1847 key='DHT - height of the diagnostic level (t)', &
1848 value=heightSfcOffset_T_r4)
1849 call msg('calcPressure_gsv_nl_vcode2100x (czp)', &
1850 'height offset for near-sfc momentum level is:'//str(heightSfcOffset_M_r4)//' meters'&
1851 //new_line('')//'height offset for near-sfc thermo level is:'//str(heightSfcOffset_T_r4)//' meters', &
1852 verb_opt=2, mpiAll_opt=.false.)
1853 if ( .not.statevector%addHeightSfcOffset ) then
1854 call msg('calcPressure_gsv_nl_vcode2100x (czp)', new_line('') &
1855 //'--------------------------------------------------------------------------'//new_line('')&
1856 //'BUT HEIGHT OFFSET REMOVED FOR DIAGNOSTIC LEVELS FOR BACKWARD COMPATIBILITY'//new_line('')&
1857 //'--------------------------------------------------------------------------', &
1858 verb_opt=2, mpiAll_opt=.false.)
1859 end if
1860
1861 allocate(pressure_T(nlev_T))
1862 allocate(pressure_M(nlev_M))
1863
1864 if ( statevector%dataKind == 4 ) then
1865 if ( .not. (present(ZTin_r4_opt) .and. present(ZMin_r4_opt))) then
1866 call utl_abort('calcPressure_gsv_nl_vcode2100x (czp): dataKind=4: Z{T,M}in_r4_opt expected')
1867 end if
1868 height_T_ptr_r4 => ZTin_r4_opt
1869 height_M_ptr_r4 => ZMin_r4_opt
1870
1871 if ( .not. (present(PTout_r4_opt) .and. present(PMout_r4_opt))) then
1872 call utl_abort('calcPressure_gsv_nl_vcode2100x (czp): dataKind=4: P{T,M}out_r4_opt expected')
1873 end if
1874 P_M_ptr_r4 => PMout_r4_opt
1875 P_T_ptr_r4 => PTout_r4_opt
1876 call gsv_getField(statevector,hu_ptr_r4,'HU')
1877 call gsv_getField(statevector,tt_ptr_r4,'TT')
1878 call gsv_getField(statevector,P0_ptr_r4,'P0')
1879
1880 ! initialize the pressure pointer to zero
1881 P_M_ptr_r4(:,:,:,:) = 0.0
1882 P_T_ptr_r4(:,:,:,:) = 0.0
1883 else ! datakind = 8
1884 if ( .not. (present(ZTin_r8_opt) .and. present(ZMin_r8_opt))) then
1885 call utl_abort('calcPressure_gsv_nl_vcode2100x (czp): dataKind=4: Z{T,M}in_r8_opt expected')
1886 end if
1887 height_T_ptr_r8 => ZTin_r8_opt
1888 height_M_ptr_r8 => ZMin_r8_opt
1889
1890 if ( .not. (present(PTout_r8_opt) .and. present(PMout_r8_opt))) then
1891 call utl_abort('calcPressure_gsv_nl_vcode2100x (czp): dataKind=8: P{T,M}_r8_opt expected')
1892 end if
1893 P_M_ptr_r8 => PMout_r8_opt
1894 P_T_ptr_r8 => PTout_r8_opt
1895 call gsv_getField(statevector,hu_ptr_r8,'HU')
1896 call gsv_getField(statevector,tt_ptr_r8,'TT')
1897 call gsv_getField(statevector,P0_ptr_r8,'P0')
1898
1899 ! initialize the pressure pointer to zero
1900 P_M_ptr_r8(:,:,:,:) = 0.0d0
1901 P_T_ptr_r8(:,:,:,:) = 0.0d0
1902 end if
1903 HeightSfc_ptr_r8 => gsv_getHeightSfc(statevector)
1904
1905 ! Development notes (@mad001)
1906 ! if feasible, consider reusing the same code for both
1907 ! `calcPressure_{gsv,col}_nl_vcode2100x`
1908 do_computePressure_gsv_nl: do stepIndex = 1, numStep
1909 do latIndex = statevector%myLatBeg, statevector%myLatEnd
1910 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
1911
1912 pressure_T(:) = 0.0D0
1913 pressure_M(:) = 0.0D0
1914
1915 ! latitude
1916 lat_4 = statevector%hco%lat2d_4(lonIndex,latIndex)
1917 lat_8 = real(lat_4,8)
1918 sLat = sin(lat_8)
1919 cLat = cos(lat_8)
1920
1921 if ( statevector%dataKind == 4 ) then
1922 P0 = real(P0_ptr_r4(lonIndex,latIndex,1, stepIndex),8)
1923 else
1924 P0 = P0_ptr_r8(lonIndex,latIndex,1, stepIndex)
1925 end if
1926
1927 rMT = HeightSfc_ptr_r8(lonIndex,latIndex)
1928
1929 ! compute pressure on diagnostic levels
1930 if ( statevector%dataKind == 4 ) then
1931 hu = real(hu_ptr_r4(lonIndex,latIndex,nlev_T,stepIndex),8)
1932 tt = real(tt_ptr_r4(lonIndex,latIndex,nlev_T,stepIndex),8)
1933 else
1934 hu = hu_ptr_r8(lonIndex,latIndex,nlev_T,stepIndex)
1935 tt = tt_ptr_r8(lonIndex,latIndex,nlev_T,stepIndex)
1936 end if
1937 tv0 = phf_fotvt8(tt,hu)
1938
1939 ! thermo diagnostic level
1940 if ( statevector%dataKind == 4 ) then
1941 Z_T = real(height_T_ptr_r4(lonIndex,latIndex,nlev_T,stepIndex),8)
1942 else
1943 Z_T = height_T_ptr_r8(lonIndex,latIndex,nlev_T,stepIndex)
1944 end if
1945 cmp = gpscompressibility(P0,tt,hu)
1946 tv(nlev_T) = tv0*cmp
1947 dh = Z_T - rMT
1948 Rgh = phf_gravityalt(sLat, rMT+0.5D0*dh)
1949 pressure_T(nlev_T) = P0*exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(nlev_T))
1950
1951 ! momentum diagnostic level
1952 if ( statevector%dataKind == 4 ) then
1953 Z_M = real(height_M_ptr_r4(lonIndex,latIndex,nlev_M,stepIndex),8)
1954 else
1955 Z_M = height_M_ptr_r8(lonIndex,latIndex,nlev_M,stepIndex)
1956 end if
1957 dh = Z_M - rMT
1958 Rgh = phf_gravityalt(sLat, rMT+0.5D0*dh)
1959 pressure_M(nlev_M) = P0*exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(nlev_T))
1960
1961 ! compute pressure on all levels above except the last
1962 do lev_M = nlev_M-1, 1, -1
1963 lev_T = lev_M ! thermo level just below
1964 if ( statevector%dataKind == 4 ) then
1965 hu = real(hu_ptr_r4(lonIndex,latIndex,lev_T,stepIndex),8)
1966 tt = real(tt_ptr_r4(lonIndex,latIndex,lev_T,stepIndex),8)
1967 Z_M = real(height_M_ptr_r4(lonIndex,latIndex,lev_M,stepIndex),8)
1968 Z_M1 = real(height_M_ptr_r4(lonIndex,latIndex,lev_M+1,stepIndex),8)
1969 Z_T = real(height_T_ptr_r4(lonIndex,latIndex,lev_T,stepIndex),8)
1970 else
1971 hu = hu_ptr_r8(lonIndex,latIndex,lev_T,stepIndex)
1972 tt = tt_ptr_r8(lonIndex,latIndex,lev_T,stepIndex)
1973 Z_M = height_M_ptr_r8(lonIndex,latIndex,lev_M,stepIndex)
1974 Z_M1 = height_M_ptr_r8(lonIndex,latIndex,lev_M+1,stepIndex)
1975 Z_T = height_T_ptr_r8(lonIndex,latIndex,lev_T,stepIndex)
1976 end if
1977 tv0 = phf_fotvt8(tt,hu)
1978 dh = Z_M - Z_M1
1979 Rgh = phf_gravityalt(sLat, Z_M1+0.5D0*dh)
1980
1981 ! approximation of tv from pressure on previous momentum level
1982 cmp = gpscompressibility(pressure_M(lev_M+1),tt,hu)
1983 tv(lev_T) = tv0*cmp
1984 pressure_M(lev_M) = pressure_M(lev_M+1) * &
1985 exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(lev_T))
1986 ! first interpolation of thermo pressure
1987 scaleFactorBottom = (Z_T-Z_M1)/(Z_M-Z_M1)
1988 logP = (1.0D0-scaleFactorBottom)*log(pressure_M(lev_M+1)) + &
1989 scaleFactorBottom*log(pressure_M(lev_M))
1990 pressure_T(lev_T) = exp(logP)
1991
1992 ! second iteration on tv
1993 cmp = gpscompressibility(pressure_T(lev_T),tt,hu)
1994 tv(lev_T) = tv0*cmp
1995 pressure_M(lev_M) = pressure_M(lev_M+1) * &
1996 exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(lev_T))
1997
1998 ! second iteration interpolation of thermo pressure
1999 logP = (1.0D0-scaleFactorBottom)*log(pressure_M(lev_M+1)) + &
2000 scaleFactorBottom*log(pressure_M(lev_M))
2001 pressure_T(lev_T) = exp(logP)
2002
2003 end do
2004
2005 ! fill the height array
2006 if ( statevector%dataKind == 4 ) then
2007 do lev_T = 1, nlev_T
2008 P_T_ptr_r4(lonIndex,latIndex,lev_T,stepIndex) = &
2009 real(pressure_T(lev_T),4)
2010 end do
2011 do lev_M = 1, nlev_M
2012 P_M_ptr_r4(lonIndex,latIndex,lev_M,stepIndex) = &
2013 real(pressure_M(lev_M),4)
2014 end do
2015 else
2016 P_T_ptr_r8(lonIndex,latIndex,1:nlev_T,stepIndex)=pressure_T(1:nlev_T)
2017 P_M_ptr_r8(lonIndex,latIndex,1:nlev_M,stepIndex)=pressure_M(1:nlev_M)
2018 end if
2019
2020 end do ! lonIndex
2021 end do ! latIndex
2022 end do do_computePressure_gsv_nl
2023
2024 deallocate(pressure_T)
2025 deallocate(pressure_M)
2026 deallocate(tv)
2027
2028 call msg('calcPressure_gsv_nl_vcode2100x (czp)', 'END', verb_opt=4)
2029 end subroutine calcPressure_gsv_nl_vcode2100x
2030
2031 !---------------------------------------------------------
2032 ! calcPressure_gsv_nl_vcode5xxx_r8
2033 !---------------------------------------------------------
2034 subroutine calcPressure_gsv_nl_vcode5xxx_r8(statevector, P_T, P_M, Ps_in_hPa_opt)
2035 !
2036 !:Purpose: Pressure retrieval for GEM-P real(8) statevector, values
2037 ! values returned in pointers.
2038 !
2039 implicit none
2040
2041 ! Arguments:
2042 type(struct_gsv), intent(in) :: statevector
2043 real(8), pointer, intent(inout) :: P_T(:,:,:,:)
2044 real(8), pointer, intent(inout) :: P_M(:,:,:,:)
2045 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, conversion from hPa to mbar done for surface pressure
2046
2047 ! Locals:
2048 real(kind=8), allocatable :: Psfc(:,:), PsfcLS(:,:)
2049 real(kind=8), pointer :: PressureM_out(:,:,:), PressureT_out(:,:,:)
2050 real(kind=8), pointer :: field_Psfc(:,:,:,:), field_PsfcLS(:,:,:,:)
2051 integer :: stepIndex, numStep, Vcode
2052
2053 call msg('calcPressure_gsv_nl_vcode5xxx_r8 (czp)', 'START', verb_opt=4)
2054
2055 Vcode = gsv_getVco(statevector)%vcode
2056
2057 allocate(Psfc(statevector%myLonBeg:statevector%myLonEnd, &
2058 statevector%myLatBeg:statevector%myLatEnd))
2059 call gsv_getField(statevector,field_Psfc,'P0')
2060
2061 if (Vcode == 5100) then
2062 allocate(PsfcLS(statevector%myLonBeg:statevector%myLonEnd, &
2063 statevector%myLatBeg:statevector%myLatEnd))
2064 call gsv_getField(statevector,field_PsfcLS,'P0LS')
2065 end if
2066
2067 numStep = statevector%numStep
2068
2069 do stepIndex = 1, numStep
2070 Psfc(:,:) = field_Psfc(:,:,1,stepIndex)
2071 if ( present(Ps_in_hPa_opt) ) then
2072 if ( Ps_in_hPa_opt ) Psfc = Psfc * mpc_pa_per_mbar_r8
2073 end if
2074
2075 if (Vcode == 5100) then
2076 PsfcLS(:,:) = field_PsfcLS(:,:,1,stepIndex)
2077 if ( present(Ps_in_hPa_opt) ) then
2078 if ( Ps_in_hPa_opt ) PsfcLS = PsfcLS * mpc_pa_per_mbar_r8
2079 end if
2080 call fetch3DLevels_r8(statevector%vco, Psfc, sfcFldLS_opt=PsfcLS, &
2081 fldM_opt=PressureM_out, fldT_opt=PressureT_out)
2082 else
2083 call fetch3DLevels_r8(statevector%vco, Psfc, &
2084 fldM_opt=PressureM_out, fldT_opt=PressureT_out)
2085 end if
2086 P_M(:,:,:,stepIndex) = PressureM_out(:,:,:)
2087 P_T(:,:,:,stepIndex) = PressureT_out(:,:,:)
2088 deallocate(PressureM_out, PressureT_out)
2089
2090 end do
2091
2092 deallocate(Psfc)
2093 if (Vcode == 5100) deallocate(PsfcLS)
2094
2095 call msg('calcPressure_gsv_nl_vcode5xxx_r8 (czp)', 'END', verb_opt=4)
2096 end subroutine calcPressure_gsv_nl_vcode5xxx_r8
2097
2098 !---------------------------------------------------------
2099 ! calcPressure_gsv_nl_vcode5xxx_r4
2100 !---------------------------------------------------------
2101 subroutine calcPressure_gsv_nl_vcode5xxx_r4(statevector, P_T, P_M, Ps_in_hPa_opt)
2102 !
2103 !:Purpose: Pressure retrieval for GEM-P real(4) statevector, values
2104 ! values returned in pointers.
2105 !
2106 implicit none
2107
2108 ! Arguments:
2109 type(struct_gsv), intent(in) :: statevector
2110 real(4), pointer, intent(inout) :: P_T(:,:,:,:)
2111 real(4), pointer, intent(inout) :: P_M(:,:,:,:)
2112 logical, optional, intent(in) :: Ps_in_hPa_opt ! If true, conversion from hPa to mbar done for surface pressure
2113
2114 ! Locals:
2115 real(kind=4), allocatable :: Psfc(:,:), PsfcLS(:,:)
2116 real(kind=4), pointer :: PressureM_out(:,:,:), PressureT_out(:,:,:)
2117 real(kind=4), pointer :: field_Psfc(:,:,:,:), field_PsfcLS(:,:,:,:)
2118 integer :: stepIndex, numStep, Vcode
2119
2120 call msg('calcPressure_gsv_nl_vcode5xxx_r4 (czp)', 'START', verb_opt=4)
2121
2122 Vcode = gsv_getVco(statevector)%vcode
2123
2124 allocate(Psfc(statevector%myLonBeg:statevector%myLonEnd, &
2125 statevector%myLatBeg:statevector%myLatEnd))
2126 call gsv_getField(statevector,field_Psfc,'P0')
2127
2128 if (Vcode == 5100) then
2129 allocate(PsfcLS(statevector%myLonBeg:statevector%myLonEnd, &
2130 statevector%myLatBeg:statevector%myLatEnd))
2131 call gsv_getField(statevector,field_PsfcLS,'P0LS')
2132 end if
2133
2134 numStep = statevector%numStep
2135
2136 do stepIndex = 1, numStep
2137 Psfc(:,:) = field_Psfc(:,:,1,stepIndex)
2138 if ( present(Ps_in_hPa_opt) ) then
2139 if ( Ps_in_hPa_opt ) Psfc = Psfc * mpc_pa_per_mbar_r4
2140 end if
2141
2142 if (Vcode == 5100) then
2143 PsfcLS(:,:) = field_PsfcLS(:,:,1,stepIndex)
2144 if ( present(Ps_in_hPa_opt) ) then
2145 if ( Ps_in_hPa_opt ) PsfcLS = PsfcLS * mpc_pa_per_mbar_r4
2146 end if
2147
2148 call fetch3DLevels_r4(statevector%vco, Psfc, sfcFldLS_opt=PsfcLS, &
2149 fldM_opt=PressureM_out, fldT_opt=PressureT_out)
2150 else
2151 call fetch3DLevels_r4(statevector%vco, Psfc, &
2152 fldM_opt=PressureM_out, fldT_opt=PressureT_out)
2153 end if
2154 P_M(:,:,:,stepIndex) = PressureM_out(:,:,:)
2155 P_T(:,:,:,stepIndex) = PressureT_out(:,:,:)
2156 deallocate(PressureM_out, PressureT_out)
2157
2158 end do
2159
2160 deallocate(Psfc)
2161 if (Vcode == 5100) deallocate(PsfcLS)
2162
2163 call msg('calcPressure_gsv_nl_vcode5xxx_r4 (czp)', 'START', verb_opt=4)
2164 end subroutine calcPressure_gsv_nl_vcode5xxx_r4
2165
2166 !---------------------------------------------------------
2167 ! calcPressure_gsv_tl
2168 !---------------------------------------------------------
2169 subroutine calcPressure_gsv_tl( statevector, statevectorRef)
2170 !
2171 !:Purpose: Tangent of pressure computation.
2172 !
2173 implicit none
2174
2175 ! Arguments:
2176 type(struct_gsv), intent(inout) :: statevector ! statevector that will contain the P_T/P_M increments
2177 type(struct_gsv), intent(in) :: statevectorRef ! statevector containing needed reference fields
2178
2179 ! Locals:
2180 integer :: Vcode
2181
2182 call utl_tmg_start(178,'low-level--czp_calcPressure_tl')
2183 call msg('calcPressure_gsv_tl (czp)', 'START', verb_opt=2)
2184
2185 Vcode = gsv_getVco(statevectorRef)%vcode
2186 if (Vcode == 5005 .or. Vcode == 5002) then
2187 if ( .not. gsv_varExist(statevector,'P_*') ) then
2188 call utl_abort('calcPressure_gsv_tl (czp): for vcode 5xxx, variables P_T and P_M must be allocated in gridstatevector')
2189 end if
2190 if ( .not. gsv_varExist(statevector,'P0') ) then
2191 call utl_abort('calcPressure_gsv_tl (czp): for vcode 5xxx, variable P0 must be allocated in gridstatevector')
2192 end if
2193 call calcPressure_gsv_tl_vcode5xxx
2194 else if (Vcode == 21001) then
2195 ! Development notes (@mad001)
2196 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
2197 call calcPressure_gsv_tl_vcode2100x
2198 else
2199 call utl_abort('calcPressure_gsv_tl (czp): not implemented')
2200 end if
2201
2202 call msg('calcPressure_gsv_tl (czp)', 'END', verb_opt=2)
2203 call utl_tmg_stop(178)
2204
2205 contains
2206
2207 !---------------------------------------------------------
2208 ! calcPressure_gsv_tl_vcode2100x
2209 !---------------------------------------------------------
2210 subroutine calcPressure_gsv_tl_vcode2100x
2211 implicit none
2212
2213 call utl_abort('calcPressure_gsv_tl (czp): vcode 21001 not implemented yet')
2214
2215 end subroutine calcPressure_gsv_tl_vcode2100x
2216
2217 !---------------------------------------------------------
2218 ! calcPressure_gsv_tl_vcode5xxx
2219 !---------------------------------------------------------
2220 subroutine calcPressure_gsv_tl_vcode5xxx
2221 implicit none
2222
2223 ! Locals:
2224 real(8), allocatable :: Psfc(:,:)
2225 real(4), pointer :: delPsfc_r4(:,:,:,:)
2226 real(8), pointer :: delPsfc_r8(:,:,:,:)
2227 real(8), pointer :: field_Psfc(:,:,:,:)
2228 real(4), pointer :: delP_T_r4(:,:,:,:)
2229 real(8), pointer :: delP_T_r8(:,:,:,:)
2230 real(4), pointer :: delP_M_r4(:,:,:,:)
2231 real(8), pointer :: delP_M_r8(:,:,:,:)
2232 real(8), pointer :: dP_dPsfc_T(:,:,:)
2233 real(8), pointer :: dP_dPsfc_M(:,:,:)
2234 integer :: status, stepIndex,lonIndex,latIndex
2235 integer :: lev_M, lev_T, nlev_T, nlev_M, numStep
2236
2237 call msg('calcPressure_gsv_tl_vcode5xxx (czp)', 'START', verb_opt=4)
2238
2239 nullify(dP_dPsfc_T)
2240 nullify(dP_dPsfc_M)
2241 nullify(delPsfc_r4,delPsfc_r8)
2242 nullify(delP_T_r4,delP_T_r8)
2243 nullify(delP_M_r4,delP_M_r8)
2244
2245 if (gsv_getDataKind(statevector) == 4) then
2246 call gsv_getField(statevector,delP_T_r4,'P_T')
2247 call gsv_getField(statevector,delP_M_r4,'P_M')
2248 call gsv_getField(statevector,delPsfc_r4,'P0')
2249 else
2250 call gsv_getField(statevector,delP_T_r8,'P_T')
2251 call gsv_getField(statevector,delP_M_r8,'P_M')
2252 call gsv_getField(statevector,delPsfc_r8,'P0')
2253 end if
2254
2255 nullify(field_Psfc)
2256 call gsv_getField(statevectorRef,field_Psfc,'P0')
2257
2258 nlev_T = gsv_getNumLev(statevector,'TH')
2259 nlev_M = gsv_getNumLev(statevector,'MM')
2260 numStep = statevector%numstep
2261
2262 allocate(Psfc(statevector%myLonBeg:statevector%myLonEnd, &
2263 statevector%myLatBeg:statevector%myLatEnd))
2264
2265 do stepIndex = 1, numStep
2266
2267 Psfc(:,:) = field_Psfc(:,:,1,stepIndex)
2268
2269 ! dP_dPsfc_M
2270 nullify(dP_dPsfc_M)
2271 status = vgd_dpidpis(statevector%vco%vgrid, &
2272 statevector%vco%ip1_M, &
2273 dP_dPsfc_M, &
2274 Psfc)
2275 if( status .ne. VGD_OK ) then
2276 call utl_abort('calcPressure_gsv_tl (czp): ERROR with vgd_dpidpis')
2277 end if
2278 ! calculate delP_M
2279 if (gsv_getDataKind(statevector) == 4) then
2280 do lev_M = 1, nlev_M
2281 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2282 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2283 delP_M_r4(lonIndex,latIndex,lev_M,stepIndex) = &
2284 dP_dPsfc_M(lonIndex-statevector%myLonBeg+1,&
2285 latIndex-statevector%myLatBeg+1,lev_M) * &
2286 delPsfc_r4(lonIndex,latIndex,1,stepIndex)
2287 end do
2288 end do
2289 end do
2290 else
2291 do lev_M = 1, nlev_M
2292 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2293 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2294 delP_M_r8(lonIndex,latIndex,lev_M,stepIndex) = &
2295 dP_dPsfc_M(lonIndex-statevector%myLonBeg+1,&
2296 latIndex-statevector%myLatBeg+1,lev_M) * &
2297 delPsfc_r8(lonIndex,latIndex,1,stepIndex)
2298 end do
2299 end do
2300 end do
2301 end if
2302 deallocate(dP_dPsfc_M)
2303
2304 ! dP_dPsfc_T
2305 nullify(dP_dPsfc_T)
2306 status = vgd_dpidpis(statevector%vco%vgrid, &
2307 statevector%vco%ip1_T, &
2308 dP_dPsfc_T, &
2309 Psfc)
2310 if( status .ne. VGD_OK ) then
2311 call utl_abort('calcPressure_gsv_tl (czp): ERROR with vgd_dpidpis')
2312 end if
2313 ! calculate delP_T
2314 if (gsv_getDataKind(statevector) == 4) then
2315 do lev_T = 1, nlev_T
2316 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2317 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2318 delP_T_r4(lonIndex,latIndex,lev_T,stepIndex) = &
2319 dP_dPsfc_T(lonIndex-statevector%myLonBeg+1,&
2320 latIndex-statevector%myLatBeg+1,lev_T) * &
2321 delPsfc_r4(lonIndex,latIndex,1,stepIndex)
2322 end do
2323 end do
2324 end do
2325 else
2326 do lev_T = 1, nlev_T
2327 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2328 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2329 delP_T_r8(lonIndex,latIndex,lev_T,stepIndex) = &
2330 dP_dPsfc_T(lonIndex-statevector%myLonBeg+1,&
2331 latIndex-statevector%myLatBeg+1,lev_T) * &
2332 delPsfc_r8(lonIndex,latIndex,1,stepIndex)
2333 end do
2334 end do
2335 end do
2336 end if
2337 deallocate(dP_dPsfc_T)
2338
2339 end do
2340
2341 deallocate(Psfc)
2342
2343 call msg('calcPressure_gsv_tl_vcode5xxx (czp)', 'END', verb_opt=4)
2344 end subroutine calcPressure_gsv_tl_vcode5xxx
2345
2346 end subroutine calcPressure_gsv_tl
2347
2348 !---------------------------------------------------------
2349 ! calcPressure_gsv_ad
2350 !---------------------------------------------------------
2351 subroutine calcPressure_gsv_ad( statevector, statevectorRef)
2352 !
2353 !:Purpose: Adjoint of pressure computation.
2354 !
2355 implicit none
2356
2357 ! Arguments:
2358 type(struct_gsv), intent(inout) :: statevector ! statevector that will contain increment of P_T/P_M
2359 type(struct_gsv), intent(in) :: statevectorRef ! statevector containing needed reference fields
2360
2361 ! Locals:
2362 integer :: Vcode
2363
2364 call utl_tmg_start(179,'low-level--czp_calcPressure_ad')
2365 call msg('calcPressure_gsv_ad (czp)', 'START', verb_opt=2)
2366
2367 Vcode = gsv_getVco(statevectorRef)%vcode
2368 if (Vcode == 5005 .or. Vcode == 5002) then
2369 if ( .not. gsv_varExist(statevector,'P_*') ) then
2370 call utl_abort('calcPressure_gsv_ad (czp): for vcode 5xxx, variables P_M and P_T must be allocated in gridstatevector')
2371 end if
2372 if ( .not. gsv_varExist(statevector,'P0') ) then
2373 call utl_abort('calcPressure_gsv_ad (czp): for vcode 5xxx, variable P0 must be allocated in gridstatevector')
2374 end if
2375 call calcPressure_gsv_ad_vcode5xxx
2376 else if (Vcode == 21001) then
2377 ! Development notes (@mad001)
2378 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
2379 call calcPressure_gsv_ad_vcode2100x
2380 else
2381 call utl_abort('calcPressure_gsv_ad (czp): not implemented')
2382 end if
2383
2384 call msg('calcPressure_gsv_ad (czp)', 'END', verb_opt=2)
2385 call utl_tmg_stop(179)
2386
2387 contains
2388
2389 !---------------------------------------------------------
2390 ! calcPressure_gsv_ad_vcode2100x
2391 !---------------------------------------------------------
2392 subroutine calcPressure_gsv_ad_vcode2100x
2393 implicit none
2394
2395 call utl_abort('calcPressure_gsv_ad (czp): vcode 21001 not implemented yet')
2396
2397 end subroutine calcPressure_gsv_ad_vcode2100x
2398
2399 !---------------------------------------------------------
2400 ! calcPressure_gsv_ad_vcode5xxx
2401 !---------------------------------------------------------
2402 subroutine calcPressure_gsv_ad_vcode5xxx
2403 implicit none
2404
2405 ! Locals:
2406 real(8), allocatable :: Psfc(:,:)
2407 real(4), pointer :: delPsfc_r4(:,:,:,:)
2408 real(8), pointer :: delPsfc_r8(:,:,:,:)
2409 real(8), pointer :: field_Psfc(:,:,:,:)
2410 real(4), pointer :: delP_T_r4(:,:,:,:)
2411 real(8), pointer :: delP_T_r8(:,:,:,:)
2412 real(4), pointer :: delP_M_r4(:,:,:,:)
2413 real(8), pointer :: delP_M_r8(:,:,:,:)
2414 real(8), pointer :: dP_dPsfc_T(:,:,:)
2415 real(8), pointer :: dP_dPsfc_M(:,:,:)
2416 integer :: status, stepIndex,lonIndex,latIndex
2417 integer :: lev_M, lev_T, nlev_T, nlev_M, numStep
2418
2419 call msg('calcPressure_gsv_ad_vcode5xxx (czp)', 'START', verb_opt=4)
2420
2421 nullify(delPsfc_r4, delPsfc_r8)
2422 nullify(field_Psfc)
2423 nullify(delP_T_r4, delP_T_r8)
2424 nullify(delP_M_r4, delP_M_r8)
2425 nullify(dP_dPsfc_T)
2426 nullify(dP_dPsfc_M)
2427
2428 if (gsv_getDataKind(statevector) == 4) then
2429 call gsv_getField(statevector,delP_T_r4,'P_T')
2430 call gsv_getField(statevector,delP_M_r4,'P_M')
2431 call gsv_getField(statevector,delPsfc_r4,'P0')
2432 else
2433 call gsv_getField(statevector,delP_T_r8,'P_T')
2434 call gsv_getField(statevector,delP_M_r8,'P_M')
2435 call gsv_getField(statevector,delPsfc_r8,'P0')
2436 end if
2437 call gsv_getField(statevectorRef,field_Psfc,'P0')
2438
2439 nlev_T = gsv_getNumLev(statevector,'TH')
2440 nlev_M = gsv_getNumLev(statevector,'MM')
2441 numStep = statevector%numstep
2442
2443 allocate(Psfc(statevector%myLonBeg:statevector%myLonEnd, &
2444 statevector%myLatBeg:statevector%myLatEnd))
2445
2446 do stepIndex = 1, numStep
2447
2448 Psfc(:,:) = field_Psfc(:,:,1,stepIndex)
2449
2450 ! dP_dPsfc_M
2451 nullify(dP_dPsfc_M)
2452 status = vgd_dpidpis(statevector%vco%vgrid, &
2453 statevector%vco%ip1_M, &
2454 dP_dPsfc_M, &
2455 Psfc)
2456 if( status .ne. VGD_OK ) then
2457 call utl_abort('calcPressure_gsv_ad (czp): ERROR with vgd_dpidpis')
2458 end if
2459 ! calculate delP_M
2460 if (gsv_getDataKind(statevector) == 4) then
2461 do lev_M = 1, nlev_M
2462 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2463 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2464 delPsfc_r4(lonIndex,latIndex,1,stepIndex) = &
2465 delPsfc_r4(lonIndex,latIndex,1,stepIndex) + &
2466 dP_dPsfc_M(lonIndex-statevector%myLonBeg+1,&
2467 latIndex-statevector%myLatBeg+1,lev_M) * &
2468 delP_M_r4(lonIndex,latIndex,lev_M,stepIndex)
2469 end do
2470 end do
2471 end do
2472 else
2473 do lev_M = 1, nlev_M
2474 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2475 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2476 delPsfc_r8(lonIndex,latIndex,1,stepIndex) = &
2477 delPsfc_r8(lonIndex,latIndex,1,stepIndex) + &
2478 dP_dPsfc_M(lonIndex-statevector%myLonBeg+1,&
2479 latIndex-statevector%myLatBeg+1,lev_M) * &
2480 delP_M_r8(lonIndex,latIndex,lev_M,stepIndex)
2481 end do
2482 end do
2483 end do
2484 end if
2485 deallocate(dP_dPsfc_M)
2486
2487 ! dP_dPsfc_T
2488 nullify(dP_dPsfc_T)
2489 status = vgd_dpidpis(statevector%vco%vgrid, &
2490 statevector%vco%ip1_T, &
2491 dP_dPsfc_T, &
2492 Psfc)
2493 if( status .ne. VGD_OK ) then
2494 call utl_abort('calcPressure_gsv_ad (czp): ERROR with vgd_dpidpis')
2495 end if
2496 ! calculate delP_T
2497 if (gsv_getDataKind(statevector) == 4) then
2498 do lev_T = 1, nlev_T
2499 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2500 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2501 delPsfc_r4(lonIndex,latIndex,1,stepIndex) = &
2502 delPsfc_r4(lonIndex,latIndex,1,stepIndex) + &
2503 dP_dPsfc_T(lonIndex-statevector%myLonBeg+1,&
2504 latIndex-statevector%myLatBeg+1,lev_T) * &
2505 delP_T_r4(lonIndex,latIndex,lev_T,stepIndex)
2506 end do
2507 end do
2508 end do
2509 else
2510 do lev_T = 1, nlev_T
2511 do latIndex = statevector%myLatBeg, statevector%myLatEnd
2512 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
2513 delPsfc_r8(lonIndex,latIndex,1,stepIndex) = &
2514 delPsfc_r8(lonIndex,latIndex,1,stepIndex) + &
2515 dP_dPsfc_T(lonIndex-statevector%myLonBeg+1,&
2516 latIndex-statevector%myLatBeg+1,lev_T) * &
2517 delP_T_r8(lonIndex,latIndex,lev_T,stepIndex)
2518 end do
2519 end do
2520 end do
2521 end if
2522 deallocate(dP_dPsfc_T)
2523
2524 end do
2525
2526 deallocate(Psfc)
2527
2528 call msg('calcPressure_gsv_ad_vcode5xxx (czp)', 'END', verb_opt=4)
2529 end subroutine calcPressure_gsv_ad_vcode5xxx
2530
2531 end subroutine calcPressure_gsv_ad
2532
2533 !---------------------------------------------------------------------
2534 ! subroutines operating on struct_columnData
2535 !---------------------------------------------------------------------
2536
2537 !---------------------------------------------------------
2538 ! calcZandP_col_nl
2539 !---------------------------------------------------------
2540 subroutine calcZandP_col_nl(column)
2541 !
2542 !:Purpose: Compute pressure and height in the column in proper order
2543 ! depending on the vgrid kind
2544 !
2545 implicit none
2546
2547 ! Arguments:
2548 type(struct_columnData), intent(inout) :: column ! column that will contain the Z_*/P_* fields
2549
2550 ! Locals:
2551 integer :: Vcode
2552
2553 call msg('calcZandP_col_nl (czp)', 'START', verb_opt=2)
2554
2555 Vcode = column%vco%vcode
2556 if (Vcode == 5002 .or. Vcode == 5005 .or. Vcode == 5100) then
2557 ! if P_T, P_M not allocated : do nothing
2558 if (col_varExist(column,'P_*')) then
2559 call calcPressure_col_nl(column)
2560 if (col_varExist(column,'Z_*')) then
2561 call calcHeight_col_nl(column)
2562 end if
2563 end if
2564 else if (Vcode == 21001) then
2565 ! if Z_T, Z_M not allocated : do nothing
2566 if (col_varExist(column,'Z_*')) then
2567 call calcHeight_col_nl(column)
2568 if (col_varExist(column,'P_*')) then
2569 call calcPressure_col_nl(column)
2570 end if
2571 end if
2572 end if
2573
2574 call msg('calcZandP_col_nl (czp)', 'END', verb_opt=2)
2575 end subroutine calcZandP_col_nl
2576
2577 !---------------------------------------------------------
2578 ! calcZandP_col_tl
2579 !---------------------------------------------------------
2580 subroutine calcZandP_col_tl(columnInc, columnIncRef)
2581 !
2582 !:Purpose: Compute pressure and height increment in the column in proper
2583 ! order depending on the vgrid kind
2584 !
2585 implicit none
2586
2587 ! Arguments:
2588 type(struct_columnData), intent(inout) :: columnInc ! column that will contain the Z_*/P_* increments
2589 type(struct_columnData), intent(in) :: columnIncRef ! column containing needed reference fields
2590
2591 ! Locals:
2592 integer :: Vcode
2593
2594 call msg('calcZandP_col_tl (czp)', 'START', verb_opt=2)
2595
2596 Vcode = columnInc%vco%vcode
2597 if (Vcode == 5002 .or. Vcode == 5005) then
2598 ! if P_T, P_M not allocated : do nothing
2599 if (col_varExist(columnInc,'P_*')) then
2600 call calcPressure_col_tl( columnInc, columnIncRef)
2601 if (col_varExist(columnInc,'Z_*')) then
2602 call calcHeight_col_tl(columnInc, columnIncRef)
2603 end if
2604 end if
2605 else if (Vcode == 21001) then
2606 ! if Z_T, Z_M not allocated : do nothing
2607 if (col_varExist(columnInc,'Z_*')) then
2608 call calcHeight_col_tl(columnInc, columnIncRef)
2609
2610 if (col_varExist(columnInc,'P_*')) then
2611 call calcPressure_col_tl( columnInc, columnIncRef)
2612 end if
2613 end if
2614 else
2615 call utl_abort('calcZandP_col_tl (czp): not implemented')
2616 end if
2617
2618 call msg('calcZandP_col_tl (czp)', 'END', verb_opt=2)
2619 end subroutine calcZandP_col_tl
2620
2621 !---------------------------------------------------------
2622 ! calcZandP_col_ad
2623 !---------------------------------------------------------
2624 subroutine calcZandP_col_ad(columnInc, columnIncRef)
2625 !
2626 !:Purpose: Adjoint of pressure and height increment computation in the
2627 ! column in proper order depending on the vgrid kind
2628 !
2629 implicit none
2630
2631 ! Arguments:
2632 type(struct_columnData), intent(inout) :: columnInc ! column that will contain the Z_*/P_* increments
2633 type(struct_columnData), intent(in) :: columnIncRef ! column containing needed reference fields
2634
2635 ! Locals:
2636 integer :: Vcode
2637
2638 call msg('calcZandP_col_ad (czp)', 'START', verb_opt=2)
2639
2640 Vcode = columnInc%vco%vcode
2641 if (Vcode == 5002 .or. Vcode == 5005) then
2642 ! if Z_T, Z_M not allocated : do nothing
2643 if (col_varExist(columnInc,'Z_*')) then
2644 call calcHeight_col_ad(columnInc, columnIncRef)
2645 if (col_varExist(columnInc,'P_*')) then
2646 call calcPressure_col_ad( columnInc, columnIncRef)
2647 end if
2648 end if
2649 else if (Vcode == 21001) then
2650 ! if P_T, P_M not allocated : do nothing
2651 if (col_varExist(columnInc,'P_*')) then
2652 call calcPressure_col_ad( columnInc, columnIncRef)
2653 if (col_varExist(columnInc,'Z_*')) then
2654 call calcHeight_col_ad(columnInc, columnIncRef)
2655 end if
2656 end if
2657 else
2658 call utl_abort('calcZandP_col_ad (czp): not implelmented')
2659 end if
2660
2661 call msg('calcZandP_col_ad (czp)', 'END', verb_opt=2)
2662 end subroutine calcZandP_col_ad
2663
2664 !---------------------------------------------------------
2665 ! calcHeight_col_nl
2666 !---------------------------------------------------------
2667 subroutine calcHeight_col_nl(column)
2668 !
2669 !:Purpose: Compute or retrieve heights on the column.
2670 !
2671 implicit none
2672
2673 ! Arguments:
2674 type(struct_columnData), intent(inout) :: column ! column that will contain the Z_M/Z_T fields
2675
2676 ! Locals:
2677 real(8), pointer :: Z_T(:,:), Z_M(:,:)
2678
2679 call msg('calcHeight_col_nl (czp)', 'START', verb_opt=2)
2680
2681 Z_T => col_getAllColumns(column, 'Z_T')
2682 Z_M => col_getAllColumns(column, 'Z_M')
2683 call czp_calcReturnHeight_col_nl(column, Z_T, Z_M)
2684
2685 call msg('calcHeight_col_nl (czp)', &
2686 new_line('')//'Z_M = '//str(col_getColumn(column,1,'Z_M')) &
2687 //new_line('')//'Z_T = '//str(col_getColumn(column,1,'Z_T')), &
2688 verb_opt=2)
2689
2690 call msg('calcHeight_col_nl (czp)', 'END', verb_opt=2)
2691 end subroutine calcHeight_col_nl
2692
2693 !---------------------------------------------------------
2694 ! czp_calcReturnHeight_col_nl
2695 !---------------------------------------------------------
2696 subroutine czp_calcReturnHeight_col_nl(column, Z_T, Z_M)
2697 !
2698 !:Purpose: Return heights on the column.
2699 !
2700 implicit none
2701
2702 ! Arguments:
2703 type(struct_columnData), intent(in) :: column ! reference column containing temperature and geopotential
2704 real(8), pointer, intent(inout) :: Z_T(:,:) ! computed column height values on thermodynamic levels
2705 real(8), pointer, intent(inout) :: Z_M(:,:) ! computed column height values on momentum levels
2706
2707 ! Locals:
2708 integer :: vcode
2709
2710 call msg('czp_calcReturnHeight_col_nl (czp)', 'START', verb_opt=2)
2711
2712 Vcode = col_getVco(column)%vcode
2713 if (Vcode == 5005 .or. Vcode == 5002 .or. Vcode == 5100) then
2714 if ( .not. (col_varExist(column,'P0') .and. col_varExist(column,'TT') .and. &
2715 col_varExist(column,'HU')) ) then
2716 call utl_abort('czp_calcReturnHeight_col_nl: for vcode 5xxx, variables P0, TT and HU must be allocated in column')
2717 end if
2718 call calcHeight_col_nl_vcode5xxx(column, Z_T, Z_M)
2719 else if (Vcode == 21001) then
2720 call calcHeight_col_nl_vcode2100x(column, Z_T, Z_M)
2721
2722 end if
2723
2724 call msg('czp_calcReturnHeight_col_nl (czp)', 'END', verb_opt=2)
2725 end subroutine czp_calcReturnHeight_col_nl
2726
2727 !---------------------------------------------------------
2728 ! calcHeight_col_nl_vcode2100x
2729 !---------------------------------------------------------
2730 subroutine calcHeight_col_nl_vcode2100x(column, Z_T, Z_M)
2731 !
2732 !:Purpose: Return heights on a GEM-H column.
2733 !
2734 implicit none
2735
2736 ! Arguments:
2737 type(struct_columnData), intent(in) :: column ! reference column containing temperature and geopotential
2738 real(8), pointer, intent(inout) :: Z_T(:,:) ! computed column height values on thermodynamic levels
2739 real(8), pointer, intent(inout) :: Z_M(:,:) ! computed column height values on momentum levels
2740
2741 ! Locals:
2742 real(8), allocatable :: hSfc(:,:)
2743 real(8), pointer :: hPtrM(:,:,:), hPtrT(:,:,:)
2744 integer :: numCol, colIndex
2745
2746 call msg('calcHeight_col_nl_vcode2100x (czp)', 'START', verb_opt=4)
2747 if ( col_getNumCol(column) <= 0 ) then
2748 call msg('calcHeight_col_nl_vcode2100x (czp)',&
2749 'END (number of columns <= 0)', verb_opt=2)
2750 return
2751 end if
2752
2753 numCol = col_getNumCol(column)
2754 allocate(hSfc(1, numCol))
2755 do colIndex = 1, numCol
2756 hSfc(1,colIndex) = col_getHeight(column,1,colIndex, 'SF')
2757 end do
2758
2759 call fetch3DLevels_r8(column%vco, hSfc, fldM_opt=hPtrM, fldT_opt=hPtrT)
2760 Z_M(:,:) = hPtrM(1,:,:)
2761 Z_T(:,:) = hPtrT(1,:,:)
2762 deallocate(hPtrM, hPtrT)
2763
2764 deallocate(hSfc)
2765 call msg('calcHeight_col_nl_vcode2100x (czp)', 'END', verb_opt=4)
2766 end subroutine calcHeight_col_nl_vcode2100x
2767
2768 !---------------------------------------------------------
2769 ! calcHeight_col_nl_vcode5xxx
2770 !---------------------------------------------------------
2771 subroutine calcHeight_col_nl_vcode5xxx(column, Z_T, Z_M)
2772 !
2773 !:Purpose: Compute heights for GEM-P columns, return height values
2774 ! in pointer arguments.
2775 !
2776 implicit none
2777
2778 ! Arguments:
2779 type(struct_columnData), intent(in) :: column ! reference column containing temperature and geopotential
2780 real(8), pointer, intent(inout) :: Z_T(:,:) ! computed column height values on thermodynamic levels
2781 real(8), pointer, intent(inout) :: Z_M(:,:) ! computed column height values on momentum levels
2782
2783 ! Developement notes (@mad001)
2784 ! Null subroutine, no computation needed at time of writing.
2785 ! The code is traversed because of `calcZandP_nl` call in `cvt` (agnostic if
2786 ! dealing with GEM-P or GEM-H), but the results for heights are not used
2787 ! at this time.
2788 ! We keep that stub however for future functionalities.
2789 call msg('calcHeight_col_nl_vcode5xxx (czp)', 'END (nothing done)', verb_opt=4)
2790 return
2791
2792 ! to prevent 'variable not used' remark
2793 if (.false.) then
2794 write(*,*) column%nk
2795 Z_T = 0.0
2796 Z_M = 0.0
2797 end if
2798
2799 end subroutine calcHeight_col_nl_vcode5xxx
2800
2801 !---------------------------------------------------------
2802 ! calcHeight_col_tl
2803 !---------------------------------------------------------
2804 subroutine calcHeight_col_tl(columnInc,columnIncRef)
2805 !
2806 !:Purpose: Tangent height computation on the column.
2807 !
2808 implicit none
2809
2810 ! Arguments:
2811 type(struct_columnData), intent(inout) :: columnInc
2812 type(struct_columnData), intent(in) :: columnIncRef
2813
2814 ! Locals:
2815 integer :: Vcode
2816
2817 call utl_tmg_start(173,'low-level--czp_calcHeight_tl')
2818 call msg('calcHeight_col_tl (czp)', 'START', verb_opt=2)
2819
2820 Vcode = col_getVco(columnIncRef)%vcode
2821 if (Vcode == 5005 .or. Vcode == 5002) then
2822 if ( .not. col_varExist(columnInc,'P_*') ) then
2823 call utl_abort('calcHeight_col_tl (czp): for vcode 5xxx, variables P_M and P_T must be allocated in column')
2824 end if
2825 if ( .not. col_varExist(columnInc,'Z_*') ) then
2826 call utl_abort('calcHeight_col_tl (czp): for vcode 5xxx, variables Z_M and Z_T must be allocated in column')
2827 end if
2828 if ( .not. col_varExist(columnInc,'TT') ) then
2829 call utl_abort('calcHeight_col_tl (czp): for vcode 5xxx, variable TT must be allocated in column')
2830 end if
2831 if ( .not. col_varExist(columnInc,'HU') ) then
2832 call utl_abort('calcHeight_col_tl (czp): for vcode 5xxx, variable HU must be allocated in column')
2833 end if
2834 if ( .not. col_varExist(columnInc,'P0') ) then
2835 call utl_abort('calcHeight_col_tl (czp): for vcode 5xxx, variable P0 must be allocated in column')
2836 end if
2837 call calcHeight_col_tl_vcode5xxx
2838 else if (Vcode == 21001) then
2839 ! Development notes (@mad001)
2840 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
2841 call calcHeight_col_tl_vcode2100x
2842 else
2843 call utl_abort('calcHeight_col_tl (czp): not implemented')
2844 end if
2845
2846 call msg('calcHeight_col_tl (czp)', 'END', verb_opt=2)
2847 call utl_tmg_stop(173)
2848
2849 contains
2850 !---------------------------------------------------------
2851 ! calcHeight_col_tl_vcode2100x
2852 !---------------------------------------------------------
2853 subroutine calcHeight_col_tl_vcode2100x
2854 implicit none
2855
2856 call utl_abort('calcHeight_col_tl: vcode 21001 not implemented yet')
2857
2858 end subroutine calcHeight_col_tl_vcode2100x
2859
2860 !---------------------------------------------------------
2861 ! calcHeight_col_tl_vcode5xxx
2862 !---------------------------------------------------------
2863 subroutine calcHeight_col_tl_vcode5xxx
2864 implicit none
2865
2866 ! Locals:
2867 integer :: lev_M,lev_T,nlev_M,nlev_T,colIndex,numColumns
2868 real(8) :: ScaleFactorBottom, ScaleFactorTop
2869 real(8), allocatable :: delThick(:,:)
2870 real(8), pointer :: height_T_ptr(:,:),height_M_ptr(:,:)
2871 real(8), pointer :: P_T(:,:), P_M(:,:)
2872 real(8), pointer :: delHeight_M_ptr(:,:),delHeight_T_ptr(:,:)
2873 real(8), pointer :: delTT(:,:),delHU(:,:),delP0(:,:)
2874 real(8), pointer :: delP_T(:,:), delP_M(:,:)
2875
2876 call msg('calcHeight_col_tl_vcode5xxx (czp)', 'START', verb_opt=4)
2877
2878 nlev_T = col_getNumLev(columnIncRef,'TH')
2879 nlev_M = col_getNumLev(columnIncRef,'MM')
2880
2881 numColumns = col_getNumCol(columnInc)
2882
2883 allocate(delThick(nlev_T,numColumns))
2884 delThick(:,:) = 0.0d0
2885
2886 ! generate the height coefficients on the grid
2887 call calcHeightCoeff_col(columnIncRef)
2888
2889 ! loop over all lat/lon/step
2890
2891 height_M_ptr => col_getAllColumns(columnIncRef,'Z_M')
2892 height_T_ptr => col_getAllColumns(columnIncRef,'Z_T')
2893 P_M => col_getAllColumns(columnIncRef,'P_M')
2894 P_T => col_getAllColumns(columnIncRef,'P_T')
2895
2896 delHeight_M_ptr => col_getAllColumns(columnInc,'Z_M')
2897 delHeight_T_ptr => col_getAllColumns(columnInc,'Z_T')
2898 delTT => col_getAllColumns(columnInc,'TT')
2899 delHU => col_getAllColumns(columnInc,'HU')
2900 delP0 => col_getAllColumns(columnInc,'P0')
2901 delP_M => col_getAllColumns(columnInc,'P_M')
2902 delP_T => col_getAllColumns(columnInc,'P_T')
2903
2904 ! ensure increment at sfc is zero (fixed height level)
2905 delHeight_M_ptr(nlev_M,:) = 0.0d0
2906 delHeight_T_ptr(nlev_T,:) = 0.0d0
2907
2908 if_computeHeight_col_tl_vcodes : if (Vcode == 5002) then
2909
2910 ! compute increment to thickness for each layer between the two momentum levels
2911 do colIndex = 1, numColumns
2912 do lev_T = 2, (nlev_T-1)
2913 delThick(lev_T,colIndex) = &
2914 coeff_M_TT_col(lev_T,colIndex) * delTT(lev_T,colIndex) + &
2915 coeff_M_HU_col(lev_T,colIndex) * delHU(lev_T,colIndex) + &
2916 coeff_M_P0_delPM_col(lev_T,colIndex) * &
2917 ( delP_M(lev_T ,colIndex) / P_M(lev_T ,colIndex) - &
2918 delP_M(lev_T-1,colIndex) / P_M(lev_T-1,colIndex) ) + &
2919 coeff_M_P0_dP_delPT_col(lev_T,colIndex) * &
2920 delP_T(lev_T,colIndex) + &
2921 coeff_M_P0_dP_delP0_col(lev_T,colIndex) * delP0(1,colIndex)
2922 end do
2923 end do
2924
2925 ! compute height increment on momentum levels above the surface
2926 do colIndex = 1, numColumns
2927 do lev_M = (nlev_M-1), 1, -1
2928 lev_T = lev_M + 1 ! thermo level just below momentum level being computed
2929 delHeight_M_ptr(lev_M,colIndex) = &
2930 delHeight_M_ptr(lev_M+1,colIndex) + delThick(lev_T,colIndex)
2931 end do
2932 end do
2933
2934 ! compute height increment on thermo levels using weighted average of height increment of momentum levels
2935 do colIndex = 1, numColumns
2936 do lev_T = 1, (nlev_T-1)
2937 if ( lev_T == 1) then
2938 ! compute height increment for top thermo level (from top momentum level)
2939 delHeight_T_ptr(1,colIndex) = delHeight_M_ptr(1,colIndex) + &
2940 coeff_T_TT_col(colIndex) * delTT(1,colIndex) + &
2941 coeff_T_HU_col(colIndex) * delHU(1,colIndex) + &
2942 coeff_T_P0_delP1_col(colIndex) * &
2943 ( delP_M(1,colIndex) / P_M(1,colIndex) - &
2944 delP_T(1,colIndex) / P_T(1,colIndex) ) + &
2945 coeff_T_P0_dP_delPT_col(colIndex) * delP_T(1,colIndex) + &
2946 coeff_T_P0_dP_delP0_col(colIndex) * delP0(1,colIndex)
2947 else
2948 lev_M = lev_T ! momentum level just below thermo level being computed
2949 ScaleFactorBottom = (height_T_ptr(lev_T,colIndex) - &
2950 height_M_ptr(lev_M-1,colIndex)) / &
2951 (height_M_ptr(lev_M,colIndex) - height_M_ptr(lev_M-1,colIndex))
2952 ScaleFactorTop = 1 - ScaleFactorBottom
2953 delHeight_T_ptr(lev_T,colIndex) = &
2954 ScaleFactorBottom * delHeight_M_ptr(lev_M ,colIndex) + &
2955 ScaleFactorTop * delHeight_M_ptr(lev_M-1,colIndex)
2956 end if
2957 end do
2958 end do
2959
2960 else if(Vcode == 5005) then if_computeHeight_col_tl_vcodes
2961
2962 ! compute increment to thickness for each layer between the two momentum levels
2963 do colIndex = 1, numColumns
2964 do lev_T = 1, (nlev_T-1)
2965 delThick(lev_T,colIndex) = &
2966 coeff_M_TT_col(lev_T,colIndex) * delTT(lev_T,colIndex) + &
2967 coeff_M_HU_col(lev_T,colIndex) * delHU(lev_T,colIndex) + &
2968 coeff_M_P0_delPM_col(lev_T,colIndex) * &
2969 ( delP_M(lev_T+1,colIndex) / P_M(lev_T+1,colIndex) - &
2970 delP_M(lev_T ,colIndex) / P_M(lev_T ,colIndex) ) + &
2971 coeff_M_P0_dP_delPT_col(lev_T,colIndex) * &
2972 delP_T(lev_T,colIndex) + &
2973 coeff_M_P0_dP_delP0_col(lev_T,colIndex) * delP0(1,colIndex)
2974 end do
2975 end do
2976
2977 ! compute height increment on momentum levels above the surface
2978 do colIndex = 1, numColumns
2979 do lev_M = (nlev_M-1), 1, -1
2980 lev_T = lev_M ! thermo level just below momentum level being computed
2981 delHeight_M_ptr(lev_M,colIndex) = &
2982 delHeight_M_ptr(lev_M+1,colIndex) + delThick(lev_T,colIndex)
2983 end do
2984 end do
2985
2986 ! compute height increment on thermo levels using weighted average of height increment of momentum levels
2987 do colIndex = 1, numColumns
2988 do lev_T = 1, (nlev_T-1)
2989 lev_M = lev_T + 1 ! momentum level just below thermo level being computed
2990 ScaleFactorBottom = &
2991 (height_T_ptr(lev_T,colIndex) - height_M_ptr(lev_M-1,colIndex)) / &
2992 (height_M_ptr(lev_M,colIndex) - height_M_ptr(lev_M-1,colIndex))
2993 ScaleFactorTop = 1 - ScaleFactorBottom
2994 delHeight_T_ptr(lev_T,colIndex) = &
2995 ScaleFactorBottom * delHeight_M_ptr(lev_M ,colIndex) + &
2996 ScaleFactorTop * delHeight_M_ptr(lev_M-1,colIndex)
2997 end do
2998 end do
2999 else
3000 call utl_abort('calcHeight_col_tl_vcode5xxx (czp): not implemented')
3001 end if if_computeHeight_col_tl_vcodes
3002
3003 deallocate(delThick)
3004
3005 call msg('calcHeight_col_tl_vcode5xxx (czp)', 'END', verb_opt=4)
3006 end subroutine calcHeight_col_tl_vcode5xxx
3007
3008 end subroutine calcHeight_col_tl
3009
3010 !---------------------------------------------------------
3011 ! calcHeight_col_ad
3012 !---------------------------------------------------------
3013 subroutine calcHeight_col_ad(columnInc,columnIncRef)
3014 !
3015 !:Purpose: Adjoint of height computation on the column.
3016 !
3017 !
3018 implicit none
3019
3020 ! Arguments:
3021 type(struct_columnData), intent(inout) :: columnInc
3022 type(struct_columnData), intent(in) :: columnIncRef
3023
3024 ! Locals:
3025 integer :: Vcode
3026
3027 call utl_tmg_start(174,'low-level--czp_calcHeight_ad')
3028 call msg('calcHeight_col_ad (czp)', 'START', verb_opt=2)
3029
3030 Vcode = col_getVco(columnIncRef)%vcode
3031 if (Vcode == 5005 .or. Vcode == 5002) then
3032 if ( .not. col_varExist(columnInc,'P_*') ) then
3033 call utl_abort('calcHeight_col_ad (czp): for vcode 5xxx, variables P_M and P_T must be allocated in column')
3034 end if
3035 if ( .not. col_varExist(columnInc,'Z_*') ) then
3036 call utl_abort('calcHeight_col_ad (czp): for vcode 5xxx, variables Z_M and Z_T must be allocated in column')
3037 end if
3038 if ( .not. col_varExist(columnInc,'TT') ) then
3039 call utl_abort('calcHeight_col_ad (czp): for vcode 5xxx, variable TT must be allocated in column')
3040 end if
3041 if ( .not. col_varExist(columnInc,'HU') ) then
3042 call utl_abort('calcHeight_col_ad (czp): for vcode 5xxx, variable HU must be allocated in column')
3043 end if
3044 if ( .not. col_varExist(columnInc,'P0') ) then
3045 call utl_abort('calcHeight_col_ad (czp): for vcode 5xxx, variable P0 must be allocated in column')
3046 end if
3047 call calcHeight_col_ad_vcode5xxx
3048 else if (Vcode == 21001) then
3049 ! Development notes (@mad001)
3050 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
3051 call calcHeight_col_ad_vcode2100x
3052 else
3053 call utl_abort('calcHeight_col_ad (czp): not implemented')
3054 end if
3055
3056 call msg('calcHeight_col_ad (czp)', 'END', verb_opt=2)
3057 call utl_tmg_stop(174)
3058
3059 contains
3060 !---------------------------------------------------------
3061 ! calcHeight_col_ad_vcode2100x
3062 !---------------------------------------------------------
3063 subroutine calcHeight_col_ad_vcode2100x
3064 implicit none
3065
3066 call utl_abort('calcHeight_col_ad (czp): vcode 21001 not implemented yet')
3067
3068 end subroutine calcHeight_col_ad_vcode2100x
3069
3070 !---------------------------------------------------------
3071 ! calcHeight_col_ad_vcode5xxx
3072 !---------------------------------------------------------
3073 subroutine calcHeight_col_ad_vcode5xxx
3074 implicit none
3075
3076 ! Locals:
3077 integer :: lev_M,lev_T,nlev_M,nlev_T,numColumns,colIndex
3078 real(8) :: ScaleFactorBottom, ScaleFactorTop
3079 real(8), allocatable :: delThick(:,:)
3080 real(8), pointer :: height_M_ptr(:,:),height_T_ptr(:,:)
3081 real(8), allocatable :: delHeight_M(:,:)
3082 real(8), pointer :: P_M(:,:),P_T(:,:)
3083 real(8), pointer :: delHeight_M_ptr(:,:),delHeight_T_ptr(:,:)
3084 real(8), pointer :: delTT(:,:),delHU(:,:),delP0(:,:)
3085 real(8), pointer :: delP_M(:,:),delP_T(:,:)
3086
3087 call msg('calcHeight_col_ad_vcode5xxx (czp)', 'START', verb_opt=4)
3088
3089 nlev_T = col_getNumLev(columnIncRef,'TH')
3090 nlev_M = col_getNumLev(columnIncRef,'MM')
3091 numColumns = col_getNumCol(columnIncRef)
3092
3093 allocate(delHeight_M(nlev_M,numColumns))
3094 allocate(delThick(0:nlev_T,numColumns))
3095
3096 ! generate the height coefficients on the grid
3097 call calcHeightCoeff_col(columnIncRef)
3098
3099 height_M_ptr => col_getAllColumns(columnIncRef,'Z_M')
3100 height_T_ptr => col_getAllColumns(columnIncRef,'Z_T')
3101 P_M => col_getAllColumns(columnIncRef,'P_M')
3102 P_T => col_getAllColumns(columnIncRef,'P_T')
3103
3104 delHeight_M_ptr => col_getAllColumns(columnInc,'Z_M')
3105 delHeight_T_ptr => col_getAllColumns(columnInc,'Z_T')
3106 delTT => col_getAllColumns(columnInc,'TT')
3107 delHU => col_getAllColumns(columnInc,'HU')
3108 delP0 => col_getAllColumns(columnInc,'P0')
3109 delP_M => col_getAllColumns(columnInc,'P_M')
3110 delP_T => col_getAllColumns(columnInc,'P_T')
3111
3112 delHeight_M(:,:) = delHeight_M_ptr(:,:)
3113
3114 if_computeHeight_col_ad_vcodes : if(Vcode == 5002) then
3115
3116 ! adjoint of compute height increment on thermo levels by simple averaging
3117 do colIndex = 1, numColumns
3118 do lev_T = 1, (nlev_T-1)
3119 lev_M = lev_T ! momentum level just below thermo level being computed
3120
3121 ! adjoint of compute height increment on top thermo level (from top momentum level)
3122 if (lev_T == 1) then
3123 delHeight_M(1,colIndex) = &
3124 delHeight_M(1,colIndex) + &
3125 delHeight_T_ptr(1,colIndex)
3126
3127 delTT(1,colIndex) = &
3128 delTT(1,colIndex) + &
3129 coeff_T_TT_col(colIndex) * delHeight_T_ptr(1,colIndex)
3130
3131 delHU(1,colIndex) = &
3132 delHU(1,colIndex) + &
3133 coeff_T_HU_col (colIndex) * delHeight_T_ptr(1,colIndex)
3134
3135 delP_M(1,colIndex) = &
3136 delP_M(1,colIndex) + &
3137 coeff_T_P0_delP1_col(colIndex) / P_M(1,colIndex) * &
3138 delHeight_T_ptr(1,colIndex)
3139
3140 delP_T(1,colIndex) = &
3141 delP_T(1,colIndex) - &
3142 coeff_T_P0_delP1_col(colIndex) / P_T(1,colIndex) * &
3143 delHeight_T_ptr(1,colIndex)
3144
3145 delP_T(1,colIndex) = &
3146 delP_T(1,colIndex) + &
3147 coeff_T_P0_dP_delPT_col(colIndex) * delHeight_T_ptr(1,colIndex)
3148
3149 delP0(1,colIndex) = &
3150 delP0(1,colIndex) + &
3151 coeff_T_P0_dp_delP0_col(colIndex) * delHeight_T_ptr(1,colIndex)
3152 else
3153 ScaleFactorBottom = &
3154 (height_T_ptr(lev_T,colIndex) - &
3155 height_M_ptr(lev_M-1,colIndex)) / &
3156 (height_M_ptr(lev_M,colIndex) - &
3157 height_M_ptr(lev_M-1,colIndex))
3158 ScaleFactorTop = 1 - ScaleFactorBottom
3159
3160 delHeight_M(lev_M-1,colIndex) = &
3161 delHeight_M(lev_M-1,colIndex) + &
3162 ScaleFactorTop * delHeight_T_ptr(lev_T,colIndex)
3163
3164 delHeight_M(lev_M,colIndex) = &
3165 delHeight_M(lev_M ,colIndex) + &
3166 ScaleFactorBottom * delHeight_T_ptr(lev_T,colIndex)
3167 end if
3168 end do
3169 end do
3170
3171 ! adjoint of compute height increment on momentum levels above the surface
3172 delThick(0:1,:) = 0.0d0
3173 do colIndex = 1, numColumns
3174 do lev_M = 1, (nlev_M-1)
3175 lev_T = lev_M + 1 ! thermo level just below momentum level being computed
3176 delThick(lev_T,colIndex) = &
3177 delThick(lev_T-1,colIndex) + &
3178 delHeight_M(lev_M ,colIndex)
3179 end do
3180 end do
3181
3182 ! adjoint of compute increment to thickness for each layer between the two momentum levels
3183 do colIndex = 1, numColumns
3184 do lev_T = 2, nlev_T-1
3185 delTT(lev_T,colIndex) = &
3186 delTT (lev_T,colIndex) + &
3187 coeff_M_TT_col (lev_T,colIndex) * delThick(lev_T,colIndex)
3188
3189 delHU(lev_T,colIndex) = &
3190 delHU (lev_T,colIndex) + &
3191 coeff_M_HU_col (lev_T,colIndex) * delThick(lev_T,colIndex)
3192
3193 delP_M(lev_T,colIndex)= &
3194 delP_M(lev_T,colIndex) + &
3195 coeff_M_P0_delPM_col(lev_T,colIndex) / P_M(lev_T,colIndex) * &
3196 delThick(lev_T,colIndex)
3197
3198 delP_M(lev_T-1,colIndex) = &
3199 delP_M(lev_T-1,colIndex) - &
3200 coeff_M_P0_delPM_col(lev_T,colIndex) / P_M(lev_T-1,colIndex)*&
3201 delThick(lev_T,colIndex)
3202
3203 delP_T(lev_T,colIndex) = &
3204 delP_T(lev_T,colIndex) + &
3205 coeff_M_P0_dP_delPT_col(lev_T,colIndex) * &
3206 delThick(lev_T,colIndex)
3207
3208 delP0(1,colIndex) = &
3209 delP0(1,colIndex) + &
3210 coeff_M_P0_dP_delP0_col(lev_T,colIndex) * &
3211 delThick(lev_T,colIndex)
3212 end do
3213 end do
3214
3215 else if(Vcode == 5005) then if_computeHeight_col_ad_vcodes
3216
3217 ! adjoint of compute height increment on thermo levels by simple averaging
3218 do colIndex = 1, numColumns
3219 do lev_T = 1, (nlev_T-1)
3220 lev_M = lev_T+1 ! momentum level just below thermo level being computed
3221 ScaleFactorBottom = (height_T_ptr(lev_T,colIndex) - &
3222 height_M_ptr(lev_M-1,colIndex)) / &
3223 (height_M_ptr(lev_M,colIndex) - &
3224 height_M_ptr(lev_M-1,colIndex))
3225 ScaleFactorTop = 1 - ScaleFactorBottom
3226 delHeight_M(lev_M-1,colIndex) = delHeight_M(lev_M-1,colIndex) + &
3227 ScaleFactorTop * delHeight_T_ptr(lev_T ,colIndex)
3228 delHeight_M(lev_M,colIndex) = delHeight_M(lev_M ,colIndex) + &
3229 ScaleFactorBottom * delHeight_T_ptr(lev_T ,colIndex)
3230 end do
3231 end do
3232
3233 ! adjoint of compute height increment on momentum levels
3234 delThick(0,:) = 0.0d0
3235 do colIndex = 1, numColumns
3236 do lev_M = 1, (nlev_M-1)
3237 lev_T = lev_M ! thermo level just below momentum level being computed
3238 delThick(lev_T,colIndex) = delThick(lev_T-1,colIndex) + &
3239 delHeight_M (lev_M ,colIndex)
3240 end do
3241 end do
3242
3243 do colIndex = 1, numColumns
3244 do lev_T = 1, nlev_T-1
3245 delTT(lev_T,colIndex) = &
3246 delTT(lev_T,colIndex) + &
3247 coeff_M_TT_col(lev_T,colIndex) * delThick(lev_T,colIndex)
3248
3249 delHU(lev_T,colIndex) = &
3250 delHU(lev_T,colIndex) + &
3251 coeff_M_HU_col(lev_T,colIndex) * delThick(lev_T,colIndex)
3252
3253 delP_M(lev_T+1,colIndex) = &
3254 delP_M(lev_T+1,colIndex) + &
3255 coeff_M_P0_delPM_col(lev_T,colIndex) / P_M(lev_T+1,colIndex)*&
3256 delThick(lev_T,colIndex)
3257
3258 delP_M(lev_T ,colIndex) = &
3259 delP_M(lev_T ,colIndex) - &
3260 coeff_M_P0_delPM_col(lev_T,colIndex) / P_M(lev_T ,colIndex)*&
3261 delThick(lev_T,colIndex)
3262
3263 delP_T(lev_T ,colIndex) = &
3264 delP_T(lev_T ,colIndex) + &
3265 coeff_M_P0_dP_delPT_col(lev_T,colIndex) * &
3266 delThick(lev_T,colIndex)
3267
3268 delP0(1,colIndex) = &
3269 delP0(1,colIndex) + &
3270 coeff_M_P0_dP_delP0_col(lev_T,colIndex) * &
3271 delThick(lev_T,colIndex)
3272 end do
3273 end do
3274 else
3275 call utl_abort('calcHeight_col_ad_vcode5xxx (czp): not implemented')
3276 end if if_computeHeight_col_ad_vcodes
3277
3278 deallocate(delThick)
3279 deallocate(delHeight_M)
3280
3281 call msg('calcHeight_col_ad_vcode5xxx (czp)', 'END', verb_opt=4)
3282 end subroutine calcHeight_col_ad_vcode5xxx
3283
3284 end subroutine calcHeight_col_ad
3285
3286 !---------------------------------------------------------
3287 ! calcPressure_col_nl
3288 !---------------------------------------------------------
3289 subroutine calcPressure_col_nl(column)
3290 !
3291 !:Purpose: Pressure computation on the column.
3292 !
3293 implicit none
3294
3295 ! Arguments:
3296 type(struct_columnData), intent(inout) :: column
3297
3298 ! Locals:
3299 real(8), pointer :: P_T(:,:), P_M(:,:)
3300
3301 call msg('calcPressure_col_nl (czp)', 'START', verb_opt=2)
3302
3303 P_T => col_getAllColumns(column, 'P_T')
3304 P_M => col_getAllColumns(column, 'P_M')
3305 call czp_calcReturnPressure_col_nl(column, P_T, P_M)
3306
3307 call msg('calcPressure_col_nl (czp)', &
3308 new_line('')//'P_M = '//str(col_getColumn(column,1,'P_M')) &
3309 //new_line('')//'P_T = '//str(col_getColumn(column,1,'P_T')), &
3310 verb_opt=2)
3311
3312 call msg('calcPressure_col_nl (czp)', 'END', verb_opt=2)
3313 end subroutine calcPressure_col_nl
3314
3315 !---------------------------------------------------------
3316 ! czp_calcReturnPressure_col_nl
3317 !---------------------------------------------------------
3318 subroutine czp_calcReturnPressure_col_nl(column, P_T, P_M)
3319 !
3320 !:Purpose: Pressure computation on the column, return values in pointers.
3321 !
3322 implicit none
3323
3324 ! Arguments:
3325 type(struct_columnData), intent(in) :: column ! reference column
3326 real(8), pointer, intent(inout) :: P_T(:,:) ! computed column pressure values on thermodynamic levels
3327 real(8), pointer, intent(inout) :: P_M(:,:) ! computed column pressure values on momentum levels
3328
3329 ! Locals:
3330 integer :: Vcode
3331
3332 call msg('czp_calcReturnPressure_col_nl (czp)', 'START', verb_opt=2)
3333
3334 Vcode = col_getVco(column)%vcode
3335 if (Vcode == 5005 .or. Vcode == 5002 .or. Vcode == 5100) then
3336 if ( .not. col_varExist(column,'P0') ) then
3337 call utl_abort('czp_calcReturnPressure_col_nl (czp): for vcode 5xxx, variable P0 must be allocated in column')
3338 end if
3339 call calcPressure_col_nl_vcode5xxx(column, P_T, P_M)
3340 else if (Vcode == 21001) then
3341 if ( .not. (col_varExist(column,'P0') .and. col_varExist(column,'TT') .and. &
3342 col_varExist(column,'HU')) ) then
3343 call utl_abort('czp_calcReturnPressure_col_nl (czp): for vcode 2100x, variables P0, TT and HU must be allocated in column')
3344 end if
3345 call calcPressure_col_nl_vcode2100x(column, P_T, P_M)
3346 end if
3347
3348 call msg('czp_calcReturnPressure_col_nl (czp)', 'END', verb_opt=2)
3349 end subroutine czp_calcReturnPressure_col_nl
3350
3351 !---------------------------------------------------------
3352 ! calcPressure_col_nl_vcode2100x
3353 !---------------------------------------------------------
3354 subroutine calcPressure_col_nl_vcode2100x(column, P_T, P_M)
3355 !
3356 !:Purpose: Compute pressure and return values in pointer arguments.
3357 ! GEM-H column input.
3358 !
3359 implicit none
3360 ! Development notes (@mad001)
3361 ! if feasible, consider reusing the same code for both
3362 ! `calcPressure_{gsv,col}_nl_vcode2100x`
3363 ! (@mab001) Also should remove need for `Z_T/M` to be allocated in column
3364 ! and use local array instead.
3365
3366 ! Arguments:
3367 type(struct_columnData), intent(in) :: column ! reference column
3368 real(8), pointer, intent(inout) :: P_T(:,:) ! computed column pressure values on thermodynamic levels
3369 real(8), pointer, intent(inout) :: P_M(:,:) ! computed column pressure values on momentum levels
3370
3371 ! Locals:
3372 real(8), allocatable :: tv(:)
3373 integer :: numCol, nLev_T, nLev_M
3374 integer :: colIndex, lev_T, lev_M
3375 real(8) :: lat, sLat, cLat
3376 real(8) :: P0, rMT, hu, tt, tv0, cmp, dh, Rgh
3377 real(8) :: scaleFactorBottom, logP
3378 real(8) :: Z_T, Z_M, Z_M1
3379
3380 call msg('calcPressure_col_nl_vcode2100x (czp)', 'START', verb_opt=4)
3381
3382 numCol = col_getNumCol(column)
3383 nLev_M = col_getNumLev(column, 'MM')
3384 nLev_T = col_getNumLev(column, 'TH')
3385
3386 allocate(tv(nLev_T))
3387
3388 do_onAllcolumns: do colIndex = lbound(column%lat,1), ubound(column%lat,1)
3389 ! column%lat populated in innovation_mod from obsSpaceData latitudes
3390 lat = column%lat(colIndex)
3391 sLat = sin(lat)
3392 cLat = cos(lat)
3393
3394 ! surface values
3395 P0 = col_getElem( column, 1, colIndex, 'P0') ! surface pressure
3396 rMT = col_getHeight(column, 1, colIndex, 'SF') ! surface height
3397
3398 ! compute pressure on diagnostic (nLev_{T,M}) levels
3399 hu = col_getElem(column, nLev_T, colIndex, 'HU')
3400 tt = col_getElem(column, nLev_T, colIndex, 'TT')
3401 tv0 = phf_fotvt8(tt,hu)
3402
3403 ! thermo diagnostic level
3404 Z_T = col_getHeight(column, nLev_T, colIndex, 'TH')
3405 cmp = gpscompressibility(P0,tt,hu)
3406 tv(nlev_T) = tv0*cmp
3407 dh = Z_T - rMT
3408 Rgh = phf_gravityalt(sLat, rMT+0.5D0*dh)
3409 P_T(colIndex, nlev_T) = P0*exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(nlev_T))
3410
3411 ! momentum diagnostic level
3412 Z_M = col_getHeight(column,nLev_M,colIndex,'MM')
3413 dh = Z_M - rMT
3414 Rgh = phf_gravityalt(sLat, rMT+0.5D0*dh)
3415 P_M(colIndex, nlev_M) = P0*exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(nlev_T))
3416
3417 call msg('calcPressure_col_nl_vcode2100x (czp)', &
3418 'Column index '//str(colIndex) //': lat='//str(lat)&
3419 //' Surface: height='//str(rMT)//', P0='//str(P0) &
3420 //' TH_diag_lvl ('//str(nLev_T)//'): height='//str(Z_T)&
3421 //', P_T='//str(P_T(colIndex, nlev_T)) &
3422 //' MM_diag_lvl ('//str(nLev_M)//'): height='//str(Z_M)&
3423 //', P_M='//str(P_M(colIndex, nlev_M)), &
3424 verb_opt=6)
3425
3426 ! compute pressure on all levels above except the last
3427 do lev_M = nlev_M-1, 1, -1
3428 lev_T = lev_M ! thermo level just below
3429 hu = col_getElem( column, lev_T, colIndex, 'HU')
3430 tt = col_getElem( column, lev_T, colIndex, 'TT')
3431 Z_M = col_getHeight(column, lev_M, colIndex, 'MM')
3432 Z_M1 = col_getHeight(column, lev_M+1, colIndex, 'MM')
3433 Z_T = col_getHeight(column, lev_T, colIndex, 'TH')
3434
3435 tv0 = phf_fotvt8(tt,hu)
3436 dh = Z_M - Z_M1
3437 Rgh = phf_gravityalt(sLat, Z_M1+0.5D0*dh)
3438
3439 ! approximation of tv from pressure on previous momentum level
3440 cmp = gpscompressibility(P_M(colIndex, lev_M+1),tt,hu)
3441 tv(lev_T) = tv0*cmp
3442 P_M(colIndex, lev_M) = P_M(colIndex, lev_M+1) * &
3443 exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(lev_T))
3444 ! first interpolation of thermo pressure
3445 scaleFactorBottom = (Z_T-Z_M1)/(Z_M-Z_M1)
3446 logP = (1.0D0-scaleFactorBottom)*log(P_M(colIndex, lev_M+1)) + &
3447 scaleFactorBottom*log(P_M(colIndex, lev_M))
3448 P_T(colIndex, lev_T) = exp(logP)
3449
3450 ! second iteration on tv
3451 cmp = gpscompressibility(P_T(colIndex, lev_T),tt,hu)
3452 tv(lev_T) = tv0*cmp
3453 P_M(colIndex, lev_M) = P_M(colIndex, lev_M+1) * &
3454 exp(-Rgh*dh/MPC_RGAS_DRY_AIR_R8/tv(lev_T))
3455
3456 ! second iteration interpolation of thermo pressure
3457 logP = (1.0D0-scaleFactorBottom)*log(P_M(colIndex, lev_M+1)) + &
3458 scaleFactorBottom*log(P_M(colIndex, lev_M))
3459 P_T(colIndex, lev_T) = exp(logP)
3460
3461 end do
3462 end do do_onAllColumns
3463
3464 deallocate(tv)
3465
3466 call msg('calcPressure_col_nl_vcode2100x (czp)', 'END', verb_opt=4)
3467 end subroutine calcPressure_col_nl_vcode2100x
3468
3469 !---------------------------------------------------------
3470 ! calcPressure_col_nl_vcode5xxx
3471 !---------------------------------------------------------
3472 subroutine calcPressure_col_nl_vcode5xxx(column, P_T, P_M)
3473 implicit none
3474
3475 ! Arguments:
3476 type(struct_columnData), intent(in) :: column ! reference column
3477 real(8), pointer, intent(inout) :: P_T(:,:) ! computed column pressure values on thermodynamic levels
3478 real(8), pointer, intent(inout) :: P_M(:,:) ! computed column pressure values on momentum levels
3479
3480 ! Locals:
3481 real(kind=8), allocatable :: Psfc(:,:), PsfcLS(:,:)
3482 real(kind=8), pointer :: zppobsM(:,:,:), zppobsT(:,:,:)
3483 integer :: headerIndex, Vcode
3484
3485 call msg('calcPressure_col_nl_vcode5xxx (czp)', 'START', verb_opt=4)
3486 if ( col_getNumCol(column) <= 0 ) then
3487 call msg('calcPressure_col_nl_vcode5xxx (czp)',&
3488 'END (number of columns <= 0)', verb_opt=2)
3489 return
3490 end if
3491
3492 Vcode = col_getVco(column)%vcode
3493
3494 if (.not.col_varExist(column,'P0')) then
3495 call utl_abort('calcPressure_col_nl (czp): P0 must be present as an analysis variable!')
3496 end if
3497
3498 allocate(Psfc(1,col_getNumCol(column)))
3499 do headerIndex = 1,col_getNumCol(column)
3500 Psfc(1,headerIndex) = col_getElem(column,1,headerIndex,'P0')
3501 end do
3502
3503 if (Vcode == 5100) then
3504 if (.not.col_varExist(column,'P0LS')) then
3505 call utl_abort('calcPressure_col_nl (czp): P0LS must be present as an analysis variable!')
3506 end if
3507
3508 allocate(PsfcLS(1,col_getNumCol(column)))
3509 do headerIndex = 1,col_getNumCol(column)
3510 PsfcLS(1,headerIndex) = col_getElem(column,1,headerIndex,'P0LS')
3511 end do
3512
3513 call fetch3DLevels_r8(column%vco, Psfc ,sfcFldLS_opt=PsfcLS, &
3514 fldM_opt=zppobsM, fldT_opt=zppobsT)
3515 deallocate(PsfcLS)
3516 else
3517 call fetch3DLevels_r8(column%vco, Psfc ,fldM_opt=zppobsM, fldT_opt=zppobsT)
3518 end if
3519 P_M(:,:) = zppobsM(1,:,:)
3520 P_T(:,:) = zppobsT(1,:,:)
3521 deallocate(zppobsM, zppobsT)
3522 deallocate(Psfc)
3523
3524 call msg('calcPressure_col_nl_vcode5xxx (czp)', 'END', verb_opt=4)
3525 end subroutine calcPressure_col_nl_vcode5xxx
3526
3527
3528 !---------------------------------------------------------
3529 ! calcPressure_col_tl
3530 !---------------------------------------------------------
3531 subroutine calcPressure_col_tl(columnInc, columnIncRef)
3532 !
3533 !:Purpose: Tangent pressure computation on the column.
3534 !
3535 implicit none
3536
3537 ! Arguments:
3538 type(struct_columnData), intent(inout) :: columnInc ! column that will contain the P_T/P_M increments
3539 type(struct_columnData), intent(in) :: columnIncRef ! column containing needed reference fields
3540
3541 ! Locals:
3542 integer :: Vcode
3543
3544 call msg('calcPressure_col_tl (czp)', 'START', verb_opt=2)
3545
3546 Vcode = col_getVco(columnIncRef)%vcode
3547 if (Vcode == 5005 .or. Vcode == 5002) then
3548 if ( .not. col_varExist(columnInc,'P_*') ) then
3549 call utl_abort('calcPressure_col_tl (czp): for vcode 5xxx, variables P_M and P_T must be allocated in column')
3550 end if
3551 if ( .not. col_varExist(columnInc,'P0') ) then
3552 call utl_abort('calcPressure_col_tl (czp): for vcode 5xxx, variable P0 must be allocated in column')
3553 end if
3554 call calcPressure_col_tl_vcode5xxx
3555 else if (Vcode == 21001) then
3556 ! Development notes (@mad001)
3557 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
3558 call calcPressure_col_tl_vcode2100x
3559 else
3560 call utl_abort('calcPressure_col_tl (czp): not implemented')
3561 end if
3562
3563 call msg('calcPressure_col_tl (czp)', 'END', verb_opt=2)
3564
3565 contains
3566 !---------------------------------------------------------
3567 ! calcPressure_col_tl_vcode2100x
3568 !---------------------------------------------------------
3569 subroutine calcPressure_col_tl_vcode2100x
3570 implicit none
3571
3572 call utl_abort('calcPressure_col_tl (czp): vcode 21001 not implemented yet')
3573
3574 end subroutine calcPressure_col_tl_vcode2100x
3575
3576 !---------------------------------------------------------
3577 ! calcPressure_col_tl_vcode5xxx
3578 !---------------------------------------------------------
3579 subroutine calcPressure_col_tl_vcode5xxx
3580 implicit none
3581
3582 ! Locals:
3583 real(8) :: Psfc
3584 real(8), pointer :: delPsfc(:,:), PsfcRef(:,:)
3585 real(8), pointer :: delP_T(:,:), delP_M(:,:)
3586 real(8), pointer :: dP_dPsfc_T(:), dP_dPsfc_M(:)
3587 integer :: status, colIndex
3588 integer :: lev_M, lev_T, nlev_T, nlev_M, numColumns
3589
3590 call msg('calcPressure_col_tl_vcode5xxx (czp)', 'START', verb_opt=4)
3591
3592 nullify(dP_dPsfc_T)
3593 nullify(dP_dPsfc_M)
3594 nullify(delPsfc)
3595 nullify(delP_T)
3596 nullify(delP_M)
3597
3598 delP_M => col_getAllColumns(columnInc,'P_M')
3599 delP_T => col_getAllColumns(columnInc,'P_T')
3600 delPsfc => col_getAllColumns(columnInc,'P0')
3601 PsfcRef => col_getAllColumns(columnIncRef,'P0')
3602
3603 nlev_T = col_getNumLev(columnInc,'TH')
3604 nlev_M = col_getNumLev(columnInc,'MM')
3605 numColumns = col_getNumCol(columnInc)
3606
3607 do colIndex = 1, numColumns
3608
3609 Psfc = PsfcRef(1,colIndex)
3610
3611 ! dP_dPsfc_M
3612 nullify(dP_dPsfc_M)
3613 status = vgd_dpidpis(columnInc%vco%vgrid, &
3614 columnInc%vco%ip1_M, &
3615 dP_dPsfc_M, &
3616 Psfc)
3617 if( status .ne. VGD_OK ) then
3618 call utl_abort('calcPressure_col_tl (czp): ERROR with vgd_dpidpis')
3619 end if
3620 ! calculate delP_M
3621 do lev_M = 1, nlev_M
3622 delP_M(lev_M,colIndex) = dP_dPsfc_M(lev_M) * delPsfc(1,colIndex)
3623 end do
3624 deallocate(dP_dPsfc_M)
3625
3626 ! dP_dPsfc_T
3627 nullify(dP_dPsfc_T)
3628 status = vgd_dpidpis(columnInc%vco%vgrid, &
3629 columnInc%vco%ip1_T, &
3630 dP_dPsfc_T, &
3631 Psfc)
3632 if( status .ne. VGD_OK ) then
3633 call utl_abort('calcPressure_col_tl (czp): ERROR with vgd_dpidpis')
3634 end if
3635 ! calculate delP_T
3636 do lev_T = 1, nlev_T
3637 delP_T(lev_T,colIndex) = dP_dPsfc_T(lev_T) * delPsfc(1,colIndex)
3638 end do
3639 deallocate(dP_dPsfc_T)
3640
3641 end do
3642
3643 call msg('calcPressure_col_tl_vcode5xxx (czp)', 'END', verb_opt=4)
3644 end subroutine calcPressure_col_tl_vcode5xxx
3645
3646 end subroutine calcPressure_col_tl
3647
3648 !---------------------------------------------------------
3649 ! calcPressure_col_ad
3650 !---------------------------------------------------------
3651 subroutine calcPressure_col_ad( columnInc, columnIncRef)
3652 !
3653 !:Purpose: Adjoint of pressure computation on the column.
3654 !
3655 implicit none
3656
3657 ! Arguments:
3658 type(struct_columnData), intent(inout) :: columnInc ! column that will contain increments of P_M/P_T
3659 type(struct_columnData), intent(in) :: columnIncRef ! column containing needed reference fields
3660
3661 ! Locals:
3662 integer :: Vcode
3663
3664 call msg('calcPressure_col_ad (czp)', 'START', verb_opt=2)
3665
3666 Vcode = col_getVco(columnIncRef)%vcode
3667 if (Vcode == 5005 .or. Vcode == 5002) then
3668 if ( .not. col_varExist(columnInc,'P_*') ) then
3669 call utl_abort('calcPressure_col_ad (czp): for vcode 5xxx, variables P_M and P_T must be allocated in column')
3670 end if
3671 if ( .not. col_varExist(columnInc,'P0') ) then
3672 call utl_abort('calcPressure_col_ad (czp): for vcode 5xxx, variable P0 must be allocated in column')
3673 end if
3674 call calcPressure_col_ad_vcode5xxx
3675 else if (Vcode == 21001) then
3676 ! Development notes (@mad001)
3677 ! probably some some gsv_varExist(statevector,.) needed for GEM-H
3678 call calcPressure_col_ad_vcode2100x
3679 else
3680 call utl_abort('calcPressure_col_ad (czp): not implemented')
3681 end if
3682
3683 call msg('calcPressure_col_ad (czp)', 'END', verb_opt=2)
3684
3685 contains
3686 !---------------------------------------------------------
3687 ! calcPressure_col_ad_vcode2100x
3688 !---------------------------------------------------------
3689 subroutine calcPressure_col_ad_vcode2100x
3690 implicit none
3691
3692 call utl_abort('calcPressure_col_ad (czp): vcode 21001 not implemented yet')
3693
3694 end subroutine calcPressure_col_ad_vcode2100x
3695
3696 !---------------------------------------------------------
3697 ! calcPressure_col_ad_vcode5xxx
3698 !---------------------------------------------------------
3699 subroutine calcPressure_col_ad_vcode5xxx
3700 implicit none
3701
3702 ! Locals:
3703 real(8) :: Psfc
3704 real(8), pointer :: delPsfc(:,:), PsfcRef(:,:)
3705 real(8), pointer :: delP_T(:,:), delP_M(:,:)
3706 real(8), pointer :: dP_dPsfc_T(:), dP_dPsfc_M(:)
3707 integer :: status, colIndex
3708 integer :: lev_M, lev_T, nlev_T, nlev_M, numColumns
3709
3710 call msg('calcPressure_col_ad_vcode5xxx (czp)', 'START', verb_opt=4)
3711
3712 nullify(delPsfc)
3713 nullify(PsfcRef)
3714 nullify(delP_T)
3715 nullify(delP_M)
3716 nullify(dP_dPsfc_T)
3717 nullify(dP_dPsfc_M)
3718
3719 delP_M => col_getAllColumns(columnInc,'P_M')
3720 delP_T => col_getAllColumns(columnInc,'P_T')
3721 delPsfc => col_getAllColumns(columnInc,'P0')
3722 PsfcRef => col_getAllColumns(columnIncRef,'P0')
3723
3724 nlev_T = col_getNumLev(columnInc,'TH')
3725 nlev_M = col_getNumLev(columnInc,'MM')
3726 numColumns = col_getNumCol(columnInc)
3727
3728 do colIndex = 1, numColumns
3729
3730 Psfc = PsfcRef(1,colIndex)
3731
3732 ! dP_dPsfc_M
3733 nullify(dP_dPsfc_M)
3734 status = vgd_dpidpis(columnInc%vco%vgrid, &
3735 columnInc%vco%ip1_M, &
3736 dP_dPsfc_M, &
3737 Psfc)
3738 if( status .ne. VGD_OK ) then
3739 call utl_abort('calcPressure_col_ad (czp): ERROR with vgd_dpidpis')
3740 end if
3741 ! calculate delP_M
3742 do lev_M = 1, nlev_M
3743 delPsfc(1,colIndex) = delPsfc(1,colIndex) + &
3744 dP_dPsfc_M(lev_M) * delP_M(lev_M,colIndex)
3745 end do
3746 deallocate(dP_dPsfc_M)
3747
3748 ! dP_dPsfc_T
3749 nullify(dP_dPsfc_T)
3750 status = vgd_dpidpis(columnInc%vco%vgrid, &
3751 columnInc%vco%ip1_T, &
3752 dP_dPsfc_T, &
3753 Psfc)
3754 if( status .ne. VGD_OK ) then
3755 call utl_abort('calcPressure_col_ad (czp): ERROR with vgd_dpidpis')
3756 end if
3757 ! calculate delP_T
3758 do lev_T = 1, nlev_T
3759 delPsfc(1,colIndex) = &
3760 delPsfc(1,colIndex) + dP_dPsfc_T(lev_T) * delP_T(lev_T,colIndex)
3761 end do
3762 deallocate(dP_dPsfc_T)
3763
3764 end do
3765 call msg('calcPressure_col_ad_vcode5xxx (czp)', 'END', verb_opt=4)
3766 end subroutine calcPressure_col_ad_vcode5xxx
3767
3768 end subroutine calcPressure_col_ad
3769
3770 !---------------------------------------------------------------------
3771 ! subroutines wrapping vgd_levels and vgd_dpidpis queries
3772 !---------------------------------------------------------------------
3773
3774 !---------------------------------------------------------
3775 ! fetch3DLevels_r8
3776 !---------------------------------------------------------
3777 subroutine fetch3DLevels_r8(vco, sfcFld, sfcFldLS_opt, fldM_opt, fldT_opt)
3778 !
3779 !:Purpose: Main vgd_levels wrapper for field queries. Return vertical coordinate
3780 ! fields for both momentum and thermodynamic levels; real(8) flavor.
3781 !
3782 implicit none
3783
3784 ! Arguments:
3785 type(struct_vco), intent(in) :: vco ! Vertical descriptor
3786 real(8), intent(in) :: sfcFld(:,:) ! Surface field reference for coordinate
3787 real(8), optional, intent(in) :: sfcFldLS_opt(:,:)! Large scale surface field reference for coordinate (SLEVE)
3788 real(8), optional, pointer, intent(inout) :: fldM_opt(:,:,:) ! Momemtum levels field
3789 real(8), optional, pointer, intent(inout) :: fldT_opt(:,:,:) ! Thermodynamic levels field
3790
3791 ! Locals:
3792 integer :: status
3793
3794 if ( minval(sfcFld) <=0 ) then
3795 if ( vco%vcode == 21001 ) then
3796 call msg('fetch3DLevels_r8','WARNING negative surface height reference')
3797 else
3798 call utl_abort('fetch3DLevels_r8: negative surface reference')
3799 end if
3800 end if
3801
3802 if (present(fldM_opt)) then
3803 nullify(fldM_opt)
3804 if (vco%vcode == 5100) then
3805 if (.not. present(sfcFldLS_opt) ) then
3806 call utl_abort('fetch3DLevels_r8: require sfcFldLS_opt for SLEVE')
3807 end if
3808 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_M, &
3809 levels=fldM_opt, &
3810 sfc_field=sfcFld, sfc_field_ls=sfcFldLS_opt, &
3811 in_log=.false.)
3812 else
3813 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_M, &
3814 levels=fldM_opt, sfc_field=sfcFld, in_log=.false.)
3815 end if
3816 if ( status .ne. VGD_OK ) then
3817 call utl_abort('fetch3DLevels_r8: ERROR with vgd_levels (momentum levels)')
3818 end if
3819 end if
3820
3821 if (present(fldT_opt)) then
3822 nullify(fldT_opt)
3823 if (vco%vcode == 5100) then
3824 if (.not. present(sfcFldLS_opt) ) then
3825 call utl_abort('fetch3DLevels_r8: require sfcFldLS_opt for SLEVE')
3826 end if
3827 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_T, &
3828 levels=fldT_opt, &
3829 sfc_field=sfcFld, sfc_field_ls=sfcFldLS_opt, &
3830 in_log=.false.)
3831 else
3832 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_T, &
3833 levels=fldT_opt, sfc_field=sfcFld, in_log=.false.)
3834 end if
3835 if ( status .ne. VGD_OK ) then
3836 call utl_abort('fetch3DLevels_r8: ERROR with vgd_levels (thermodynamic levels)')
3837 end if
3838 end if
3839 end subroutine fetch3DLevels_r8
3840
3841 !---------------------------------------------------------
3842 ! fetch3DLevels_r4
3843 !---------------------------------------------------------
3844 subroutine fetch3DLevels_r4(vco, sfcFld, sfcFldLS_opt, fldM_opt, fldT_opt)
3845 !
3846 !:Purpose: Main vgd_levels wrapper for field query. Return vertical coordinate
3847 ! fields for both momentum and thermodynamic levels; real(4) flavor.
3848 !
3849 implicit none
3850
3851 ! Arguments:
3852 type(struct_vco), intent(in) :: vco ! Vertical descriptor
3853 real(4), intent(in) :: sfcFld(:,:) ! Surface field reference for coordinate
3854 real(4), optional, intent(in) :: sfcFldLS_opt(:,:)! Large scale surface field reference for coordinate (SLEVE)
3855 real(4), optional, pointer, intent(inout) :: fldM_opt(:,:,:) ! Momemtum levels field
3856 real(4), optional, pointer, intent(inout) :: fldT_opt(:,:,:) ! Thermodynamic levels field
3857
3858 ! Locals:
3859 integer :: status
3860
3861 if ( minval(sfcFld) <=0 ) then
3862 if ( vco%vcode == 21001 ) then
3863 call msg('fetch3DLevels_r4','WARNING negative surface height reference')
3864 else
3865 call utl_abort('fetch3DLevels_r4: negative surface reference')
3866 end if
3867 end if
3868
3869 if (present(fldM_opt)) then
3870 nullify(fldM_opt)
3871 if (vco%vcode == 5100) then
3872 if (.not. present(sfcFldLS_opt) ) then
3873 call utl_abort('fetch3DLevels_r4: require sfcFldLS_opt for SLEVE')
3874 end if
3875 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_M, &
3876 levels=fldM_opt, &
3877 sfc_field=sfcFld, sfc_field_ls=sfcFldLS_opt, &
3878 in_log=.false.)
3879 else
3880 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_M, &
3881 levels=fldM_opt, sfc_field=sfcFld, in_log=.false.)
3882 end if
3883 if ( status .ne. VGD_OK ) then
3884 call utl_abort('fetch3DLevels_r4: ERROR with vgd_levels (momentum levels)')
3885 end if
3886 end if
3887
3888 if (present(fldT_opt)) then
3889 nullify(fldT_opt)
3890 if (vco%vcode == 5100) then
3891 if (.not. present(sfcFldLS_opt) ) then
3892 call utl_abort('fetch3DLevels_r4: require sfcFldLS_opt for SLEVE')
3893 end if
3894 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_T, &
3895 levels=fldT_opt, &
3896 sfc_field=sfcFld, sfc_field_ls=sfcFldLS_opt, &
3897 in_log=.false.)
3898 else
3899 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_T, &
3900 levels=fldT_opt, sfc_field=sfcFld, in_log=.false.)
3901 end if
3902 if ( status .ne. VGD_OK ) then
3903 call utl_abort('fetch3DLevels_r4: ERROR with vgd_levels (thermodynamic levels)')
3904 end if
3905 end if
3906 end subroutine fetch3DLevels_r4
3907
3908 !---------------------------------------------------------
3909 ! fetch1DLevels_r8
3910 !---------------------------------------------------------
3911 subroutine fetch1DLevels_r8(vco, sfcValue, profM_opt, profT_opt)
3912 !
3913 !:Purpose: Main vgd_levels wrapper for profile query. Return vertical coordinate
3914 ! profile for both momentum and thermodynamic levels; real(8) flavor.
3915 !
3916 implicit none
3917
3918 ! Arguments:
3919 type(struct_vco), intent(in) :: vco ! Vertical descriptor
3920 real(8), intent(in) :: sfcValue ! Surface field reference for coordinate
3921 real(8), pointer, optional, intent(inout) :: profM_opt(:) ! Momemtum levels profile
3922 real(8), pointer, optional, intent(inout) :: profT_opt(:) ! Thermodynamic levels profile
3923
3924 ! Locals:
3925 integer :: status
3926
3927 if ( sfcValue <=0 ) then
3928 if ( vco%vcode == 21001 ) then
3929 call msg('fetch1DLevels_r8','WARNING negative surface height reference')
3930 else
3931 call utl_abort('fetch1DLevels_r8: negative surface reference')
3932 end if
3933 end if
3934
3935 if (present(profM_opt)) then
3936 nullify(profM_opt)
3937 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_M, &
3938 levels=profM_opt, sfc_field=sfcValue, in_log=.false.)
3939 if ( status .ne. VGD_OK ) then
3940 call utl_abort('fetch1DLevels_r8: ERROR with vgd_levels (momentum levels)')
3941 end if
3942 end if
3943
3944 if (present(profT_opt)) then
3945 nullify(profT_opt)
3946 status = vgd_levels(vco%vgrid, ip1_list=vco%ip1_T, &
3947 levels=profT_opt, sfc_field=sfcValue, in_log=.false.)
3948 if ( status .ne. VGD_OK ) then
3949 call utl_abort('fetch1DLevels_r8: ERROR with vgd_levels (thermodynamic levels)')
3950 end if
3951 end if
3952 end subroutine fetch1DLevels_r8
3953
3954 !---------------------------------------------------------
3955 ! fetch1DdPdPs_r8
3956 !---------------------------------------------------------
3957 subroutine fetch1DdPdPs_r8(vco, sfcValue, profM_opt, profT_opt)
3958 !
3959 !:Purpose: Main vgd_levels wrapper for iderivative profile query. Return vertical
3960 ! coordinate profile for both momentum and thermodynamic levels;
3961 ! real(8) flavor.
3962 !
3963 implicit none
3964
3965 ! Arguments:
3966 type(struct_vco), intent(in) :: vco ! Vertical descriptor
3967 real(8), intent(in) :: sfcValue ! Surface field reference for coordinate
3968 real(8), pointer, optional, intent(inout) :: profM_opt(:) ! Momemtum levels profile
3969 real(8), pointer, optional, intent(inout) :: profT_opt(:) ! Thermodynamic levels profile
3970
3971 ! Locals:
3972 integer :: status
3973
3974 if ( sfcValue <=0 ) then
3975 if ( vco%vcode == 21001 ) then
3976 call msg('fetch1DdPdPs_r8','WARNING negative surface height reference')
3977 else
3978 call utl_abort('fetch1DdPdPs_r8: negative surface reference')
3979 end if
3980 end if
3981
3982 if (present(profM_opt)) then
3983 nullify(profM_opt)
3984 status = vgd_dpidpis(vco%vgrid, vco%ip1_M, profM_opt, sfcValue)
3985 if ( status .ne. VGD_OK ) then
3986 call utl_abort('fetch1DdPdPs_r8: ERROR with vgd_dpidpis (momentum levels)')
3987 end if
3988 end if
3989
3990 if (present(profT_opt)) then
3991 nullify(profT_opt)
3992 status = vgd_dpidpis(vco%vgrid, vco%ip1_T, profT_opt, sfcValue)
3993 if ( status .ne. VGD_OK ) then
3994 call utl_abort('fetch1DdPdPs_r8: ERROR with vgd_dpidpis (thermodynamic levels)')
3995 end if
3996 end if
3997 end subroutine fetch1DdPdPs_r8
3998
3999 !---------------------------------------------------------------------
4000 ! other vertical coordinate related functions and subroutines
4001 !---------------------------------------------------------------------
4002
4003 !--------------------------------------------------------------------------
4004 ! czp_ensureCompatibleTops
4005 !--------------------------------------------------------------------------
4006 subroutine czp_ensureCompatibleTops(vco_sourceGrid,vco_destGrid)
4007 !
4008 !:Purpose: This function checks if the top of a destination grid
4009 ! is ~equal or lower than the top of the source grid
4010 ! the code aborts if this is not the case
4011 !
4012 implicit none
4013
4014 ! arguments:
4015 type(struct_vco), pointer, intent(in) :: vco_sourceGrid ! vertical coordinate source grid
4016 type(struct_vco), pointer, intent(in) :: vco_destGrid ! vertical coordinate destination grid
4017
4018 ! locals:
4019 integer :: nAbove, numLevSource, numLevDest
4020 real(8) :: sourceModelTop
4021 real(8) :: pSfc(1,1), pSfcLS(1,1)
4022 real(8), pointer :: sourcePressureLevels(:,:,:)
4023 real(8), pointer :: destPressureLevels(:,:,:)
4024
4025 nullify(sourcePressureLevels)
4026 nullify(destPressureLevels)
4027
4028 numLevSource = vco_getNumLev(vco_sourceGrid,'MM')
4029 numLevDest = vco_getNumLev(vco_destGrid,'MM')
4030 if (numLevSource == 1 .and. numLevDest == 1) then
4031 write(*,*) 'czp_ensureCompatibleTops: both grids only have 1 level, skip the test'
4032 return
4033 end if
4034
4035 ! dummy pressure value
4036 pSfc(1,1) = 100.0D3 !100 kPa
4037 pSfcLS(1,1) = 100.0D3 !100 kPa
4038
4039 ! pressure on momentum levels of source grid
4040 if (vco_sourceGrid%vcode == 5100) then
4041 call fetch3DLevels_r8(vco_sourceGrid, pSfc, sfcFldLS_opt=pSfcLS, &
4042 fldM_opt=sourcePressureLevels)
4043 else
4044 call fetch3DLevels_r8(vco_sourceGrid, pSfc, fldM_opt=sourcePressureLevels)
4045 end if
4046 ! pressure on momentum levels of destination grid
4047 if (vco_destGrid%vcode == 5100) then
4048 call fetch3DLevels_r8(vco_destGrid, pSfc, sfcFldLS_opt=pSfcLS, &
4049 fldM_opt=destPressureLevels)
4050 else
4051 call fetch3DLevels_r8(vco_destGrid, pSfc, fldM_opt=destPressureLevels)
4052 end if
4053
4054 ! count number of levels where output grid is higher than input grid
4055 sourceModelTop = sourcePressureLevels(1,1,1)
4056 nAbove=0
4057 do while (sourceModelTop > destPressureLevels(1,1,nAbove+1))
4058 nAbove = nAbove + 1
4059 end do
4060
4061 ! Destination grid has "nAbove" levels above source grid; tolerate one
4062 if ( nAbove > 1 ) then
4063 write(*,*) 'czp_ensureCompatibleTops: numLevSource/Dest = ', numLevSource, numLevDest
4064 write(*,*) 'czp_ensureCompatibleTops: sourcePressureLevels = ', sourcePressureLevels(1,1,:)
4065 write(*,*) 'czp_ensureCompatibleTops: destPressureLevels = ', destPressureLevels(1,1,:)
4066 call utl_abort('czp_ensureCompatibleTops: top of destination grid more than one level higher than top of source grid')
4067 end if
4068
4069 deallocate(sourcePressureLevels)
4070 deallocate(destPressureLevels)
4071
4072 end subroutine czp_ensureCompatibleTops
4073
4074 !---------------------------------------------------------------------
4075 ! helper private functions and subroutines
4076 !---------------------------------------------------------------------
4077
4078 !---------------------------------------------------------
4079 ! calcHeightCoeff_gsv
4080 !---------------------------------------------------------
4081 subroutine calcHeightCoeff_gsv(statevector)
4082 !
4083 !:Purpose: Calculating the coefficients of height for
4084 ! czp_calcHeight_tl/czp_calcHeight_ad
4085 !
4086 implicit none
4087
4088 ! Arguments:
4089 type(struct_gsv), intent(in) :: statevector
4090
4091 ! Locals:
4092 integer :: lev_T,nlev_M,nlev_T,numStep,stepIndex,latIndex,lonIndex,Vcode
4093 real(8) :: hu,tt,Pr,height_T,cmp,cmp_TT,cmp_HU,cmp_P0_1,cmp_P0_2,ratioP1
4094 real(4) :: lat_4
4095 real(8) :: Rgh, sLat, lat_8
4096 real(8), pointer :: hu_ptr(:,:,:,:),tt_ptr(:,:,:,:)
4097 real(8), pointer :: P_T_ptr(:,:,:,:),P_M_ptr(:,:,:,:)
4098 real(8), pointer :: height_T_ptr(:,:,:,:)
4099 type(struct_vco), pointer :: vco
4100 logical, save :: firstTimeHeightCoeff_gsv = .true.
4101
4102 if ( .not. firstTimeHeightCoeff_gsv ) return
4103
4104 call msg('calcHeightCoeff_gsv (czp)', 'START', verb_opt=2)
4105
4106 ! initialize and save coefficients for increased efficiency
4107 ! (assumes no relinearization)
4108 firstTimeHeightCoeff_gsv = .false.
4109
4110 vco => gsv_getVco(statevector)
4111 Vcode = vco%vcode
4112
4113 nlev_T = gsv_getNumLev(statevector,'TH')
4114 nlev_M = gsv_getNumLev(statevector,'MM')
4115 numStep = statevector%numstep
4116
4117 ! saved arrays
4118 allocate(coeff_M_TT_gsv(&
4119 statevector%myLonBeg:statevector%myLonEnd, &
4120 statevector%myLatBeg:statevector%myLatEnd, &
4121 nlev_T,numStep))
4122 allocate(coeff_M_HU_gsv(&
4123 statevector%myLonBeg:statevector%myLonEnd, &
4124 statevector%myLatBeg:statevector%myLatEnd, &
4125 nlev_T,numStep))
4126 allocate(coeff_M_P0_delPM_gsv(&
4127 statevector%myLonBeg:statevector%myLonEnd, &
4128 statevector%myLatBeg:statevector%myLatEnd, &
4129 nlev_T,numStep))
4130 allocate(coeff_M_P0_dP_delPT_gsv(&
4131 statevector%myLonBeg:statevector%myLonEnd, &
4132 statevector%myLatBeg:statevector%myLatEnd, &
4133 nlev_T,numStep))
4134 allocate(coeff_M_P0_dP_delP0_gsv(&
4135 statevector%myLonBeg:statevector%myLonEnd, &
4136 statevector%myLatBeg:statevector%myLatEnd, &
4137 nlev_T,numStep))
4138
4139 allocate(coeff_T_TT_gsv(&
4140 statevector%myLonBeg:statevector%myLonEnd, &
4141 statevector%myLatBeg:statevector%myLatEnd, &
4142 numStep))
4143 allocate(coeff_T_HU_gsv(&
4144 statevector%myLonBeg:statevector%myLonEnd, &
4145 statevector%myLatBeg:statevector%myLatEnd, &
4146 numStep))
4147 allocate(coeff_T_P0_delP1_gsv(&
4148 statevector%myLonBeg:statevector%myLonEnd, &
4149 statevector%myLatBeg:statevector%myLatEnd, &
4150 numStep))
4151 allocate(coeff_T_P0_dP_delPT_gsv(&
4152 statevector%myLonBeg:statevector%myLonEnd, &
4153 statevector%myLatBeg:statevector%myLatEnd, &
4154 numStep))
4155 allocate(coeff_T_P0_dP_delP0_gsv(&
4156 statevector%myLonBeg:statevector%myLonEnd, &
4157 statevector%myLatBeg:statevector%myLatEnd, &
4158 numStep))
4159
4160 coeff_M_TT_gsv(:,:,:,:) = 0.0D0
4161 coeff_M_HU_gsv(:,:,:,:) = 0.0D0
4162
4163 coeff_M_P0_delPM_gsv(:,:,:,:) = 0.0D0
4164
4165 coeff_M_P0_dP_delPT_gsv(:,:,:,:) = 0.0D0
4166 coeff_M_P0_dP_delP0_gsv(:,:,:,:) = 0.0D0
4167
4168 coeff_T_TT_gsv(:,:,:) = 0.0D0
4169 coeff_T_HU_gsv(:,:,:) = 0.0D0
4170
4171 coeff_T_P0_delP1_gsv(:,:,:) = 0.0D0
4172
4173 coeff_T_P0_dP_delPT_gsv(:,:,:) = 0.0D0
4174 coeff_T_P0_dP_delP0_gsv(:,:,:) = 0.0D0
4175
4176 call gsv_getField(statevector,hu_ptr,'HU')
4177 call gsv_getField(statevector,tt_ptr,'TT')
4178 call gsv_getField(statevector,P_T_ptr,'P_T')
4179 call gsv_getField(statevector,P_M_ptr,'P_M')
4180 call gsv_getField(statevector,height_T_ptr,'Z_T')
4181
4182 if_calcHeightCoeff_gsv_vcodes : if (Vcode == 5002) then
4183
4184 do stepIndex = 1, numStep
4185 do lev_T = 1, (nlev_T-1)
4186 do latIndex = statevector%myLatBeg, statevector%myLatEnd
4187 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
4188 if ( lev_T == 1 ) then
4189 ! compute height coefficients on only the top thermo level
4190 ratioP1 = log( P_M_ptr(lonIndex,latIndex,1,stepIndex) / &
4191 P_T_ptr(lonIndex,latIndex,1,stepIndex) )
4192 hu = max(hu_ptr(lonIndex,latIndex,1,stepIndex),&
4193 MPC_MINIMUM_HU_R8)
4194 tt = tt_ptr(lonIndex,latIndex,1,stepIndex)
4195 Pr = P_T_ptr(lonIndex,latIndex,1,stepIndex)
4196 height_T = height_T_ptr(lonIndex,latIndex,1,stepIndex)
4197
4198 cmp = gpscompressibility(Pr,tt,hu)
4199 cmp_TT = gpscompressibility_TT(Pr,tt,hu)
4200 cmp_HU = gpscompressibility_HU(Pr,tt,hu)
4201 cmp_P0_1 = gpscompressibility_P0_1(Pr,tt,hu,1.0d0)
4202 cmp_P0_2 = gpscompressibility_P0_2(Pr,tt,hu)
4203
4204 ! Gravity acceleration
4205 lat_4 = statevector%hco%lat2d_4(lonIndex,latIndex)
4206 lat_8 = real(lat_4,8)
4207 sLat = sin(lat_8)
4208 Rgh = phf_gravityalt(sLat, height_T)
4209
4210 coeff_T_TT_gsv(lonIndex,latIndex,stepIndex) = &
4211 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_fottva(hu,1.0D0) * cmp + &
4212 phf_fotvt8(tt,hu) * cmp_TT) * ratioP1
4213
4214 coeff_T_HU_gsv(lonIndex,latIndex,stepIndex) = &
4215 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_folnqva(hu,tt,1.0d0) / &
4216 hu * cmp + phf_fotvt8(tt,hu) * cmp_HU) * ratioP1
4217
4218 coeff_T_P0_delP1_gsv(lonIndex,latIndex,stepIndex) = &
4219 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp
4220
4221 coeff_T_P0_dP_delPT_gsv(lonIndex,latIndex,stepIndex) = &
4222 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_1 * &
4223 ratioP1
4224
4225 coeff_T_P0_dP_delP0_gsv(lonIndex,latIndex,stepIndex) = &
4226 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_2 * &
4227 ratioP1
4228 else
4229 ! compute height coefficients on momentum levels
4230 ratioP1 = log( P_M_ptr(lonIndex,latIndex,lev_T ,stepIndex) / &
4231 P_M_ptr(lonIndex,latIndex,lev_T-1,stepIndex) )
4232 hu = max( hu_ptr(lonIndex,latIndex,lev_T,stepIndex),&
4233 MPC_MINIMUM_HU_R8)
4234 tt = tt_ptr(lonIndex,latIndex,lev_T,stepIndex)
4235 Pr = P_T_ptr(lonIndex,latIndex,lev_T,stepIndex)
4236 height_T = height_T_ptr(lonIndex,latIndex,lev_T,stepIndex)
4237
4238 cmp = gpscompressibility(Pr,tt,hu)
4239 cmp_TT = gpscompressibility_TT(Pr,tt,hu)
4240 cmp_HU = gpscompressibility_HU(Pr,tt,hu)
4241 cmp_P0_1 = gpscompressibility_P0_1(Pr,tt,hu,1.0d0)
4242 cmp_P0_2 = gpscompressibility_P0_2(Pr,tt,hu)
4243
4244 ! Gravity acceleration
4245 lat_4 = statevector%hco%lat2d_4(lonIndex,latIndex)
4246 lat_8 = real(lat_4,8)
4247 sLat = sin(lat_8)
4248 Rgh = phf_gravityalt(sLat, height_T)
4249
4250 coeff_M_TT_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4251 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_fottva(hu,1.0D0) * cmp + &
4252 phf_fotvt8(tt,hu) * cmp_TT) * ratioP1
4253
4254 coeff_M_HU_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4255 (MPC_RGAS_DRY_AIR_R8 / RgH) * (phf_folnqva(hu,tt,1.0d0) / hu * &
4256 cmp + phf_fotvt8(tt,hu) * cmp_HU) * ratioP1
4257
4258 coeff_M_P0_delPM_gsv (lonIndex,latIndex,lev_T,stepIndex) = &
4259 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp
4260
4261 coeff_M_P0_dP_delPT_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4262 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_1 * &
4263 ratioP1
4264
4265 coeff_M_P0_dP_delP0_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4266 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_2 * &
4267 ratioP1
4268 end if
4269 end do
4270 end do
4271 end do
4272 end do
4273
4274 else if (Vcode == 5005) then if_calcHeightCoeff_gsv_vcodes
4275
4276 do stepIndex = 1, numStep
4277 do lev_T = 1, (nlev_T-1)
4278 do latIndex = statevector%myLatBeg, statevector%myLatEnd
4279 do lonIndex = statevector%myLonBeg, statevector%myLonEnd
4280 ! compute height coefficients on momentum levels
4281 ratioP1 = log( P_M_ptr(lonIndex,latIndex,lev_T+1,stepIndex) / &
4282 P_M_ptr(lonIndex,latIndex,lev_T ,stepIndex) )
4283 hu = max( hu_ptr(lonIndex,latIndex,lev_T,stepIndex),&
4284 MPC_MINIMUM_HU_R8)
4285 tt = tt_ptr(lonIndex,latIndex,lev_T,stepIndex)
4286 Pr = P_T_ptr(lonIndex,latIndex,lev_T,stepIndex)
4287 height_T = height_T_ptr(lonIndex,latIndex,lev_T,stepIndex)
4288
4289 cmp = gpscompressibility(Pr,tt,hu)
4290 cmp_TT = gpscompressibility_TT(Pr,tt,hu)
4291 cmp_HU = gpscompressibility_HU(Pr,tt,hu)
4292 cmp_P0_1 = gpscompressibility_P0_1(Pr,tt,hu,1.0d0)
4293 cmp_P0_2 = gpscompressibility_P0_2(Pr,tt,hu)
4294
4295 ! Gravity acceleration
4296 lat_4 = statevector%hco%lat2d_4(lonIndex,latIndex)
4297 lat_8 = real(lat_4,8)
4298 sLat = sin(lat_8)
4299 Rgh = phf_gravityalt(sLat, height_T)
4300
4301 coeff_M_TT_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4302 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_fottva(hu,1.0D0) * cmp + &
4303 phf_fotvt8(tt,hu) * cmp_TT) * ratioP1
4304
4305 coeff_M_HU_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4306 (MPC_RGAS_DRY_AIR_R8 / RgH) * (phf_folnqva(hu,tt,1.0d0) / hu * &
4307 cmp + phf_fotvt8(tt,hu) * cmp_HU) * ratioP1
4308
4309 coeff_M_P0_delPM_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4310 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp
4311
4312 coeff_M_P0_dP_delPT_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4313 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_1 * &
4314 ratioP1
4315
4316 coeff_M_P0_dP_delP0_gsv(lonIndex,latIndex,lev_T,stepIndex) = &
4317 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_2 * &
4318 ratioP1
4319 end do
4320 end do
4321 end do
4322 end do
4323
4324 else if_calcHeightCoeff_gsv_vcodes
4325
4326 call utl_abort('calcHeightCoeff_gsv (czp): only vcode 5002 and 5005 implemented')
4327
4328 end if if_calcHeightCoeff_gsv_vcodes
4329
4330 call msg('calcHeightCoeff_gsv (czp)', 'END', verb_opt=2)
4331 end subroutine calcHeightCoeff_gsv
4332
4333 !---------------------------------------------------------
4334 ! calcHeightCoeff_col
4335 !---------------------------------------------------------
4336 subroutine calcHeightCoeff_col(column)
4337 !
4338 !:Purpose: Calculating the coefficients of height for
4339 ! czp_calcHeight_tl/czp_calcHeight_ad
4340 !
4341 implicit none
4342
4343 ! Arguments:
4344 type(struct_columnData), intent(in) :: column
4345
4346 ! Locals:
4347 integer :: lev_T,nlev_M,nlev_T,numColumns,colIndex,Vcode
4348 real(8) :: hu,tt,Pr,height_T,cmp,cmp_TT,cmp_HU,cmp_P0_1,cmp_P0_2,ratioP1
4349 real(8) :: Rgh, sLat, lat_8
4350 real(8), pointer :: hu_ptr(:,:),tt_ptr(:,:)
4351 real(8), pointer :: P_T_ptr(:,:),P_M_ptr(:,:)
4352 real(8), pointer :: height_T_ptr(:,:)
4353 type(struct_vco), pointer :: vco
4354 logical, save :: firstTimeHeightCoeff_col = .true.
4355
4356 if ( .not. firstTimeHeightCoeff_col ) return
4357
4358 call msg('calcHeightCoeff_col (czp)', 'START', verb_opt=2)
4359
4360 ! initialize and save coefficients for increased efficiency
4361 ! (assumes no relinearization)
4362 firstTimeHeightCoeff_col = .false.
4363
4364 vco => col_getVco(column)
4365 Vcode = vco%vcode
4366
4367 nlev_T = col_getNumLev(column,'TH')
4368 nlev_M = col_getNumLev(column,'MM')
4369 numColumns = col_getNumCol(column)
4370
4371 ! saved arrays
4372 allocate(coeff_M_TT_col (nlev_T,numColumns))
4373 allocate(coeff_M_HU_col (nlev_T,numColumns))
4374 allocate(coeff_M_P0_delPM_col (nlev_T,numColumns))
4375 allocate(coeff_M_P0_dP_delPT_col(nlev_T,numColumns))
4376 allocate(coeff_M_P0_dP_delP0_col(nlev_T,numColumns))
4377
4378 allocate(coeff_T_TT_col (numColumns))
4379 allocate(coeff_T_HU_col (numColumns))
4380 allocate(coeff_T_P0_delP1_col (numColumns))
4381 allocate(coeff_T_P0_dP_delPT_col(numColumns))
4382 allocate(coeff_T_P0_dP_delP0_col(numColumns))
4383
4384 coeff_M_TT_col(:,:) = 0.0D0
4385 coeff_M_HU_col(:,:) = 0.0D0
4386
4387 coeff_M_P0_delPM_col(:,:) = 0.0D0
4388
4389 coeff_M_P0_dP_delPT_col(:,:) = 0.0D0
4390 coeff_M_P0_dP_delP0_col(:,:) = 0.0D0
4391
4392 coeff_T_TT_col(:) = 0.0D0
4393 coeff_T_HU_col(:) = 0.0D0
4394
4395 coeff_T_P0_delP1_col(:) = 0.0D0
4396
4397 coeff_T_P0_dP_delPT_col(:) = 0.0D0
4398 coeff_T_P0_dP_delP0_col(:) = 0.0D0
4399
4400 hu_ptr => col_getAllColumns(column,'HU')
4401 tt_ptr => col_getAllColumns(column,'TT')
4402 P_T_ptr => col_getAllColumns(column,'P_T')
4403 P_M_ptr => col_getAllColumns(column,'P_M')
4404 height_T_ptr => col_getAllColumns(column,'Z_T')
4405
4406 if_calcHeightCoeff_col_vcodes : if (Vcode == 5002) then
4407
4408 do colIndex = 1, numColumns
4409 do lev_T = 1, (nlev_T-1)
4410 if ( lev_T == 1 ) then
4411 ! compute height coefficients on only the top thermo level
4412 ratioP1 = log( P_M_ptr(1,colIndex) / &
4413 P_T_ptr(1,colIndex) )
4414 hu = max(hu_ptr(1,colIndex),MPC_MINIMUM_HU_R8)
4415 tt = tt_ptr(1,colIndex)
4416 Pr = P_T_ptr(1,colIndex)
4417 height_T = height_T_ptr(1,colIndex)
4418
4419 cmp = gpscompressibility(Pr,tt,hu)
4420 cmp_TT = gpscompressibility_TT(Pr,tt,hu)
4421 cmp_HU = gpscompressibility_HU(Pr,tt,hu)
4422 cmp_P0_1 = gpscompressibility_P0_1(Pr,tt,hu,1.0d0)
4423 cmp_P0_2 = gpscompressibility_P0_2(Pr,tt,hu)
4424
4425 ! Gravity acceleration
4426 lat_8 = column%lat(colIndex)
4427 sLat = sin(lat_8)
4428 Rgh = phf_gravityalt(sLat, height_T)
4429
4430 coeff_T_TT_col(colIndex) = &
4431 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_fottva(hu,1.0D0) * cmp + &
4432 phf_fotvt8(tt,hu) * cmp_TT) * ratioP1
4433
4434 coeff_T_HU_col(colIndex) = &
4435 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_folnqva(hu,tt,1.0d0) / hu * &
4436 cmp + phf_fotvt8(tt,hu) * cmp_HU) * ratioP1
4437
4438 coeff_T_P0_delP1_col (colIndex) = &
4439 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp
4440
4441 coeff_T_P0_dP_delPT_col(colIndex) = &
4442 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_1 * ratioP1
4443
4444 coeff_T_P0_dP_delP0_col(colIndex) = &
4445 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_2 * ratioP1
4446 else
4447 ! compute height coefficients on momentum levels
4448 ratioP1 = log( P_M_ptr(lev_T ,colIndex) / &
4449 P_M_ptr(lev_T-1,colIndex) )
4450 hu = max(hu_ptr(lev_T,colIndex),MPC_MINIMUM_HU_R8)
4451 tt = tt_ptr(lev_T,colIndex)
4452 Pr = P_T_ptr(lev_T,colIndex)
4453 height_T = height_T_ptr(lev_T,colIndex)
4454
4455 cmp = gpscompressibility(Pr,tt,hu)
4456 cmp_TT = gpscompressibility_TT(Pr,tt,hu)
4457 cmp_HU = gpscompressibility_HU(Pr,tt,hu)
4458 cmp_P0_1 = gpscompressibility_P0_1(Pr,tt,hu,1.0d0)
4459 cmp_P0_2 = gpscompressibility_P0_2(Pr,tt,hu)
4460
4461 ! Gravity acceleration
4462 lat_8 = column%lat(colIndex)
4463 sLat = sin(lat_8)
4464 Rgh = phf_gravityalt(sLat, height_T)
4465
4466 coeff_M_TT_col(lev_T,colIndex) = &
4467 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_fottva(hu,1.0D0) * cmp + &
4468 phf_fotvt8(tt,hu) * cmp_TT) * ratioP1
4469
4470 coeff_M_HU_col(lev_T,colIndex) = &
4471 (MPC_RGAS_DRY_AIR_R8 / RgH) * (phf_folnqva(hu,tt,1.0d0) / hu * &
4472 cmp + phf_fotvt8(tt,hu) * cmp_HU) * ratioP1
4473
4474 coeff_M_P0_delPM_col(lev_T,colIndex) = &
4475 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp
4476
4477 coeff_M_P0_dP_delPT_col(lev_T,colIndex) = &
4478 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_1 * ratioP1
4479
4480 coeff_M_P0_dP_delP0_col(lev_T,colIndex) = &
4481 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_2 * ratioP1
4482 end if
4483 end do
4484 end do
4485
4486 else if (Vcode == 5005) then if_calcHeightCoeff_col_vcodes
4487
4488 do colIndex = 1, numColumns
4489 do lev_T = 1, (nlev_T-1)
4490 ! compute height coefficients on momentum levels
4491 ratioP1 = log( P_M_ptr(lev_T+1,colIndex) / &
4492 P_M_ptr(lev_T ,colIndex) )
4493 hu = max(hu_ptr(lev_T,colIndex),MPC_MINIMUM_HU_R8)
4494 tt = tt_ptr(lev_T,colIndex)
4495 Pr = P_T_ptr(lev_T,colIndex)
4496 height_T = height_T_ptr(lev_T,colIndex)
4497
4498 cmp = gpscompressibility(Pr,tt,hu)
4499 cmp_TT = gpscompressibility_TT(Pr,tt,hu)
4500 cmp_HU = gpscompressibility_HU(Pr,tt,hu)
4501 cmp_P0_1 = gpscompressibility_P0_1(Pr,tt,hu,1.0d0)
4502 cmp_P0_2 = gpscompressibility_P0_2(Pr,tt,hu)
4503
4504 ! Gravity acceleration
4505 lat_8 = column%lat(colIndex)
4506 sLat = sin(lat_8)
4507 Rgh = phf_gravityalt(sLat, height_T)
4508
4509 coeff_M_TT_col(lev_T,colIndex) = &
4510 (MPC_RGAS_DRY_AIR_R8 / Rgh) * (phf_fottva(hu,1.0D0) * cmp + &
4511 phf_fotvt8(tt,hu) * cmp_TT) * ratioP1
4512
4513 coeff_M_HU_col(lev_T,colIndex) = &
4514 (MPC_RGAS_DRY_AIR_R8 / RgH) * (phf_folnqva(hu,tt,1.0d0) / hu * cmp + &
4515 phf_fotvt8(tt,hu) * cmp_HU) * ratioP1
4516
4517 coeff_M_P0_delPM_col (lev_T,colIndex) = &
4518 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp
4519
4520 coeff_M_P0_dP_delPT_col(lev_T,colIndex) = &
4521 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_1 * ratioP1
4522
4523 coeff_M_P0_dP_delP0_col(lev_T,colIndex) = &
4524 (MPC_RGAS_DRY_AIR_R8 / Rgh) * phf_fotvt8(tt,hu) * cmp_P0_2 * ratioP1
4525 end do
4526 end do
4527
4528 else
4529
4530 call utl_abort('calcHeightCoeff_col (czp): only vcode 5002 and 5005 implemented')
4531
4532 end if if_calcHeightCoeff_col_vcodes
4533
4534 call msg('calcHeightCoeff_col (czp)', 'END', verb_opt=2)
4535 end subroutine calcHeightCoeff_col
4536
4537 !---------------------------------------------------------
4538 ! gpscompressibility
4539 !---------------------------------------------------------
4540 function gpscompressibility(p,t,q)
4541 implicit none
4542
4543 ! Arguments:
4544 real(8), intent(in) :: p
4545 real(8), intent(in) :: t
4546 real(8), intent(in) :: q
4547 ! Result:
4548 real(8) :: gpscompressibility
4549
4550 ! Locals:
4551 real(8), parameter :: a0= 1.58123D-6
4552 real(8), parameter :: a1=-2.9331D-8
4553 real(8), parameter :: a2= 1.1043D-10
4554 real(8), parameter :: b0= 5.707D-6
4555 real(8), parameter :: b1=-2.051D-8
4556 real(8), parameter :: c0= 1.9898D-4
4557 real(8), parameter :: c1=-2.376D-6
4558 real(8), parameter :: d = 1.83D-11
4559 real(8), parameter :: e =-0.765D-8
4560 real(8) :: x,tc,pt,tc2,x2
4561
4562 if ( t <= 0 ) call utl_abort('gpscompressibility: t <= 0')
4563
4564 x = gps_p_wa * q / (1.D0 + gps_p_wb * q)
4565 ! Estimate, from CIPM, Picard (2008)
4566 tc = t - MPC_K_C_DEGREE_OFFSET_R8
4567 pt = p / t
4568 tc2= tc * tc
4569 x2 = x * x
4570 gpscompressibility = 1.D0 - pt * &
4571 (a0+a1*tc+a2*tc2+(b0+b1*tc)*x+(c0+c1*tc)*x2) + pt*pt*(d+e*x2)
4572 end function gpscompressibility
4573
4574 !---------------------------------------------------------
4575 ! gpscompressibility_TT
4576 !---------------------------------------------------------
4577 function gpscompressibility_TT(p,t,q)
4578 implicit none
4579
4580 ! Arguments:
4581 real(8), intent(in) :: p
4582 real(8), intent(in) :: t
4583 real(8), intent(in) :: q
4584 ! Result:
4585 real(8) :: gpscompressibility_TT
4586
4587 ! Locals:
4588 real(8), parameter :: a0= 1.58123D-6
4589 real(8), parameter :: a1=-2.9331D-8
4590 real(8), parameter :: a2= 1.1043D-10
4591 real(8), parameter :: b0= 5.707D-6
4592 real(8), parameter :: b1=-2.051D-8
4593 real(8), parameter :: c0= 1.9898D-4
4594 real(8), parameter :: c1=-2.376D-6
4595 real(8), parameter :: d = 1.83D-11
4596 real(8), parameter :: e =-0.765D-8
4597 real(8) :: x,tc,pt,tc2,x2
4598 real(8) :: d_x,d_tc,d_pt,d_tc2,d_x2
4599
4600 if ( t <= 0 ) call utl_abort('gpscompressibility_TT: t <= 0')
4601
4602 x = gps_p_wa * q / (1.D0 + gps_p_wb * q)
4603 d_x = 0.0D0
4604 ! Estimate, from CIPM, Picard (2008)
4605 tc = t - MPC_K_C_DEGREE_OFFSET_R8
4606 d_tc = 1.D0
4607 pt = p / t
4608 d_pt = - p / t**2
4609 tc2= tc * tc
4610 d_tc2= 2 * tc * d_tc
4611 x2 = x * x
4612 d_x2 = 2 * x * d_x
4613 gpscompressibility_TT = &
4614 -d_pt * (a0+a1*tc+a2*tc2+(b0+b1*tc)*x+(c0+c1*tc)*x2) - pt * &
4615 (a1*d_tc + a2*d_tc2 + b1*d_tc*x + (b0+b1*tc)*d_x + c1*d_tc*x2 + &
4616 (c0+c1*tc)*d_x2) + 2*pt*d_pt*(d+e*x2) + pt*pt*e*d_x2
4617 end function gpscompressibility_TT
4618
4619 !---------------------------------------------------------
4620 ! gpscompressibility_HU
4621 !---------------------------------------------------------
4622 function gpscompressibility_HU(p,t,q)
4623 implicit none
4624
4625 ! Arguments:
4626 real(8), intent(in) :: p
4627 real(8), intent(in) :: t
4628 real(8), intent(in) :: q
4629 ! Result:
4630 real(8) :: gpscompressibility_HU
4631
4632 ! Locals:
4633 real(8), parameter :: a0= 1.58123D-6
4634 real(8), parameter :: a1=-2.9331D-8
4635 real(8), parameter :: a2= 1.1043D-10
4636 real(8), parameter :: b0= 5.707D-6
4637 real(8), parameter :: b1=-2.051D-8
4638 real(8), parameter :: c0= 1.9898D-4
4639 real(8), parameter :: c1=-2.376D-6
4640 real(8), parameter :: d = 1.83D-11
4641 real(8), parameter :: e =-0.765D-8
4642 real(8) :: x,tc,pt,tc2,x2
4643 real(8) :: d_x,d_tc,d_pt,d_tc2,d_x2
4644
4645 if ( t <= 0 ) call utl_abort('gpscompressibility_HU: t <= 0')
4646
4647 x = gps_p_wa * q / (1.D0+gps_p_wb*q)
4648 d_x = gps_p_wa * (1.0D0 / (1.D0+gps_p_wb*q) - q / (1.D0+gps_p_wb*q)**2 * gps_p_wb * 1.0D0)
4649 ! Estimate, from CIPM, Picard (2008)
4650 tc = t - MPC_K_C_DEGREE_OFFSET_R8
4651 d_tc = 0.0D0
4652 pt = p / t
4653 d_pt = 0.0D0
4654 tc2= tc * tc
4655 d_tc2= 2 * tc * d_tc
4656 x2 = x * x
4657 d_x2 = 2 * x * d_x
4658 gpscompressibility_HU = &
4659 -d_pt * (a0+a1*tc+a2*tc2+(b0+b1*tc)*x+(c0+c1*tc)*x2) - &
4660 pt * (a1*d_tc + a2*d_tc2 + b1*d_tc*x + (b0+b1*tc)*d_x + c1*d_tc*x2 + &
4661 (c0+c1*tc)*d_x2) + 2*pt*d_pt*(d+e*x2) + pt*pt*e*d_x2
4662 end function gpscompressibility_HU
4663
4664 !---------------------------------------------------------
4665 ! gpscompressibility_P0_1
4666 !---------------------------------------------------------
4667 function gpscompressibility_P0_1(p,t,q,dpdp0)
4668 ! gpscompressibility_P0_1 has dpdp0 dependency
4669 implicit none
4670
4671 ! Arguments:
4672 real(8), intent(in) :: p
4673 real(8), intent(in) :: t
4674 real(8), intent(in) :: q
4675 real(8), intent(in) :: dpdp0
4676 ! Result:
4677 real(8) :: gpscompressibility_P0_1
4678
4679 ! Locals:
4680 real(8), parameter :: a0= 1.58123D-6
4681 real(8), parameter :: a1=-2.9331D-8
4682 real(8), parameter :: a2= 1.1043D-10
4683 real(8), parameter :: b0= 5.707D-6
4684 real(8), parameter :: b1=-2.051D-8
4685 real(8), parameter :: c0= 1.9898D-4
4686 real(8), parameter :: c1=-2.376D-6
4687 real(8), parameter :: d = 1.83D-11
4688 real(8), parameter :: e =-0.765D-8
4689 real(8) :: x,tc,pt,tc2,x2
4690 real(8) :: d_x,d_tc,d_pt,d_tc2,d_x2
4691
4692 if ( t <= 0 ) call utl_abort('gpscompressibility_P0_1: t <= 0')
4693
4694 x = gps_p_wa * q / (1.D0+gps_p_wb*q)
4695 d_x = 0.0D0
4696 ! Estimate, from CIPM, Picard (2008)
4697 tc = t - MPC_K_C_DEGREE_OFFSET_R8
4698 d_tc = 0.0D0
4699 pt = p / t
4700 d_pt = dpdp0 / t
4701 tc2= tc * tc
4702 d_tc2= 2 * tc * d_tc
4703 x2 = x * x
4704 d_x2 = 2 * x * d_x
4705 gpscompressibility_P0_1 = &
4706 -d_pt * (a0+a1*tc+a2*tc2+(b0+b1*tc)*x+(c0+c1*tc)*x2) + &
4707 2*pt*d_pt*(d+e*x2)
4708 end function gpscompressibility_P0_1
4709
4710 !---------------------------------------------------------
4711 ! gpscompressibility_P0_2
4712 !---------------------------------------------------------
4713 function gpscompressibility_P0_2(p,t,q)
4714 ! gpscompressibility_P0_2 has NO dpdp0 dependency
4715 implicit none
4716
4717 ! Arguments:
4718 real(8), intent(in) :: p
4719 real(8), intent(in) :: t
4720 real(8), intent(in) :: q
4721 ! Result:
4722 real(8) :: gpscompressibility_P0_2
4723
4724 ! Locals:
4725 real(8), parameter :: a0= 1.58123D-6
4726 real(8), parameter :: a1=-2.9331D-8
4727 real(8), parameter :: a2= 1.1043D-10
4728 real(8), parameter :: b0= 5.707D-6
4729 real(8), parameter :: b1=-2.051D-8
4730 real(8), parameter :: c0= 1.9898D-4
4731 real(8), parameter :: c1=-2.376D-6
4732 real(8), parameter :: d = 1.83D-11
4733 real(8), parameter :: e =-0.765D-8
4734 real(8) :: x,tc,pt,tc2,x2
4735 real(8) :: d_x,d_tc,d_tc2,d_x2
4736
4737 if ( t <= 0 ) call utl_abort('gpscompressibility_P0_2: t <= 0')
4738
4739 x = gps_p_wa * q / (1.D0+gps_p_wb*q)
4740 d_x = 0.0D0
4741 ! Estimate, from CIPM, Picard (2008)
4742 tc = t - MPC_K_C_DEGREE_OFFSET_R8
4743 d_tc = 0.0D0
4744 pt = p / t
4745 tc2= tc * tc
4746 d_tc2= 2 * tc * d_tc
4747 x2 = x * x
4748 d_x2 = 2 * x * d_x
4749 gpscompressibility_P0_2 = &
4750 -pt * (a1*d_tc + a2*d_tc2 + b1*d_tc*x + (b0+b1*tc)*d_x + c1*d_tc*x2 &
4751 + (c0+c1*tc)*d_x2) + pt*pt*e*d_x2
4752 end function gpscompressibility_P0_2
4753
4754end module calcHeightAndPressure_mod
