1module presProfileOperators_mod
2 ! MODULE presProfileOperators_mod (prefix='ppo' category='8. Low-level utilities and constants')
3 !
4 !:Purpose: Vertical interpolation, integration, and layer averaging subroutines.
5 ! Includes the special routines designed to interpolate to the
6 ! (widely spaced) RTTOV pressure levels.
7 !
8 use utilities_mod
9
10 implicit none
11 save
12 private
13
14 ! public procedures
15
16 ! Linear interpolation routine
17 public :: ppo_lintv
18
19 ! Stand-alone interpolator calling routine providing interpolation weights
20 ! For weights W(:,:) and initial vector X(:), the interpolated vector would be W*X
21 ! Will call one of the following prodedures:
22 ! ppo_layeravgInterpWgts - weights for linear piecewise weighted averaging interpolation
23 ! ppo_piecewiseLinearWgts - weights for piecewise linear interpolation
24 ! The first of routines is similar in content to those used with RTTOV pressure levels.
25 public :: ppo_vertInterpWgts ! Main call routine
26
27 ! Routines common to ppo_vertIntegWgts and ppo_vertAvgWgts indicated below.
28 ! ppo_getLevelIndex - get the vertical input level index for level
29 ! within target layer and nearest specified layer boundary
30 public :: ppo_vertLayersSetup
31
32 ! Work arrays for ppo_vertInteg* and ppo_vertAvg*
33 real(8), allocatable :: boundaries(:), weights(:,:)
34
35 ! Stand-alone integration routines providing integration weights
36 ! For weights W(:,:) and initial vector X(:), the integrated values would be W*X
37 public :: ppo_vertIntegWgts
38
39 ! Stand-alone layer averaging routine (in lp(P)) providing weights
40 ! For weights W(:,:) and initial vector X(:), the average values would be W*X
41 public :: ppo_vertAvgWgts
42
43 contains
44
45 !--------------------------------------------------------------------------
46 ! PPO_LINTV
47 !--------------------------------------------------------------------------
48 SUBROUTINE PPO_LINTV (PVLEV,PVI,KNI, KNPROF,KNO,PPO,PVO)
49 !
50 !:Purpose: To perform the vertical interpolation in log of pressure and
51 ! constant-value extrapolation of one-dimensional vectors. Input
52 ! pressure levels can vary for each profile.
53 !
54 implicit none
55
56 ! Arguments:
57 real(8) ,intent(in) :: pvlev(kni,knprof) ! Vertical levels, pressure (source)
58 real(8) ,intent(in) :: pvi(kni,knprof) ! Vector to be interpolated (source)
59 integer ,intent(in) :: kni ! Number of input levels (source)
60 integer ,intent(in) :: knprof ! Number of profiles
61 integer ,intent(in) :: kno ! Number of output levels (destination)
62 real(8) ,intent(in) :: ppo(kno,knprof) ! Vertical levels, pressure (destination)
63 real(8) ,intent(out) :: pvo(kno,knprof) ! Interpolated profiles (destination)
64
65 ! Locals:
66 INTEGER JI, JK, JO, profileIndex, IK, IORDER
67 REAL(8) ZPI (0:KNI+1,KNPROF)
68 REAL(8) ZPVI(0:KNI+1,KNPROF)
69 INTEGER IL (KNO ,KNPROF)
70 REAL(8) ZW1, ZW2
71 REAL(8) ZP, XI, ZRT, ZP1, ZP2
72
73 !
74 !** 1. Initialization for vertical extrapolation (extra dummy levels)
75 ! . --------------------------------------------------------------
76 !
77
78 ZPI(0,:)=2000.D0
79 ZPI(KNI+1,:)=2000.D0
80
81 !
82 !** 1.1 Determine if input pressure levels are in ascending or
83 ! . descending order.
84 ! . -------------------------------------------------------
85 !
86 IF ( PVLEV(1,1) < PVLEV(KNI,1) ) THEN
87 IORDER = 1
88 ELSE
89 IORDER = -1
90 END IF
91 !
92 !** 2. Compute pressure levels pressure
93 ! . ------------------------------------------------
94 !
95
96 !
97 !** 2.1 Source levels
98 ! . -------------
99 !
100
101 ZPI(1:KNI,:) = PVLEV(1:KNI,:)
102
103 !
104 !
105 !* 3. Interpolate in log of pressure or extrapolate with constant value
106 !* . for each destination pressure level
107 ! . -----------------------------------------------------------------
108 !
109
110 !
111 !
112 !* . 3.1 Find the adjacent level below
113 ! . -----------------------------
114 !
115 !
116
117 IL(:,:)=0
118 !
119 DO JI=1,KNI
120 DO profileIndex = 1, KNPROF
121 DO JO=1,KNO
122 ZRT = PPO(JO,profileIndex)
123 ZP = ZPI(JI,profileIndex)
124 XI = SIGN(1.0D0,IORDER*(ZRT-ZP))
125 IL(JO,profileIndex) = IL(JO,profileIndex) + MAX(0.0D0,XI)
126 END DO
127 END DO
128 END DO
129 !
130 !
131 !* . 3.2 Fill extra levels, for constant value extrapolation
132 ! . ---------------------------------------------------
133 !
134
135 DO profileIndex = 1, KNPROF
136 DO JK = 1, KNI
137 ZPVI(JK,profileIndex) = PVI(JK,profileIndex)
138 END DO
139 END DO
140 DO profileIndex = 1, KNPROF
141 ZPVI(0 ,profileIndex) = PVI(1 ,profileIndex)
142 ZPVI(KNI+1,profileIndex) = PVI(KNI,profileIndex)
143 END DO
144 !
145 !
146 !* . 3.3 Interpolation/extrapolation
147 ! . ---------------------------
148 !
149
150 DO profileIndex = 1, KNPROF
151 DO JO=1,KNO
152 ZP = PPO(JO,profileIndex)
153 IK = IL(JO,profileIndex)
154 ZP1 = ZPI(IK ,profileIndex)
155 ZP2 = ZPI(IK+1,profileIndex)
156 ZW1 = LOG(ZP/ZP2)/LOG(ZP1/ZP2)
157 ZW2 = 1.D0 - ZW1
158 PVO(JO,profileIndex) = ZW1*ZPVI(IK,profileIndex) + ZW2*ZPVI(IK+1,profileIndex)
159 END DO
160 END DO
161
162 END SUBROUTINE PPO_LINTV
163
164 !==========================================================================
165 !---- Stand-alone interpolation routine providing interpolation weights ---
166 ! For weights W(:,:) and initial vector X(:), the integrated values would be W*X.
167
168 !--------------------------------------------------------------------------
169 ! ppo_vertInterpWgts
170 !--------------------------------------------------------------------------
171 subroutine ppo_vertInterpWgts(pressInput,pressTarget,numInputLevs,numTargetLevs, &
172 wgts,kstart,kend,method_opt,skipType_opt,outbound_opt,success_opt)
173 !
174 !:Purpose: Determination of interpolation weights for interpolation to points
175 ! in a profile. Applies interpolation in log(Pressure).
176 !
177 !:Input:
178 ! :pressInput: pressure on reference column levels assumed to be in ascending order
179 ! :pressTarget: target pressure levels assumed to be in ascending order
180 ! :numInputLevs: number of input/reference levels
181 ! :numTargetLevs: number of target levels
182 ! :method_opt: Specified interpolation method
183 ! :skipType_opt: Skipping processing of specific target levels depending on case:
184 ! 'default' - extrapolation allowed and skipping application via input success_opt only
185 ! 'noExtrap' - no extrapolation as well as additional skipping via input success_opt
186 ! 'doAll&noExtrap' - no extrapolation only (all other levels processed)
187 !:InOut:
188 ! :outbound_opt: Flag set when beyond range of reference levels
189 ! :success_opt: LOgical indicating a valid target level
190 !
191 !:Output:
192 !
193 ! :wgts: Interpolation coefficients/weights
194 ! :kstart: Index of first relevant referenence/input level for each target level
195 ! :kend: Index of last relevant referenence/input level for each target levbel
196 !
197 implicit none
198
199 ! Arguments:
200 integer, intent(in) :: numInputLevs ! # of original/input vertical levels
201 integer, intent(in) :: numTargetLevs ! # of target vertical levels
202 real(8), intent(in) :: pressInput(numInputLevs) ! pressure on reference column levels assumed in ascending order
203 real(8), intent(in) :: pressTarget(numTargetLevs) ! target pressure levels assumed to be in ascending order
204 character(len=*), optional, intent(in) :: method_opt ! Specified interpolation method
205 character(len=*), optional, intent(in) :: skipType_opt ! Skipping processing of specific target levels depending on case
206 integer, optional, intent(inout) :: outbound_opt(numTargetLevs) ! Flag indicating if obs outside model vertical range (0 for no)
207 logical, optional, intent(inout) :: success_opt(numTargetLevs) ! success of interpolation
208 real(8), intent(out) :: wgts(numTargetLevs,numInputLevs) ! Averaging weights
209 integer, intent(out) :: kstart(numTargetLevs) ! Index of first relevant original/input level for each target level
210 integer, intent(out) :: kend(numTargetLevs) ! Index of last relevant original/input level for each target level
211
212 ! Locals:
213 integer :: TargetIndex
214 real(8) :: logPressInput(numInputLevs),logPressTarget(numTargetLevs)
215 real(8) :: pieceLinearWgts(numTargetLevs,2)
216 logical :: success(numTargetLevs)
217 character(len=20) :: method,skipType
218
219 if (present(method_opt)) then
220 method = method_opt
221 else
222 method = 'default'
223 end if
224
225 if (present(skipType_opt)) then
226 skipType = skipType_opt
227 else
228 skipType = 'default'
229 end if
230
231 if (present(success_opt)) then
232 if ( trim(skipType) == 'doAll&noExtrap' ) then
233 success(:)=.true.
234 else
235 success(:) = success_opt(:)
236 end if
237 else
238 success(:) = .true.
239 end if
240
241 ! Flag to prevent extrapolation beyond the range of the reference column levels
242
243 if ( trim(skipType) == 'noExtrap' .or. trim(skipType) == 'doAll&noExtrap' ) then
244
245 do TargetIndex=1,numTargetLevs
246
247 ! Check if target level is outside reference column range
248 if ( PressTarget(TargetIndex) < PressInput(1) .or. &
249 PressTarget(TargetIndex) > PressInput(numInputLevs) ) then
250 success(TargetIndex)=.false.
251 success_opt(TargetIndex)=.false.
252 if (.not.present(outbound_opt)) cycle
253 if (PressTarget(TargetIndex) < PressInput(1)) then
254 outbound_opt(TargetIndex)=1
255 else
256 outbound_opt(TargetIndex)=2
257 end if
258 end if
259 end do
260 end if
261
262 if ( trim(method) == 'default' ) then
263
264 ! Piecewise log-linear interpolation
265
266 logPressInput(:) = log(PressInput(:))
267 logPressTarget(:) = log(PressTarget(:))
268
269 call ppo_piecewiseLinearWgts(logPressTarget,logPressInput,numTargetLevs,numInputLevs, &
270 pieceLinearWgts,kstart,success)
271
272 kend(:) = kstart(:) + 1
273 do TargetIndex=1,numTargetLevs
274 if (success(TargetIndex)) then
275 kend(TargetIndex) = kstart(TargetIndex) + 1
276 wgts(TargetIndex,kstart(TargetIndex))=pieceLinearWgts(TargetIndex,1)
277 wgts(TargetIndex,kend(TargetIndex))=pieceLinearWgts(TargetIndex,2)
278 else
279 kstart(TargetIndex) = 1
280 kend(TargetIndex) = 1
281 wgts(TargetIndex,1)=0.0d0
282 end if
283 end do
284
285 else if ( trim(method) == 'wgtAvg' .and. numTargetLevs > 1) then
286
287 ! Piecewise weighted averaging according to distance
288 ! Involves all model levels within the profile range
289
290 logPressInput(:) = log(PressInput(:))
291 logPressTarget(:) = log(PressTarget(:))
292
293 call ppo_layeravgInterpWgts(logPressTarget,logPressInput,numTargetLevs,numInputLevs, &
294 wgts,kstart,kend,validLevel_opt=success)
295
296 else
297 call utl_abort('ppo_vertInterpWgts: This interpolation observation operator is not recognized')
298 end if
299
300 end subroutine ppo_vertInterpWgts
301
302 !------------ routines for interface ppo_linearInterpWgts ----------------
303
304 !--------------------------------------------------------------------------
305 ! ppo_piecewiseLinearWgts
306 !--------------------------------------------------------------------------
307 subroutine ppo_piecewiseLinearWgts(pvo,pvi,kno,kni,wgts,kstart,validLevel_opt)
308 !
309 !:Purpose: To obtain peacewise linear interpolation weigths.
310 ! Assumes pv*(i) < pv*(i+1). Constant values extrapolation is applied.
311 !
312 implicit none
313
314 ! Arguments:
315 real(8), intent(in) :: pvi(kni) ! Vertical levels, pressure (source)
316 integer, intent(in) :: kni ! Number of input levels (source)
317 integer, intent(in) :: kno ! Number of output levels (destination)
318 real(8), intent(in) :: pvo(kno) ! Vertical levels, pressure (destination)
319 real(8), intent(out) :: wgts(kno,2) ! Interpolation weights (destination)
320 integer, intent(out) :: kstart(kno) ! Index i of pvlev level associated to pvo(j,1); pvo(j,2) is for pvlev level i+1
321 logical, optional, intent(out) :: validLevel_opt(kno)
322
323 ! Locals:
324 integer :: pviIndex, pvoIndex, ik
325 integer :: lowerlevel(0:KNO)
326 real(8) :: zw1, zp2
327 logical :: validLevel(kno)
328
329 if (present(validLevel_opt)) then
330 validLevel(:) = validLevel_opt(:)
331 else
332 validLevel(:) = .true.
333 end if
334
335 ! Find the adjacent level below
336
337 lowerlevel(:)=1
338 lowerlevel(1:kno) = 0
339 do pvoIndex=1,kno
340 if (.not.validLevel(pvoIndex)) cycle
341 do pviIndex=max(1,lowerlevel(pvoIndex-1)),kni
342 if ( pvo(pvoIndex) < pvi(pviIndex) ) then
343 lowerlevel(pvoIndex) = pviIndex
344 exit
345 end if
346 end do
347 if (pviIndex > kni .or. lowerlevel(pvoIndex) == 0 ) then
348 lowerlevel(pvoIndex:kno) = kni+1
349 exit
350 end if
351 end do
352
353 ! Determine interpolation/extrapolation weights
354
355 !$OMP PARALLEL DO PRIVATE(pvoIndex,ik,zp2,zw1)
356 do pvoIndex=1,kno
357 if (.not.validLevel(pvoIndex)) then
358 kstart(pvoIndex)=1
359 wgts(pvoIndex,1) = 0.0d0
360 wgts(pvoIndex,2) = 0.0d0
361 cycle
362 end if
363
364 ik = lowerlevel(pvoIndex)
365 if (ik <=1 ) then
366 kstart(pvoIndex)=1
367 wgts(pvoIndex,1) = 1.0d0
368 wgts(pvoIndex,2) = 0.0d0
369 else if ( ik > kni ) then
370 kstart(pvoIndex)=kni-1
371 wgts(pvoIndex,1)=0.0d0
372 wgts(pvoIndex,2)=1.0d0
373 else
374 zp2 = pvi(ik-1)
375 zw1 = (pvo(pvoIndex)-zp2)/(pvi(ik)-zp2)
376 wgts(pvoIndex,2) =zw1
377 wgts(pvoIndex,1) = 1.0d0 - zw1
378 kstart(pvoIndex) = ik - 1
379 end if
380 end do
381 !$OMP END PARALLEL DO
382
383 end subroutine ppo_piecewiseLinearWgts
384
385 !--------------------------------------------------------------------------
386 ! ppo_layeravgInterpWgts
387 !--------------------------------------------------------------------------
388 subroutine ppo_layeravgInterpWgts(PX1,PX2,KN1,KN2,PZ,kstart,kend,PZS1_opt, &
389 PZS2_opt,PZDPS_opt,rttov_opt,validLevel_opt)
390 !
391 !:Purpose: Determine profile interpolation weights by considering integrations
392 ! over of series of segments [PX1(KI-1),PX1(KI+1)] using piecewise
393 ! linear weighting with weights of zero at KI-1 and KI+1 and
394 ! max weight at KI. KI ranges from 1 to KN1.
395 !
396 ! Can also provide gradient contributions from both
397 ! linear and non-linear components of interpolator. The non-linear
398 ! components (case PZS* is present) stem from vertical coordinate
399 ! independent variables (e.g. dependency on Ps). The gradients
400 ! contributions from the linear components are the interpolation weights.
401 !
402 ! For the interpolation model f(x) where
403 !
404 ! f(v,x) = F(v)*x
405 ! ~ F(vo)*xo + F(vo)*(x-xo) + (dF/dv)*xo*(v-vo)
406 ! = F(vo)*x + (dF/dv)*xo*(v-vo) (eqn 1)
407 !
408 !
409 ! F(vo): array of interpolation weights
410 ! = array of gradients from the linear component
411 ! (dF/dv)*xo:
412 ! array of gradients from linearized component.
413 ! (dF/dv) or (dF/dv)*(v-vo) provided when pzs* is present
414 !
415 ! Method:
416 ! - Piecewise weighted interpolation in ln(P).
417 ! - Journal reference:
418 ! Rochon, Y.J., L. Garand, D.S. Turner, and S. Polavarapu.
419 ! Jacobian mapping between coordinate systems in data assimilation,
420 ! Q. J. of the Royal Met. Soc., vol 133, 1547-1558, 2007.
421 ! (www.interscience.wiley.com) DOI:10.1002/qj.117
422 !
423 ! URL:
424 ! http://www3.interscience.wiley.com/cgi-bin/fulltext/116320452/PDFSTART
425 !
426 !:Comments:
427 !
428 ! Assumption: PX1(i)<PX1(i+1) & PX2(i)<PX2(i+1)
429 !
430 ! 1) The input profile is now extrapolated at constant value.
431 !
432 ! The impact is of most practical significance for instruments where
433 ! the weighting function peaks at or near the surface such as SSMI.
434 !
435 ! This approach increases the weights of contribution from
436 ! the lowest and highest input domain levels for output
437 ! layers intersecting these input domain boundaries. It consists
438 ! of applying contant value extrapolation by introducing
439 ! 'fake' or 'virtual' layers. For the lowest level of the input domain,
440 ! as example, this implies creating a virtual surface layer which
441 ! extends to the lower boundary of the output domain layer which
442 ! contains the input domain surface. This increases the
443 ! contributing weight of the surface which would be otherwise
444 ! understimated in the original code due to the interpolator
445 ! actually doing piecewise weighted averaging.
446 !
447 ! 2) COmment out use of 'zb' for consistency with RTTOV-9 when
448 ! rttov_opt = .true.
449 ! See the four lines ending with !C1 and version 7 comment above.
450 !
451 ! 3) A major reduction in computational time results from only
452 ! assigning values to the non-zero ranges of the 2D output
453 ! arrays. These ranges of the 2D areas are identified by 'kstart'
454 ! and 'kend'. Initialization to zero for values within these ranges
455 ! is done using 1D work arrays 'zpz' and 'zpzd', with the resulting
456 ! values then being assigned to the related elements of 'PZ' and
457 ! 'PZDPS'.
458 !
459 ! Therefore, elements of 'PZ' and 'PZDPS' outside these ranges
460 ! could be undefined (i.e. NaN if not 0.0) and should not be used.
461 !
462 ! Other notable reductions in computational time stem (a) from inlining
463 ! of 'sublayer' code (applied to a reduced degree in 'layeravg' as
464 ! compared to 'rttov_layeravg_*) and (b) from updating the start
465 ! position 'istart' of the loops over J. The latter was faciliated
466 ! by moving the loop over KI inside 'layeravg'.
467 !
468 ! These improvements were originally devised and implemented by
469 ! Deborah Salmond and Mats Hamrud (ECMWF) in the RTTOV-9 routines
470 ! 'rttov_layeravg*'.
471 !
472 ! 4) Contributors to improvements and changes to the original version
473 ! of 'layeravg' and to 'rttov_layeravg*': members of the RTTOV9
474 ! development team, namely Niels Bormann, Alan Geer, Deborah Salmond,
475 ! and Mats Hamrud of ECMWF, Peter Rayer and Roger Saunders of the
476 ! Met Office and Pascal Brunel of Meteo-France, and Y.J. Rochon (EC).
477 !
478 implicit none
479
480 ! Arguments:
481 integer, intent(in) :: KN1 ! Dimension of PX1
482 integer, intent(in) :: KN2 ! Dimension of other arrays
483 real(8), intent(in) :: PX1(KN1) ! Levels of output domain (e.g. lnP; in increasing values)
484 real(8), intent(in) :: PX2(KN2) ! Levels of input domain (e.g. lnP; in increasing values)
485 real(8), intent(out) :: PZ(KN1,KN2) ! Resultant accumulated weighting factors for interp from input to output domain
486 integer, intent(out) :: KSTART(KN1) ! Start index for relevant PZ row
487 integer, intent(out) :: KEND(KN1) ! End index of relevant PZ row
488 real(8), optional, intent(in) :: PZS1_opt(KN1) ! dPX1/dv or perturbation dPX1/dv * delta(v) where PX1(v), e.g. v=Ps
489 real(8), optional, intent(in) :: PZS2_opt(KN2) ! dPX2/dv or perturbation dPX2/dv * delta(v) where PX2(v), e.g. v=Ps
490 real(8), optional, intent(out) :: PZDPS_opt(KN1,KN2) ! dF/dv or perturbations (dF/dv)*(v-vo): Resultant accumulated factors for the gradients w.r.t. v (or perturbations) associated to coordinates
491 logical, optional, intent(in) :: RTTOV_OPT ! Commented out use of 'zb' when .true. for consistency with RTTOV-9
492 logical, optional, intent(in) :: validLevel_opt(kn1) ! Logical indicating validity of each output level
493
494 ! Locals:
495 logical :: LGRADP,lgradp1,lgradp2,validLevel(kn1)
496 integer :: J,IC,ISTART,KI
497 real(8) :: Z1(0:KN1+1),ZW1,ZW2,ZSUM,ZB,ZBPS,ZWPS1,ZWPS2,ZDXD
498 real(8) :: PZ1(0:KN1+1),PZS2(KN2)
499 real(8) :: DZ(KN1+1),DZD(KN1+1),DXD(KN2+1)
500 real(8) :: ZPZ(KN2),ZPZD(KN2)
501 real(8) :: WX1,WX2,Y1,Y2,DY,DAD
502 real(8), parameter :: WGT1=0.0d0, WGT2=1.0d0, WGT3=0.0d0
503
504 !- Set integration/averaging range boundaries for output domain.
505 ! Range of integration for each layer ki is z1(ki-1) to z1(ki+1).
506 ! Weighting function is linear with weights of 0.0 at z1(ki-1) and
507 ! z1(ki+1) and a weight of 1.0 at z1(ki).
508
509 z1(1:kn1)=px1(1:kn1)
510 z1(0)=2.0*px1(1)-px1(2)
511 z1(kn1+1)=2.0*px1(kn1)-px1(kn1-1)
512
513 if (present(pzs1_opt)) then
514 lgradp1=.true.
515 lgradp=.true.
516 pz1(1:kn1)=pzs1_opt(1:kn1)
517 pz1(0)=2.0*pzs1_opt(1)-pzs1_opt(2)
518 pz1(kn1+1)=2.0*pzs1_opt(kn1)-pzs1_opt(kn1-1)
519 else
520 lgradp1=.false.
521 lgradp=.false.
522 pz1(0:kn1+1)=0.0d0
523 end if
524 if (present(pzs2_opt)) then
525 lgradp2=.true.
526 lgradp=.true.
527 pzs2(:)=pzs2_opt(:)
528 else
529 lgradp2=.false.
530 pzs2(:)=0.0d0
531 end if
532 if (present(validLevel_opt)) then
533 validLevel(:) = validLevel_opt(:)
534 else
535 validLevel(:) = .true.
536 end if
537
538 !- Pre-calculate values (dzd and dxd) used by subroutine sublayer
539
540 dz(1:kn1+1)=z1(1:kn1+1)-z1(0:kn1)
541 dzd(1:kn1+1)=1.0d0/dz(1:kn1+1)
542
543 dxd(2:kn2)=1.0d0/(px2(2:kn2)-px2(1:kn2-1))
544 dxd(1)=dxd(2)
545 dxd(kn2+1)=dxd(kn2)
546
547 !- Determine forward interpolator or TLM coefficients
548
549 !- Loop over output domain levels for determining
550 ! contributing weights of input domain levels over
551 ! segments [PX1(KI-1),PX1(KI+1)]
552
553 istart=1
554 do ki=1,kn1
555
556 if (.not.validLevel(ki)) then
557 kstart(ki) = 1
558 kend(ki) = 1
559 pz(ki,:) = 0.0d0
560 if (lgradp) pzdps_opt(ki,:)=0.0d0
561 cycle
562 end if
563
564 ! -- Consider constant value extrapolations cases for output domain
565 ! layers entirely outside the input domain.
566
567 if (z1(ki+1).le.px2(1)) then
568
569 ! pz(ki,1) is set to 1.0 to force constant value extrapolation
570 ! of field (e.g. temperature) above highest input level.
571
572 pz(ki,1)=1.0d0
573 if (lgradp) pzdps_opt(ki,1)=0.0
574 kstart(ki)=1
575 kend(ki)=1
576 CYCLE
577
578 else if (z1(ki-1).ge.px2(kn2)) then
579
580 ! pz(ki:kni,kn2) is set to 1.0 to force constant value extrapolation
581 ! of field (e.g. temperature) below lowest input level.
582
583 pz(ki:kn1,kn2)=1.0d0
584 if (lgradp) pzdps_opt(ki:kn1,kn2)=0.0
585 kstart(ki:kn1)=kn2
586 kend(ki:kn1)=kn2
587 return
588 end if
589
590 ! -- Consider piecewise averaging interpolation for output domain
591 ! layers and layer segments within the input domain.
592 !
593 ! -- Loop over input layers within the (z1(ki-1),z1(ki)) and
594 ! (z1(ki),z1(ki+1)) integration ranges.
595 !
596 ! Accumulate contributions to integration components over the
597 ! different segments.
598
599 ic=0
600 zpz(istart)=0.0d0
601 if (lgradp) then
602 zpzd(istart)=0.0d0
603 do j=istart,kn2-1
604 zpz(j+1)=0.0d0
605 zpzd(j+1)=0.0d0
606
607 if (px2(j).ge.z1(ki+1)) exit
608
609 if (px2(j).lt.z1(ki).and.px2(j+1).gt.z1(ki-1)) then
610
611 ! Integration over the segment of the range (z1(ki-1),z1(ki))
612 ! intersecting with the range (px2(j),px2(j+1))=(x1,x2)
613
614 call ppo_sublayerInterpWgts(z1(ki-1),z1(ki),dzd(ki),wgt1,wgt2, &
615 px2(j),px2(j+1),dxd(j+1),zpz(j),zpz(j+1), &
616 pz1(ki-1),pz1(ki), &
617 pzs2(j),pzs2(j+1),zpzd(j),zpzd(j+1),lgradp1,lgradp2)
618 if (ic.eq.0) ic=j
619 endif
620
621 if (px2(j+1).gt.z1(ki)) then
622
623 ! Integration over the segment of the range (z1(ki),z1(ki+1))
624 ! intersecting with the range (px2(j),px2(j+1))=(x1,x2).
625
626 call ppo_sublayerInterpWgts(z1(ki),z1(ki+1),dzd(ki+1),wgt2,wgt3, &
627 px2(j),px2(j+1),dxd(j+1),zpz(j),zpz(j+1), &
628 pz1(ki),pz1(ki+1), &
629 pzs2(j),pzs2(j+1),zpzd(j),zpzd(j+1),lgradp1,lgradp2)
630 if (ic.eq.0) ic=j
631 endif
632 enddo
633
634 if (ic.eq.0) then
635 j=kn2
636 pz(ki,j)=1.0d0
637 pzdps_opt(ki,j)=0.0
638 kstart(ki)=j
639 kend(ki)=j
640 CYCLE
641 else
642 istart=ic
643 kstart(ki)=istart
644 kend(ki)=j
645 end if
646
647 else
648
649 ! Same as above but with a compressed subset of 'sublayer' inlined
650 ! for improved speed at least in setting interplation weights.
651 ! This follows the corresponding change in 'rttov_layeravg'
652 ! by Deborah Salmond and Mats Hamrud (ECMWF) and takes advantage
653 ! of the known wgt* values in each case.
654 !
655 ! See 'ppo_sublayerInterpWgts' for information on applied equations.
656
657 do j=istart,kn2-1
658 zpz(j+1)=0.0d0
659
660 if (px2(j).ge.z1(ki+1)) exit
661
662 if (px2(j).lt.z1(ki).and.px2(j+1).gt.z1(ki-1)) then
663
664 ! Integration over the segment of the range (z1(ki-1),z1(ki))
665 ! intersecting with the range (px2(j),px2(j+1))=(x1,x2)
666
667 if (z1(ki-1).lt.px2(j)) then
668 y1=px2(j)
669 wx1=(y1-z1(ki-1))
670 else
671 y1=z1(ki-1)
672 wx1=0.0d0
673 end if
674 if (z1(ki).gt.px2(j+1)) then
675 wx2=(px2(j+1)-z1(ki-1))
676 dy=(px2(j+1)-y1)*dzd(ki)
677 zpz(j)=zpz(j)+dy*wx1
678 zpz(j+1)=zpz(j+1)+dy*wx2
679 else
680 dad=dxd(j+1)*dz(ki)*(z1(ki)-px2(j+1))
681 dy=(z1(ki)-y1)*dzd(ki)
682 zpz(j)=zpz(j)+dy*(wx1-dad)
683 zpz(j+1)=zpz(j+1)+dy*(dz(ki)+dad)
684 end if
685
686 if (ic.eq.0) ic=j
687 endif
688
689 if (px2(j+1).gt.z1(ki)) then
690
691 ! Integration over the segment of the range (z1(ki),z1(ki+1))
692 ! intersecting with the range (px2(j),px2(j+1))=(x1,x2).
693
694 if (z1(ki+1).gt.px2(j+1)) then
695 y2=px2(j+1)
696 wx2=z1(ki+1)-y2
697 else
698 y2=z1(ki+1)
699 wx2=0.0d0
700 end if
701 if (z1(ki).lt.px2(j)) then
702 wx1=z1(ki+1)-px2(j)
703 dy=(y2-px2(j))*dzd(ki+1)
704 zpz(j)=zpz(j)+dy*wx1
705 zpz(j+1)=zpz(j+1)+dy*wx2
706 else
707 dad=dxd(j+1)*dz(ki+1)*(z1(ki)-px2(j))
708 dy=(y2-z1(ki))*dzd(ki+1)
709 zpz(j)=zpz(j)+dy*(dz(ki+1)-dad)
710 zpz(j+1)=zpz(j+1)+dy*(wx2+dad)
711 end if
712
713 if (ic.eq.0) ic=j
714 endif
715 enddo
716
717 if (ic.eq.0) then
718 j=kn2
719 pz(ki,j)=1.0d0
720 kstart(ki)=j
721 kend(ki)=j
722 CYCLE
723 else
724 istart=ic
725 kstart(ki)=istart
726 kend(ki)=j
727 end if
728
729 end if
730
731 ! -- Consider constant value extrapolation contribution for output domain
732 ! layers crossing lower or upper input domain boundaries.
733 !
734 ! -- Extend to below (and/or above) lowest (and/or highest) input levels
735 ! if output layer intersects the corresponding input domain boundary.
736 !
737 ! A virtual layer covering the region between the input domain boundary
738 ! and the boundary of the output layer is added. This increases the
739 ! contributing in weight of the input domain boundary roughly by the
740 ! relative thickness of the virtual layer to the entire output layer
741 ! thickness.
742 !
743 ! This follows updates by Alan Geer (ECMWF) in rttov_layeravg
744 ! of RTTOV-9. It is done for two reasons:
745 !
746 ! 1) Piecewise weighted averaging using only the region
747 ! within the input domain will tend to underestimate the impact
748 ! of the input level boundary.
749 !
750 ! 2) The regression coefficients of the output domain
751 ! model are based on complete output domain layers, this being
752 ! consistent with RTTOV convention.
753
754 if (z1(ki+1).gt.px2(kn2).and.z1(ki-1).lt.px2(kn2)) then
755
756 ! Increase contribution from lowest input level.
757
758 zw1=0.0d0
759 zw2=0.0d0
760 zwps1=0.0d0
761 zwps2=0.0d0
762 zb=z1(ki+1)
763 zbps=pz1(ki+1)
764 if (present(rttov_opt)) then
765 if (rttov_opt) then
766 if (z1(ki+1).gt.2*px2(kn2)-px2(kn2-1).and.ki.eq.kn1) then !C1
767 zb=2*px2(kn2)-px2(kn2-1) !C1
768 zbps=2*pzs2(kn2)-pzs2(kn2-1) !C1
769 end if !C1
770 end if
771 end if
772 zdxd=1.0d0/(zb-px2(kn2))
773 if (z1(ki).lt.zb) then
774 call ppo_sublayerInterpWgts(z1(ki),z1(ki+1),dzd(ki+1),wgt2,wgt3, &
775 px2(kn2),zb,zdxd,zw1,zw2, &
776 pz1(ki),pz1(ki+1), &
777 pzs2(kn2),zbps,zwps1,zwps2,lgradp1,lgradp2)
778 end if
779 if (z1(ki).gt.px2(kn2)) then
780 call ppo_sublayerInterpWgts(z1(ki-1),z1(ki),dzd(ki),wgt1,wgt2, &
781 px2(kn2),zb,zdxd,zw1,zw2, &
782 pz1(ki-1),pz1(ki), &
783 pzs2(kn2),zbps,zwps1,zwps2,lgradp1,lgradp2)
784 end if
785 zpz(kn2)=zpz(kn2)+zw1+zw2
786 if (lgradp) zpzd(kn2)=zpzd(kn2)+zwps1+zwps2
787 end if
788 if (z1(ki-1).lt.px2(1).and.z1(ki+1).gt.px2(1)) then
789
790 ! Increase contribution from highest input level.
791
792 zw1=0.0d0
793 zw2=0.0d0
794 zwps1=0.0d0
795 zwps2=0.0d0
796 zb=z1(ki-1)
797 zbps=pz1(ki-1)
798 if (present(rttov_opt)) then
799 if (rttov_opt) then
800 if (z1(ki-1).lt.2*px2(1)-px2(2).and.ki.eq.1) then !C1
801 zb=2*px2(1)-px2(2) !C1
802 zbps=2*pzs2(1)-pzs2(2) !C1
803 end if !C1
804 end if
805 end if
806 zdxd=1.0d0/(px2(1)-zb)
807 if (z1(ki).gt.zb) then
808 call ppo_sublayerInterpWgts(z1(ki-1),z1(ki),dzd(ki),wgt1,wgt2, &
809 zb,px2(1),zdxd,zw1,zw2, &
810 pz1(ki-1),pz1(ki), &
811 zbps,pzs2(1),zwps1,zwps2,lgradp1,lgradp2)
812 end if
813 if (z1(ki).lt.px2(1)) then
814 call ppo_sublayerInterpWgts(z1(ki),z1(ki+1),dzd(ki+1),wgt2,wgt3, &
815 zb,px2(1),zdxd,zw1,zw2, &
816 pz1(ki),pz1(ki+1), &
817 zbps,pzs2(1),zwps1,zwps2,lgradp1,lgradp2)
818 end if
819 zpz(1)=zpz(1)+zw1+zw2
820 if (lgradp) zpzd(1)=zpzd(1)+zwps1+zwps2
821 end if
822
823 ! -- Normalize sum to unity (instead of calculating and dividing by
824 ! weighting denominator)
825
826 zsum=1.0d0/sum(zpz(kstart(ki):kend(ki)))
827 pz(ki,kstart(ki):kend(ki))=zpz(kstart(ki):kend(ki))*zsum
828
829 if (lgradp) then
830
831 ! Normalize and account for denominator gradients, i.e.
832
833 ! d[sum1(w*t)/sum2(w)]/dv = sum1[t*(dw/dv)]/sum2(w)
834 ! -sum1[t*w]*sum2[(dw/dv)]/sum2(w)^2
835
836 pzdps_opt(ki,kstart(ki):kend(ki))=(zpzd(kstart(ki):kend(ki)) &
837 -pz(ki,kstart(ki):kend(ki)) &
838 *sum(zpzd(kstart(ki):kend(ki))))*zsum
839 end if
840
841 end do
842
843 end subroutine ppo_layeravgInterpWgts
844
845 !--------------------------------------------------------------------------
846 ! ppo_sublayerInterpWgts
847 !--------------------------------------------------------------------------
848 subroutine ppo_sublayerInterpWgts(z1,z2,dzd,wgt1,wgt2,x1,x2,dxd,w1,w2, &
849 pzs1,pzs2,pxs1,pxs2,wps1,wps2,lgradpx,lgradpz)
850 !
851 !:Purpose: Determine weight coefficient contributions to w1 and w2 to assign
852 ! to input domain (e.g. NWP model) variables at x1 and x2. Weights
853 ! are determined from integration over the intersecting segment (y1,y2)
854 ! of the ranges (x1,x2) for the input domain and (z1,z2) for the
855 ! output domain. Integrals are approximated via the trapezoidal rule:
856 !
857 ! integral of f(x)=w(x)*t(x) from y1 to y2
858 !
859 ! = (f(y1)+f(y2))/2*(y2-y1)
860 ! = w(y1)*t(y1)+w(y2)*t(y2)
861 ! = w1*t(x1)+w2*t(x2)
862 ! = w1*t1+w2*t2
863 !
864 ! This is synonomous to having an integrand linear in x.
865 !
866 ! In the above (and below) equation(s), w1 and w2 are contributions to
867 ! the input values.
868 !
869 ! The weights for linearized contributions of non-linear interpolator
870 ! components, i.e. gradient w.r.t. the vertical coordinate
871 ! independent variable (e.g. v*=Ps), are calculated
872 ! when LGRADP* is .true.:
873 !
874 ! pzps = pzps + (df/dx1)*(dx1/dvx1)+(df/dx2)*(dx2/dvx2)
875 ! + (df/dz1)*(dz1/dvz1)+(df/dz2)*(dz2/dvz2)
876 !
877 ! = pzpz + (dw1/dx1*t1+dw2/dx1*t2)*pxs1
878 ! + (dw1/dx2*t1+dw2/dx2*t2)*pxs2
879 ! + (dw1/dz1*t1+dw2/dz1*t2)*pzs1
880 ! + (dw1/dz2*t1+dw2/dz2*t2)*pzs2
881 !
882 ! This routine provides terms on the right-hand-side.
883 !
884 ! Note: pxs* and pzs* can be provided either as gradients or
885 ! perturbations.
886 !
887 ! Method:
888 ! - Piecewise weighted interpolation in ln(P).
889 ! - Journal reference:
890 ! Rochon, Y.J., L. Garand, D.S. Turner, and S. Polavarapu.
891 ! Jacobian mapping between coordinate systems in data assimilation,
892 ! Q. J. of the Royal Met. Soc., vol 133, 1547-1558, 2007.
893 ! (www.interscience.wiley.com) DOI:10.1002/qj.117
894 !
895 ! URL:
896 ! http://www3.interscience.wiley.com/cgi-bin/fulltext/116320452/PDFSTART
897 !
898 !:Comments:
899 !
900 ! Assumptions:
901 !
902 ! 1) x1<x2
903 !
904 ! 2) z1<z2
905 !
906 ! 3) The ranges (z1,z2) and (x1,x2) overlap. The overlap region
907 ! will be identified as (y1,y2) with y1<y2.
908 !
909 ! 1) w(y1) and w(y2) are obtained by linear interpolation of the linear
910 ! weighting function w with w(z1)=wgt1 and w(z2)=wgt2.
911 !
912 ! 2) The w1 and w2 above are determined by expanding t(y1) and t(y2)
913 ! as linear functions of t(x1)=t1 and t(x2)=t2.
914 !
915 ! 3) The factor of 1/2 in
916 !
917 ! (f(y1)+f(y2))/2*(y2-y1)
918 ! = w(y1)*t(y1)+w(y2)*t(y2)
919 !
920 ! is omitted as normalization is performed in the calling routine
921 ! LAYERAVG.
922 !
923 ! 4) The code version of the interpolator part of 'int_sublayerInterpWgts'
924 ! provided for RTTOV-9 assumed the following conditions:
925 !
926 ! (wgt1,wgt2)=(0,1), d1=(y1-x1)=0 from y1=x1 or
927 ! wy1=0 from wgt1=0 and y1=z1,
928 ! or
929 !
930 ! (wgt1,wgt2)=(1,0), d2=(y2-x2)=0 when y2=x2 or
931 ! wy2=0 from wgt2=0 and y2=z2
932 !
933 ! and took account of the implications on d* and wy*.
934 !
935 ! The version presented here has each step accompanied by
936 ! related equations. It does not assume the above restrictions on
937 ! wgt1 and wgt2. This version is provided in the
938 ! comments section of the RTTOV-9 module 'rttov_sublayer'.
939 !
940 ! 5) When LGRADP*=.true., this routine provides terms needed for
941 ! the gradients w.r.t.the vertical coordinate independent variable,
942 ! e.g. Ps.
943 !
944 implicit none
945
946 ! Arguments:
947 logical, intent(in) :: LGRADPz ! Output domain logical indicating if gradient w.r.t. vertical coordinate independent variable if required i.e. d/dv where P(v) (e.g. v=Ps).
948 logical, intent(in) :: LGRADPx ! Input domain logical indicating if gradient w.r.t. vertical coordinate independent variable if required i.e. d/dv where P(v) (e.g. v=Ps).
949 real(8), intent(in) :: z1 ! Upper level of output domain half-layer (z1<z2)
950 real(8), intent(in) :: z2 ! Lower level of output domain half-layer
951 real(8), intent(in) :: x1 ! Upper boundary of input layer (x1<x2)
952 real(8), intent(in) :: x2 ! Lower boundary of input layer
953 real(8), intent(in) :: wgt1 ! Weight at z1 (0.0 or 1.0 or ...)
954 real(8), intent(in) :: wgt2 ! Weight at z2 (1.0 or 0.0 or ...)
955 real(8), intent(in) :: pzs1 ! dz1/dvz or perturbation dz1/dvz * delta(vz)
956 real(8), intent(in) :: pzs2 ! dz2/dvz or perturbation dz2/dvz * delta(vz)
957 real(8), intent(in) :: pxs1 ! dx1/dvx or perturbation dx1/dvx * delta(vx) (required when LGRADPx=.true.)
958 real(8), intent(in) :: pxs2 ! dx2/dvx or perturbation dx2/dvx * delta(vx)
959 real(8), intent(in) :: dxd ! 1.0/(x2-x1)=1.0/dx
960 real(8), intent(in) :: dzd ! 1.0/(z2-z1)=1.0/dz
961 real(8), intent(inout) :: w1 ! Starting (in) and updated (out) weight assigned to variable at upper level x1
962 real(8), intent(inout) :: w2 ! Starting (in) and updated (out) weight assigned to variable at upper level x2
963 real(8), intent(inout) :: wps1 ! Starting (in) and updated (out) value of (pxs1*dw1/dx1 + pxs2*dw1/dx2 +pzs1*dw1/dz1 + pzs2*dw1/dz2)
964 real(8), intent(inout) :: wps2 ! Starting (in) and updated (out) value of pxs1*dw2/dx1 + pxs2*dw2/dx2 +pzs1*dw2/dz1 + pzs2*dw2/dz2) (required when LGRADP*=.true.)
965
966 ! Locals:
967 real(8) :: y1 ! Upper boundary of integral range (y1<y2)
968 real(8) :: y2 ! Lower boundary of integral range
969 real(8) :: d1,d2,wx1,wx2,wy1,wy2,dy,dzdd
970 real(8) :: a1,a2,dxy1,dxy2,dxyd1,dxyd2,c1,c2
971 real(8) :: dydzdd,dy1z1,dy2z2,d1d1,d1d2
972 integer :: ibot,itop
973 real(8) :: zthreshold
974
975 ! 1. Initialization
976
977 !- Set upper and lower boundaries of integration layer (y1,y2):
978 ! (y1,y2) = ( max(x1,z1), min(x2,z2) )
979 itop=0
980 ibot=0
981 y1=z1
982 y2=z2
983 if (y1.lt.x1) then
984 y1=x1
985 itop=1
986 end if
987 if (y2.gt.x2) then
988 y2=x2
989 ibot=1
990 end if
991 dy=y2-y1
992
993 !- Verify for negative and zero (and near-zero) y2-y1 values.
994
995 zthreshold=epsilon(1.0d0)*100.0
996 c1=abs(y2)*zthreshold
997
998 if (dy.lt.-c1) then
999 write(*,*) 'y1,y2 = ',y1,y2
1000 write(*,*) 'z1,z2 = ',z1,z2
1001 write(*,*) 'x1,x2 = ',x1,x2
1002 call utl_abort("ppo_sublayerInterpWgts: dy is negative")
1003 else if (dy.lt.c1) then
1004 ! dy~0; weights can be taken as having values of zero.
1005 ! w1 and w2 unchanged.
1006 return
1007 end if
1008
1009 ! 2. Interpolation weights (equivalent to gradient contributions from
1010 ! linear component of interpolator)
1011
1012 !- Determine w(y1) and w(y2) of the integral f = w(y1)*t(y1)+w(y2)*t(y2)
1013 ! = wy1*t(y1)+wy2*t(y2)
1014
1015 ! Weights are obtained by linear interpolation of the linear
1016 ! weighting function w with w(z1)=wgt1 and w(z2)=wgt2.
1017
1018 dzdd=dzd*(wgt2-wgt1)
1019 dy1z1=y1-z1
1020 wx1=wgt1+dy1z1*dzdd
1021
1022 dydzdd=dy*dzdd
1023 wx2=wx1+dydzdd ! wx2=wgt1+(y2-z1)*dzdd=wgt2+(y2-z2)*dzdd
1024
1025 wy1=dy*wx1
1026 wy2=dy*wx2
1027
1028 !- Determine contribution of w1 and w2 for f=w1*t(x1)+w2*t(x2).
1029
1030 ! The w1 and w2 contributions above are determined by expanding t(y1)
1031 ! and t(y2) in f = w(y1)*t(y1)+w(y2)*t(y2) as linear functions of
1032 ! t(x1)=t1 and t(x2)=t2.
1033 !
1034 ! t(y1) = t1+(y1-x1)*(t2-t1)/(x2-x1) = t1+(y1-x1)*dxd*(t2-t1)
1035 ! t(y2) = t2+(y2-x2)*(t2-t1)/(x2-x1) = t2+(y2-x2)*dxd*(t2-t1)
1036 !
1037 ! Therefore,
1038 !
1039 ! f = wy1*[t1+(y1-x1)*dxd*(t2-t1)]
1040 ! +wy2*[t2+(y2-x2)*dxd*(t2-t1)]
1041 !
1042 ! = t1*[wy1-wy1*(y1-x1)*dxd-wy2*(y2-x2)*dxd]
1043 ! +t2*[wy1*(y1-x1)*dxd+wy2+wy2*(y2-x2)*dxd]
1044 !
1045 ! Aside: Contribution to w1+w2 = wy1+wy2 = 2 * w[(y2+y1)/2]
1046
1047 dxy1=y1-x1
1048 d1=dxy1*dxd
1049
1050 dxy2=y2-x2
1051 d2=dxy2*dxd
1052
1053 d1d1=1.0d0-d1
1054 d1d2=1.0d0+d2
1055
1056 w1=w1+wy1*d1d1-wy2*d2
1057 w2=w2+wy1*d1+wy2*d1d2
1058
1059 ! Aside used below: wy1*d1d1-wy2*d2 = wy1 + wy2 - (wy1*d1d1-wy2*d2)
1060
1061 ! 3. Provide gradient contributions from linearized terms of non-linear
1062 ! component of interpolator (i.e. vertical coordinate dependency)
1063
1064 if (LGRADPx) then
1065
1066 ! -- 3.1 Provide linearized contribution of non-linear interpolator
1067 ! component to TLM - this part in relation to the input vertical
1068 ! coordinate:
1069 !
1070 ! (df/dx1)*(dx1/dv)+(df/dx2)*(dx2/dv)
1071 ! = (dw1/dx1*t1+dw2/dx1*t2)*pxs1
1072 ! + (dw1/dx2*t1+dw2/dx2*t2)*pxs2
1073 !
1074 ! a1 and a2 are used below to denote dw*/dx1 and dw*/dx2.
1075 !
1076 ! Gradients of w1 and w2 with respect to x1 and x2 are done in two parts:
1077 !
1078 ! - Part I: The gradients calc excludes the cases where one
1079 ! might have y1=x1 and or y2=x2.
1080 ! - Part II: The gradient terms related to y1=x1 and or y2=x2 are
1081 ! added.
1082 !
1083 ! -- Gradient of f w.r.t x1
1084 !
1085 ! -- Gradients of w1 and w2 w.r.t. x1:
1086 !
1087 ! PART I:
1088 !
1089 ! For this case:
1090 !
1091 ! d(d1)/dx1 = dxd*(-1+dxy1*dxd)
1092 ! d(d2)/dx1 = dxd*dxy2*dxd
1093 !
1094 ! d(wy1)/dx1=0, d(wy2)/dx1=0
1095 !
1096 ! Note: d(dj)/dx* = 0.0 when dj=0.0 due to y1=x1 or y2=x2
1097
1098 dxyd1=d1*dxd
1099 dxyd2=d2*dxd
1100 a1=wy1*(dxd-dxyd1)-wy2*dxyd2
1101 a2=-a1 ! a2=wy1*(-dxd+dxyd1)+wy2*dxyd2
1102
1103 ! PART II:
1104 !
1105 ! For the extra terms: (gradients w.r.t. y1=x1)
1106 !
1107 ! d(d1)/dy1 = dxd
1108 ! d(d2)/dy1 = 0
1109 !
1110 ! d(wy1)/dy1 = -wx1+dzdd*dy
1111 ! d(wy2)/dy1 = -wx2
1112
1113 if (itop.eq.1) then
1114
1115 ! Case y1=x1 (d1=0)
1116
1117 c1=wy1*dxd
1118 c2=-wx1+dydzdd
1119 a1=a1-c1+c2+wx2*d2
1120 a2=a2+c1-wx2*d1d2
1121 end if
1122
1123 ! -- Add to accumulated gradient terms
1124
1125 wps1=wps1+a1*pxs1
1126 wps2=wps2+a2*pxs1
1127
1128 ! -- Gradient of f w.r.t. x2
1129 !
1130 ! -- Gradients of w1 and w2 w.r.t. x2:
1131 !
1132 ! PART I:
1133 !
1134 ! For this case:
1135 !
1136 ! d(d1)/dx2 = -dxd*dxy1*dxd
1137 ! d(d2)/dx2 = -dxd*(1+dxy2*dxd)
1138 !
1139 ! d(wy1)/dx2=0, d(wy2)/dx2=0
1140
1141 a1=wy1*dxyd1+wy2*(dxd+dxyd2)
1142 a2=-a1 ! a2=-wy1*dxyd1-wy2*(dxd+dxyd2)
1143
1144 ! PART II:
1145 !
1146 ! For the extra terms: (gradients w.r.t. y2=x2)
1147 !
1148 ! d(d1)/dy2 = 0
1149 ! d(d2)/dy2 = dxd
1150 !
1151 ! d(wy1)/dy2 = wx1
1152 ! d(wy2)/dy2 = wx2+dzdd*dy
1153
1154 if (ibot.eq.1) then
1155
1156 ! Case y2=x2 (d2=0)
1157
1158 c1=wy2*dxd
1159 c2=wx2+dydzdd
1160 a1=a1-c1+wx1*d1d1
1161 a2=a2+c1+c2+wx1*d1
1162 end if
1163
1164 ! -- Add to accumulated gradient terms
1165
1166 wps1=wps1+a1*pxs2
1167 wps2=wps2+a2*pxs2
1168
1169 end if
1170
1171 if (LGRADPz) then
1172
1173 ! -- 3.2 Provide linearized contribution of non-linear interpolator
1174 ! component to the TLM - this part in relation to the output vertical
1175 ! coordinate:
1176 !
1177 ! (df/dz1)*(dz1/dv)+(df/dz2)*(dz2/dv)
1178 !
1179 ! = (dw1/dz1*t1+dw2/dz1*t2)*pzs1
1180 ! + (dw1/dz2*t1+dw2/dz2*t2)*pzs2
1181 !
1182 ! a1 and a2 are used below to denote dw*/dz1 and dw*/dz2.
1183 !
1184 ! Gradients of w1 and w2 with respect to z1 and z2 are done in two parts:
1185 !
1186 ! - Part I: The gradients calc excludes the cases where one
1187 ! might have y1=z1 and or y2=z2.
1188 ! - Part II: The gradient terms related to y1=z1 and or y2=z2 are
1189 ! added.
1190 !
1191 ! -- Gradient of f w.r.t z1
1192 !
1193 ! -- Gradients of w1 and w2 w.r.t. z1:
1194 !
1195 ! PART I: Excludes y1=z1 consideration
1196 !
1197 ! d(wy1)/dz1 = dy*dzdd*((y1-z1)*dzd-1)=dy*dzdd*(y1-z2)*dzd
1198 ! d(wy2)/dz1 = dy*dzdd*(y2-z2)*dzd
1199
1200 dy2z2=y2-z2
1201 dxyd2=dydzdd*dzd
1202 dxyd1=dxyd2*dy1z1
1203 dxyd2=dxyd2*dy2z2
1204 a1=(dxyd1-dydzdd)*d1d1-dxyd2*d2
1205 a2=dxyd1+dxyd2-(dydzdd+a1)
1206
1207 if (itop.eq.0) then
1208
1209 ! PART II: y1=z1
1210 !
1211 ! d(d1)/dy1 = dxd
1212 ! d(d2)/dy1 = 0
1213 !
1214 ! d(wy1)/dy1 = -wx1+dzdd*dy
1215 ! d(wy2)/dy1 = -wx2
1216
1217 c1=wy1*dxd
1218 c2=-wx1+dydzdd
1219 a1=a1-c1+c2*d1d1+wx2*d2
1220 a2=a2+c1+c2*d1-wx2*d1d2
1221 end if
1222
1223 ! -- Add to accumulated gradient terms
1224
1225 wps1=wps1+a1*pzs1
1226 wps2=wps2+a2*pzs1
1227
1228 ! -- Gradient of f w.r.t. z2
1229 !
1230 ! -- Gradients of w1 and w2 w.r.t. z2:
1231 !
1232 ! PART I: Excludes y2=z2 consideration
1233 !
1234 ! d(wy1)/dz2 = -dy*dzdd*(y1-z1)*dzd
1235 ! d(wy2)/dz2 = -dy*dzdd*(1+(y2-z2)*dzd)=-dy*dzdd*(y2-z1)*dzd
1236
1237 a1=-dxyd1*d1d1+(dxyd2+dydzdd)*d2
1238 a2=-(dxyd1+dxyd2+dydzdd+a1)
1239
1240 if (ibot.eq.0) then
1241
1242 ! PART II: y2=z2
1243 !
1244 ! d(d1)/dy2 = 0
1245 ! d(d2)/dy2 = dxd
1246 !
1247 ! d(wy1)/dy2 = wx1
1248 ! d(wy2)/dy2 = wx2+dzdd*dy
1249
1250 c1=wy2*dxd
1251 c2=wx2+dydzdd
1252 a1=a1-c1-c2*d2+wx1*d1d1
1253 a2=a2+c1+c2*d1d2+wx1*d1
1254 end if
1255
1256 ! -- Add to accumulated gradient terms
1257
1258 wps1=wps1+a1*pzs2
1259 wps2=wps2+a2*pzs2
1260
1261 end if
1262
1263 end subroutine ppo_sublayerInterpWgts
1264
1265 !==========================================================================
1266 !--- Stand-alone integration (and averaging) routines providing weights ---
1267 ! For weights W(:,:) and initial vector X(:), the integrated (or average)
1268 ! values would be W*X.
1269
1270 !--------------------------------------------------------------------------
1271 ! ppo_vertLayersSetup
1272 !--------------------------------------------------------------------------
1273 subroutine ppo_vertLayersSetup(operatorType,pressInput,numInputLevs)
1274 !
1275 !:Purpose: Preliminary calculations for producing components required for
1276 ! vertical integration (or averaging) w.r.t. pressure for the full
1277 ! vertical rangeor a set of target layers. To be called before
1278 ! routine ppo_vertIntegWgts or ppo_vertAvgWgts.
1279 !
1280 ! Integration calculations are performed appling quadratic interpolation
1281 ! in P between level.
1282 !
1283 !:Input:
1284 ! :operatorType: 'integ' for intergration; 'avg' for averaging
1285 ! :pressInput: Reference input levels
1286 ! :numInputLevs: # of model vertical levels
1287 !
1288 !:Output:
1289 ! :boundaries(numInputLevs+1): Input layer boundaries assuming
1290 ! provided input levels can be taken as
1291 ! mid-layer values.
1292 ! :weights: Second order Lagrangian
1293 ! interp integration weights or
1294 ! unity for averaging weights
1295 !
1296 !:Comments:
1297 ! - This subroutine does the following:
1298 !
1299 ! - Setting of layer boundaries
1300 ! - If integration, determining integration weights associated
1301 ! to second order Lagrangian interpolation. Otherwise, initialize
1302 ! weights to unity.
1303 !
1304 ! - Layer boundaries are taken as mid-point between provided levels in
1305 ! lnP coordinate. Layer values are set to be the values
1306 ! interpolated to the mid-point in P within the various layers.
1307 ! Interpolation in P is done quadratically.
1308 !
1309 implicit none
1310
1311 ! Arguments:
1312 character(len=*), intent(in) :: operatorType ! 'integ' for integration; 'avg' for averaging
1313 integer, intent(in) :: numInputLevs ! # of model vertical levels
1314 real(8), intent(in) :: pressInput(numInputLevs) ! Reference input levels
1315
1316 ! Locals:
1317 integer :: levelIndex
1318 real(8) :: zp, zp1, zp2, zp3, zr1, zr2, zr3
1319
1320 ! Determine P boundaries of layers and save weights for
1321 ! use in setting integration (or averaging) weights.
1322 ! N.B.: Boundaries of layers set to mid-point between input levels
1323
1324 ! Calculate layer boundaries
1325
1326 if (allocated(boundaries)) deallocate(boundaries)
1327 allocate(boundaries(numInputLevs+1))
1328
1329 boundaries(1)=pressInput(1)
1330 boundaries(numInputLevs+1)= pressInput(numInputLevs)
1331
1332 !$OMP PARALLEL DO PRIVATE(levelIndex)
1333 DO levelIndex = 2, numInputLevs
1334 boundaries(levelIndex)=sqrt(pressInput(levelIndex-1)*pressInput(levelIndex))
1335 END DO
1336 !$OMP END PARALLEL DO
1337
1338 if (allocated(weights)) deallocate(weights)
1339
1340 if ( trim(operatorType) == 'avg' ) then
1341
1342 ! Initialize weights as scaling factors of unity.
1343
1344 allocate(weights(numInputLevs,1))
1345 weights(:,:) = 1.0d0
1346
1347 else
1348
1349 ! Set second degree Lagrangian interpolator weights
1350
1351 allocate(weights(numInputLevs,numInputLevs))
1352 weights(:,:) = 0.0d0
1353
1354 ! Interpolation to mid-layer level in P using
1355 ! second degree Lagrangian interpolator.
1356 ! N.B.: Integration is w.r.t. P
1357
1358 ! Calculating for levelIndex=1
1359
1360 zp1= pressInput(1)
1361 zp2= pressInput(2)
1362 zp3= pressInput(3)
1363 zp = (boundaries(2)+boundaries(1))/2.0
1364 zr1=(zp-zp2)*(zp-zp3)/(zp1-zp2)/(zp1-zp3)
1365 zr2=(zp-zp1)*(zp-zp3)/(zp2-zp1)/(zp2-zp3)
1366 zr3=(zp-zp2)*(zp-zp1)/(zp3-zp2)/(zp3-zp1)
1367 weights(1,1)=zr1
1368 weights(2,1)=zr2
1369 weights(3,1)=zr3
1370
1371 !$OMP PARALLEL DO PRIVATE(levelIndex,zp1,zp2,zp3,zp,zr1,zr2,zr3)
1372 DO levelIndex=2,numInputLevs-1
1373 zp1=pressInput(levelIndex-1)
1374 zp2=pressInput(levelIndex)
1375 zp3=pressInput(levelIndex+1)
1376 zp=(boundaries(levelIndex+1)+boundaries(levelIndex))/2.0
1377 zr1=(zp-zp2)*(zp-zp3)/(zp1-zp2)/(zp1-zp3)
1378 zr2=(zp-zp1)*(zp-zp3)/(zp2-zp1)/(zp2-zp3)
1379 zr3=(zp-zp2)*(zp-zp1)/(zp3-zp2)/(zp3-zp1)
1380 weights(levelIndex-1,levelIndex)=zr1
1381 weights(levelIndex,levelIndex)=zr2
1382 weights(levelIndex+1,levelIndex)=zr3
1383 ENDDO
1384 !$OMP END PARALLEL DO
1385
1386 ! Calculating for levelIndex=numInputLevs
1387
1388 zp1= pressInput(numInputLevs-2)
1389 zp2= pressInput(numInputLevs-1)
1390 zp3= pressInput(numInputLevs)
1391 zp = (boundaries(numInputLevs+1)+boundaries(numInputLevs))/2.0
1392 zr1=(zp-zp2)*(zp-zp3)/(zp1-zp2)/(zp1-zp3)
1393 zr2=(zp-zp1)*(zp-zp3)/(zp2-zp1)/(zp2-zp3)
1394 zr3=(zp-zp2)*(zp-zp1)/(zp3-zp2)/(zp3-zp1)
1395 weights(numInputLevs-2,numInputLevs)=zr1
1396 weights(numInputLevs-1,numInputLevs)=zr2
1397 weights(numInputLevs,numInputLevs)=zr3
1398
1399 end if
1400
1401 end subroutine ppo_vertLayersSetup
1402
1403 !--------------------------------------------------------------------------
1404 ! ppo_vertIntegWgts
1405 !--------------------------------------------------------------------------
1406 subroutine ppo_vertIntegWgts(targetLayersTop,targetLayersBot,numInputLevs, &
1407 numTargetLevs,kstart,kend,wgts,wgts_opt,skipType_opt, &
1408 outbound_opt,success_opt,dealloc_opt)
1409 !
1410 !:Purpose: To calculate integration weights "wgts" required for vertical integration w.r.t.
1411 ! pressure for the full vertical range or a set of target layers.
1412 ! Given the calculated weights and a user intergrand array vector X, the integral
1413 ! for a given layer i would be given by sum(wgts(i,:)*X(:))
1414 !
1415 ! Integration calculations are performed applying quadratic interpolation
1416 ! in P between level.
1417 !
1418 ! Routine ppo_vertLayersSetup to be called beforehand to generate Lagrangian weights
1419 ! and related layer boundaries (arrays 'weights' and 'boundaries')
1420 !
1421 !:Input:
1422 ! :targetLayersTop: top of target layers
1423 ! :targetLayersBot: bottom of target layers
1424 ! :numInputLevs: # of original/input vertical levels
1425 ! :numTargetLevs: # of target vertical levels
1426 ! :kstart: Index of first relevant original/input level for each target level
1427 ! :kend: Index of last relevant original/input level for each target level
1428 ! If kstart and kend are non-zero on input,
1429 ! the input are initial estimates of the values.
1430 ! :weights: See routine ppo_vertLayersSetup
1431 ! :boundaries: Boundaries of input layers
1432 ! :skipType_opt: Skipping processing of specific target layers depending on case:
1433 ! 'default' - skipping application via input success_opt only
1434 ! 'doAll&noExtrap' - application of both success_opt and outbound_opt
1435 ! :outbound_opt: Flag set when beyond range of reference/input levels
1436 ! :success_opt: Logical indicating a valid target layer
1437 ! :dealloc_opt: Logical indicating if deallocation is desired when done. (default: .false.)
1438 !
1439 !:Output:
1440 ! :wgts(numTargetLevs,numInputLevs): Integration weights
1441 ! :wgts_opt(numTargetLevs,numInputLevs): Part of integrtation weights not related to resolution
1442 !
1443 implicit none
1444
1445 ! Arguments:
1446 integer, intent(in) :: numInputLevs ! # of original/input vertical levels
1447 integer, intent(in) :: numTargetLevs ! # of target vertical levels
1448 real(8), intent(in) :: targetLayersTop(numTargetLevs) ! top of target layers
1449 real(8), intent(in) :: targetLayersBot(numTargetLevs) ! bottom of target layers
1450 real(8), intent(out) :: wgts(numTargetLevs,numInputLevs) ! Averaging weights
1451 integer, intent(inout) :: kstart(numTargetLevs) ! Index of first relevant original/input level for each target level
1452 integer, intent(inout) :: kend(numTargetLevs) ! Index of last relevant original/input level for each target level
1453 character(len=*), optional, intent(in) :: skipType_opt ! Skipping processing of specific target layers depending on case
1454 logical, optional, intent(in) :: dealloc_opt ! Logical indicating if deallocation is desired when done
1455 real(8), optional, intent(out) :: wgts_opt(numTargetLevs,numInputLevs) ! Part of averaging weights not related to resolution
1456 integer, optional, intent(inout) :: outbound_opt(numTargetLevs) ! Flag indicating if obs outside input vertical range (0 for no)
1457 logical, optional, intent(inout) :: success_opt(numTargetLevs) ! success of interpolation
1458
1459 ! Locals:
1460 integer :: TargetIndex
1461 logical :: success(numTargetLevs)
1462 character(len=20) :: skipType
1463 integer, parameter :: ivweights=2 ! Order of Lagrangian interpolation.
1464 integer :: levelIndex,JK,ILMAX2,ILMIN2
1465 integer :: ILMIN, ILMAX
1466 real(8) :: zp, zp1, zp2, zp3, zr1, zr2, zr3, ptop, pbtm
1467
1468 if (present(skipType_opt)) then
1469 skipType = skipType_opt
1470 else
1471 skipType = 'default'
1472 end if
1473
1474 if (present(success_opt)) then
1475 if ( trim(skipType) == 'doAll&noExtrap') then
1476 if (present(outbound_opt)) then
1477 where (outbound_opt(:) == 0)
1478 success(:) = .true.
1479 elsewhere
1480 success(:) = .false.
1481 end where
1482 else
1483 success(:) = .true.
1484 end if
1485 else
1486 success(:) = success_opt(:)
1487 end if
1488 else
1489 success(:) = .true.
1490 end if
1491
1492 do TargetIndex=1,numTargetLevs
1493
1494 if ( .not.success(TargetIndex) ) then
1495 wgts(TargetIndex,:) = 0.0D0
1496 wgts_opt(TargetIndex,:) = 0.0D0
1497 kstart(TargetIndex)=1
1498 kend(TargetIndex)=1
1499 cycle
1500 end if
1501
1502 ptop = targetLayersTop(TargetIndex)
1503 pbtm = targetLayersBot(TargetIndex)
1504
1505 ! Find the range of vertical levels over which to perform the integration
1506 ! and set the integration weights over this range.
1507
1508 ilmin=1
1509 ilmax=numInputLevs
1510 if (ptop <= boundaries(1)*1.01 .and. &
1511 pbtm >= boundaries(numInputLevs+1)*0.99) then
1512
1513 ! Total column integration part
1514
1515 !$OMP PARALLEL DO PRIVATE(jk,levelIndex)
1516 do jk = 1,numInputLevs
1517 do levelIndex=max(1,jk-ivweights),min(numInputLevs,jk+ivweights)
1518 wgts(TargetIndex,jk)=wgts(TargetIndex,jk) &
1519 +(boundaries(levelIndex+1) &
1520 -boundaries(levelIndex))*weights(jk,levelIndex)
1521 if (present(wgts_opt)) &
1522 wgts_opt(TargetIndex,jk)=wgts_opt(TargetIndex,jk)+ &
1523 weights(jk,levelIndex)
1524 end do
1525 end do
1526 !$OMP END PARALLEL DO
1527
1528 else
1529
1530 ! Partial column integration part (special treatment at boundaries)
1531
1532 ! Identify input layer boundaries just within the target layer.
1533
1534 ilmin = ppo_getLevelIndex(ptop, boundaries, 'top', numInputLevs+1)
1535 ilmax = ppo_getLevelIndex(pbtm, boundaries, 'btm', numInputLevs+1)
1536
1537 if (ilmin == ilmax+1) then
1538
1539 ! Entire target layer within one input layer
1540
1541 levelIndex=ilmax
1542 if (levelIndex < 3) levelIndex=3
1543 if (levelIndex > numInputLevs) levelIndex=numInputLevs
1544 zp1=boundaries(levelIndex-2)
1545 zp2=boundaries(levelIndex-1)
1546 zp3=boundaries(levelIndex)
1547 zp=(ptop+pbtm)/2.0
1548 zr1=(zp-zp2)*(zp-zp3)/(zp1-zp2)/(zp1-zp3)
1549 zr2=(zp-zp1)*(zp-zp3)/(zp2-zp1)/(zp2-zp3)
1550 zr3=(zp-zp2)*(zp-zp1)/(zp3-zp2)/(zp3-zp1)
1551
1552 wgts(TargetIndex,levelIndex-2)=(pbtm-ptop)*zr1
1553 wgts(TargetIndex,levelIndex-1)=(pbtm-ptop)*zr2
1554 wgts(TargetIndex,levelIndex)=(pbtm-ptop)*zr3
1555 if (present(wgts_opt)) then
1556 wgts_opt(TargetIndex,levelIndex-2)=zr1
1557 wgts_opt(TargetIndex,levelIndex-1)=zr2
1558 wgts_opt(TargetIndex,levelIndex)=zr3
1559 end if
1560 ilmin=levelIndex-2
1561 ilmax=levelIndex
1562
1563 else
1564
1565 ! Determine terms from the inner layers (excluding the lower and upper
1566 ! boundary layers when these layers not covering entire input layers)
1567
1568 if (pbtm >= boundaries(numInputLevs)*0.99) then
1569 ilmax2=numInputLevs
1570 else
1571 ilmax2=ilmax-1
1572 end if
1573 if (ptop <= boundaries(1)*1.01) then
1574 ilmin=1
1575 ilmin2=ilmin
1576 else
1577 ilmin2=ilmin
1578 end if
1579 if (ilmin2 <= ilmax2) then
1580 !$OMP PARALLEL DO PRIVATE(jk,levelIndex)
1581 do jk = ilmin2,ilmax2
1582 do levelIndex=max(1,jk-ivweights),min(numInputLevs,jk+ivweights)
1583 wgts(TargetIndex,jk)=wgts(TargetIndex,jk)+(boundaries(levelIndex+1) &
1584 -boundaries(levelIndex))*weights(jk,levelIndex)
1585 if (present(wgts_opt)) &
1586 wgts_opt(TargetIndex,jk)=wgts_opt(TargetIndex,jk)+weights(jk,levelIndex)
1587 end do
1588 end do
1589 !$OMP END PARALLEL DO
1590 end if
1591
1592 ! Determine terms from the lower and upper boundary layers
1593 ! when these layers do not cover entire input layers.
1594
1595 if (pbtm < boundaries(numInputLevs)*0.99) then
1596
1597 levelIndex=ilmax+1
1598 if (levelIndex > numInputLevs) levelIndex=numInputLevs
1599 if (levelIndex < 3) levelIndex=3
1600 zp1=boundaries(levelIndex-2)
1601 zp2=boundaries(levelIndex-1)
1602 zp3=boundaries(levelIndex)
1603 zp=(boundaries(ilmax)+pbtm)/2.0
1604 zr1=(zp-zp2)*(zp-zp3)/(zp1-zp2)/(zp1-zp3)
1605 zr2=(zp-zp1)*(zp-zp3)/(zp2-zp1)/(zp2-zp3)
1606 zr3=(zp-zp2)*(zp-zp1)/(zp3-zp2)/(zp3-zp1)
1607
1608 wgts(TargetIndex,levelIndex-2)=wgts(TargetIndex,levelIndex-2)+(pbtm - boundaries(ilmax))*zr1
1609 wgts(TargetIndex,levelIndex-1)=wgts(TargetIndex,levelIndex-1)+(pbtm - boundaries(ilmax))*zr2
1610 wgts(TargetIndex,levelIndex)=wgts(TargetIndex,levelIndex)+(pbtm - boundaries(ilmax))*zr3
1611
1612 if (present(wgts_opt)) then
1613 wgts_opt(TargetIndex,levelIndex-2)=wgts_opt(TargetIndex,levelIndex-2)+zr1
1614 wgts_opt(TargetIndex,levelIndex-1)=wgts_opt(TargetIndex,levelIndex-1)+zr2
1615 wgts_opt(TargetIndex,levelIndex)=wgts_opt(TargetIndex,levelIndex)+zr3
1616 end if
1617 ilmax=levelIndex
1618
1619 end if
1620
1621 if (ptop > boundaries(1)*1.01) then
1622
1623 levelIndex=ilmin-1
1624 if (levelIndex < 1) levelIndex=1
1625 if (levelIndex > numInputLevs-2) levelIndex=numInputLevs-2
1626 zp1= boundaries(levelIndex)
1627 zp2= boundaries(levelIndex+1)
1628 zp3= boundaries(levelIndex+2)
1629 zp = (boundaries(ilmin)+ptop)/2.0
1630 zr1=(zp-zp2)*(zp-zp3)/(zp1-zp2)/(zp1-zp3)
1631 zr2=(zp-zp1)*(zp-zp3)/(zp2-zp1)/(zp2-zp3)
1632 zr3=(zp-zp2)*(zp-zp1)/(zp3-zp2)/(zp3-zp1)
1633
1634 wgts(TargetIndex,levelIndex)=wgts(TargetIndex,levelIndex)+(boundaries(ilmin)-ptop)*zr1
1635 wgts(TargetIndex,levelIndex+1)=wgts(TargetIndex,levelIndex+1)+(boundaries(ilmin)-ptop)*zr2
1636 wgts(TargetIndex,levelIndex+2)=wgts(TargetIndex,levelIndex+2)+(boundaries(ilmin)-ptop)*zr3
1637
1638 if (present(wgts_opt)) then
1639 wgts_opt(TargetIndex,levelIndex)=wgts_opt(TargetIndex,levelIndex)+zr1
1640 wgts_opt(TargetIndex,levelIndex+1)=wgts_opt(TargetIndex,levelIndex+1)+zr2
1641 wgts_opt(TargetIndex,levelIndex+2)=wgts_opt(TargetIndex,levelIndex+2)+zr3
1642 end if
1643 ilmin=levelIndex
1644 if (ilmax < levelIndex+2) ilmax=levelIndex+2
1645
1646 end if
1647 if (ilmin > ilmax-2) ilmin=ilmax-2
1648 end if
1649 end if
1650
1651 if (kstart(TargetIndex) > 0 .and. kend(TargetIndex) > 0) then
1652 if (abs(kstart(TargetIndex)-ilmin) > 1 .or. &
1653 abs(kend(TargetIndex)-ilmax) > 1) then
1654
1655 write(*,*) 'ppo_vertIntegWgts: Suspected error in layer', &
1656 ' identification: ',TargetIndex,kstart(TargetIndex),ilmin, &
1657 kend(TargetIndex),ilmax
1658 end if
1659 end if
1660
1661 kstart(TargetIndex)=ilmin
1662 kend(TargetIndex)=ilmax
1663
1664 end do
1665
1666 if (present(dealloc_opt)) then
1667 if (dealloc_opt) deallocate(weights,boundaries)
1668 end if
1669
1670 end subroutine ppo_vertIntegWgts
1671
1672 !--------------------------------------------------------------------------
1673 ! ppo_getLevelIndex
1674 !--------------------------------------------------------------------------
1675 integer function ppo_getLevelIndex(level, layerBoundaryLevels, topbtm, numBoundaries)
1676 !
1677 !:Purpose: To get the vertical input level index for level
1678 ! within target layer and nearest specified layer boundary.
1679 !
1680 implicit none
1681
1682 ! Arguments:
1683 integer, intent(in) :: numBoundaries ! Number of layer boundaries
1684 real(8), intent(in) :: level ! Target layer index
1685 real(8), intent(in) :: layerBoundaryLevels(numBoundaries) ! Layer boundaries
1686 character(len=*), intent(in) :: topbtm ! indicating whether we are looking for top or bottom level
1687
1688 ! Locals
1689 integer :: ilev1, ilev2
1690 integer :: levelIndex
1691
1692 ! Find the model levels adjacent to pressure level
1693
1694 ! Default values
1695
1696 if (level < 0.0d0) then
1697 if ((topbtm == 'btm') .or. (topbtm == 'BTM')) then
1698 ppo_getLevelIndex = numBoundaries
1699 endif
1700 if ((topbtm == 'top') .or. (topbtm == 'TOP')) then
1701 ppo_getLevelIndex = 1
1702 endif
1703 endif
1704
1705 ilev1=0
1706 ilev2=1
1707 do levelIndex=1,numBoundaries
1708 if (level > layerBoundaryLevels(levelIndex)) then
1709 ilev1=levelIndex
1710 ilev2=levelIndex+1
1711 else
1712 exit
1713 endif
1714 enddo
1715
1716 ! Find the input level index
1717
1718 ! If we are looking for top level, the index is the level immediately
1719 ! below. if looking for bottom level, the index is the one immediately
1720 ! above.
1721
1722 if ((topbtm == 'btm') .or. (topbtm == 'BTM')) then
1723 ppo_getLevelIndex=ilev1
1724 else if ((topbtm == 'top') .or. (topbtm == 'TOP')) then
1725 ppo_getLevelIndex=ilev2
1726 endif
1727
1728 if (ppo_getLevelIndex < 1) ppo_getLevelIndex=1
1729 if (ppo_getLevelIndex > numBoundaries) ppo_getLevelIndex=numBoundaries
1730
1731 end function ppo_getLevelIndex
1732
1733 !--------------------------------------------------------------------------
1734 ! ppo_vertAvgWgts
1735 !--------------------------------------------------------------------------
1736 subroutine ppo_vertAvgWgts(targetLayersTop,targetLayersBot,numInputLevs, &
1737 numTargetLevs,kstart,kend,wgts,wgts_opt,skipType_opt, &
1738 outbound_opt,success_opt,dealloc_opt)
1739 !
1740 !:Purpose: To calculate averaging weights "wgts" required for vertical averaging
1741 ! w.r.t. ln(pressure) for the full vertical range or a set of target layers.
1742 ! Given the calculated weights and a user input array vector X, the average
1743 ! for a given layer i would be given by sum(wgts(i,:)*X(:))
1744 !
1745 ! Routine ppo_vertLayersSetup to be called beforehand to initial weigths
1746 ! and related layer boundaries (arrays 'weights' and 'boundaries')
1747 !
1748 !:Input:
1749 ! :targetLayersTop: top of target layers
1750 ! :targetLayersBot: bottom of target layers
1751 ! :numInputLevs: # of original/input vertical levels
1752 ! :numTargetLevs: # of target vertical levels
1753 ! :kstart: Index of first relevant original/input level for each target level
1754 ! :kend: Index of last relevant original/input level for each target level
1755 ! If kstart and kend are non-zero on input,
1756 ! the input are initial estimates of the values.
1757 !
1758 ! :weights: See routine ppo_vertLayersSetup
1759 ! :boundaries: Boundaries of input layers
1760 ! :skipType_opt: Skipping processing of specific target layers depending on case:
1761 ! 'default' - skipping application via input success_opt only
1762 ! 'doAll&noExtrap' - application of both success_opt and outbound_opt
1763 ! :outbound_opt: Flag set when beyond range of reference/input levels
1764 ! :success_opt: Logical indicating a valid target layer
1765 ! :dealloc_opt: Logical indicating if deallocation is desired when done. (default: .false.)
1766 !
1767 !:Output:
1768 ! :wgts(numTargetLevs,numInputLevs): Averaging weights
1769 ! :wgts_opt(numTargetLevs,numInputLevs): Part of averaging weights not related to resolution
1770 !
1771 implicit none
1772
1773 ! Arguments:
1774 integer, intent(in) :: numInputLevs ! # of original/input vertical levels
1775 integer, intent(in) :: numTargetLevs ! # of target vertical levels
1776 real(8), intent(in) :: targetLayersTop(numTargetLevs) ! top of target layers
1777 real(8), intent(in) :: targetLayersBot(numTargetLevs) ! bottom of target layers
1778 real(8), intent(out) :: wgts(numTargetLevs,numInputLevs) ! Averaging weights
1779 integer, intent(inout) :: kstart(numTargetLevs) ! Index of first relevant original/input level for each target level
1780 integer, intent(inout) :: kend(numTargetLevs) ! Index of last relevant original/input level for each target level
1781 character(len=*), optional, intent(in) :: skipType_opt ! Skipping processing of specific target layers depending on case
1782 logical, optional, intent(in) :: dealloc_opt ! Logical indicating if deallocation is desired when done
1783 real(8), optional, intent(out) :: wgts_opt(numTargetLevs,numInputLevs) ! Part of averaging weights not related to resolution
1784 integer, optional, intent(inout) :: outbound_opt(numTargetLevs) ! Flag indicating if obs outside input vertical range (0 for no)
1785 logical, optional, intent(inout) :: success_opt(numTargetLevs) ! success of interpolation
1786
1787 ! Locals:
1788 integer :: TargetIndex
1789 logical :: success(numTargetLevs)
1790 character(len=20) :: skipType
1791 integer :: levelIndex,ILMAX2,ILMIN2
1792 integer :: ILMIN, ILMAX
1793 real(8) :: SumWeights, TargetLayerThickWgt, ptop, pbtm
1794
1795 if (present(skipType_opt)) then
1796 skipType = skipType_opt
1797 else
1798 skipType = 'default'
1799 end if
1800
1801 if (present(success_opt)) then
1802 if ( trim(skipType) == 'doAll&noExtrap') then
1803 if (present(outbound_opt)) then
1804 where (outbound_opt(:) == 0)
1805 success(:) = .true.
1806 elsewhere
1807 success(:) = .false.
1808 end where
1809 else
1810 success(:) = .true.
1811 end if
1812 else
1813 success(:) = success_opt(:)
1814 end if
1815 else
1816 success(:) = .true.
1817 end if
1818
1819 do TargetIndex=1,numTargetLevs
1820
1821 if ( .not.success(TargetIndex) ) then
1822 wgts(TargetIndex,:) = 0.0D0
1823 wgts_opt(TargetIndex,:) = 0.0D0
1824 kstart(TargetIndex)=1
1825 kend(TargetIndex)=1
1826 cycle
1827 end if
1828
1829 ptop = targetLayersTop(TargetIndex)
1830 pbtm = targetLayersBot(TargetIndex)
1831 TargetLayerThickWgt=1.0D0/(min(pbtm,boundaries(numInputLevs+1))-max(ptop,boundaries(1)))
1832
1833 ! Find the range of vertical levels over which to perform the averaging
1834 ! and set the averaging weights over this range.
1835
1836 ilmin=1
1837 ilmax=numInputLevs
1838 if (ptop <= boundaries(1)*1.01 .and. pbtm >= boundaries(numInputLevs+1)*0.99) then
1839
1840 ! Total column averaging part
1841
1842 SumWeights=1.0D0/sum(weights(1:numInputLevs,1))
1843 !$OMP PARALLEL DO PRIVATE(levelIndex)
1844 do levelIndex = 1,numInputLevs
1845 wgts(TargetIndex,levelIndex)=(boundaries(levelIndex+1) &
1846 -boundaries(levelIndex))*TargetLayerThickWgt
1847 if (present(wgts_opt)) &
1848 wgts_opt(TargetIndex,levelIndex)=SumWeights
1849 end do
1850 !$OMP END PARALLEL DO
1851
1852 else
1853
1854 ! Partial column averaging part (special treatment at boundaries)
1855
1856 ! Identify the vertical input level indices for levels
1857 ! within target layer and nearest specified layer boundary.
1858
1859 ilmin = ppo_getLevelIndex(ptop, boundaries, 'top', numInputLevs+1)
1860 ilmax = ppo_getLevelIndex(pbtm, boundaries, 'btm', numInputLevs+1)
1861
1862 if (ilmin == ilmax+1) then
1863
1864 ! Entire target layer within one input layer
1865
1866 levelIndex=ilmin
1867 if (levelIndex < 1) levelIndex=1
1868 if (levelIndex > numInputLevs) levelIndex=numInputLevs
1869
1870 wgts(TargetIndex,levelIndex)=1.0D0
1871 if (present(wgts_opt)) wgts_opt(TargetIndex,levelIndex)=1.0D0
1872
1873 ilmin=levelIndex
1874 ilmax=levelIndex+1
1875
1876 else
1877
1878 ! Determine terms from the inner layers (excluding the lower and upper
1879 ! boundary layers when these layers not covering entire input layers)
1880
1881 if (pbtm >= boundaries(numInputLevs)*0.99) then
1882 ilmax2=numInputLevs
1883 else
1884 ilmax2=ilmax-1
1885 end if
1886 if (ptop <= boundaries(1)*1.01) then
1887 ilmin=1
1888 ilmin2=ilmin
1889 else
1890 ilmin2=ilmin
1891 end if
1892
1893 SumWeights=1.0D0/sum(weights(ilmin:ilmax,1))
1894
1895 if (ilmin2 <= ilmax2) then
1896 !$OMP PARALLEL DO PRIVATE(levelIndex)
1897 do levelIndex = ilmin2,ilmax2
1898 wgts(TargetIndex,levelIndex)= &
1899 (boundaries(levelIndex+1)-boundaries(levelIndex))*TargetLayerThickWgt
1900 if (present(wgts_opt)) wgts_opt(TargetIndex,levelIndex)=SumWeights
1901 end do
1902 !$OMP END PARALLEL DO
1903 end if
1904
1905 ! Determine terms from the lower and upper boundary layers
1906 ! when these layers do not cover entire input layers.
1907
1908 if (pbtm < boundaries(numInputLevs)*0.99) then
1909
1910 levelIndex=ilmax+1
1911 if (levelIndex > numInputLevs) levelIndex=numInputLevs
1912 if (levelIndex < 1) levelIndex=1
1913
1914 wgts(TargetIndex,levelIndex)= &
1915 (pbtm - boundaries(ilmax))*TargetLayerThickWgt
1916
1917 if (present(wgts_opt)) wgts_opt(TargetIndex,levelIndex)=SumWeights
1918
1919 ilmax=levelIndex
1920
1921 end if
1922
1923 if (ptop > boundaries(1)*1.01) then
1924
1925 levelIndex=ilmin-1
1926 if (levelIndex < 1) levelIndex=1
1927 if (levelIndex > numInputLevs) levelIndex=numInputLevs
1928
1929 wgts(TargetIndex,levelIndex)= &
1930 (boundaries(ilmin)-ptop)*TargetLayerThickWgt
1931
1932 if (present(wgts_opt)) wgts_opt(TargetIndex,levelIndex)=SumWeights
1933
1934 ilmin=levelIndex
1935 if (ilmax < levelIndex+1) ilmax=levelIndex+1
1936
1937 end if
1938 if (ilmin > ilmax-1) ilmin=ilmax-1
1939 end if
1940 end if
1941
1942 kstart(TargetIndex)=ilmin
1943 kend(TargetIndex)=ilmax
1944
1945 end do
1946
1947 if (present(dealloc_opt)) then
1948 if (dealloc_opt) deallocate(weights,boundaries)
1949 end if
1950
1951 end subroutine ppo_vertAvgWgts
1952
1953end module presProfileOperators_mod
