Ticket #443: ropp_apps_utils_ROPP90.f90

File ropp_apps_utils_ROPP90.f90, 30.8 KB (added by Ian Culverwell, 10 years ago)

ropp_apps_utils.f90

Line 
1! $Id: ropp_apps_utils.f90 4707 2016-02-11 16:51:04Z idculv $
2
3MODULE ropp_apps_utils
4
5!****m* Modules/ropp_apps_utils *
6!
7! NAME
8! ropp_apps_utils - Utility routines for ROPP applications module
9!
10! SYNOPSIS
11! USE ropp_apps_utils
12!
13! DESCRIPTION
14! This module provides signal and matrix processing routines used by
15! the ROPP pre-processor ionospheric correction package.
16!
17! AUTHOR
18! Met Office, Exeter, UK and ECMWF, Reading, UK.
19! Any comments on this software should be given via the ROM SAF
20! Helpdesk at http://www.romsaf.org
21!
22! COPYRIGHT
23! (c) EUMETSAT. All rights reserved.
24! For further details please refer to the file COPYRIGHT
25! which you should have received as part of this distribution.
26!
27!****
28
29CONTAINS
30
31! ==============================================================================
32
33!****s* APPSUtils/ropp_apps_calc_geop *
34!
35! NAME
36! ropp_apps_calc_geop - calculate geopotential for ECMWF-like bgr fields
37! from {T, q, Ak, Bk}
38!
39! SYNOPSIS
40! CALL ropp_apps_calc_geop(ro_data, geop)
41!
42!****
43
44 SUBROUTINE ropp_apps_calc_geop(ro_data, geop)
45!
46! Calculate geopotential for ECMWF-like bgr files from {T, q, Ak, Bk}.
47! Slightly modified code lifted from relevant section of
48! ropp_fm/common/ropp_fm_roprof2state.f90,
49! to which the reader is referred for explanatory comments.
50
51 USE typesizes, ONLY: wp => EightByteReal
52 USE ropp_utils
53 USE ropp_io_types, ONLY: ROprof
54 USE ropp_apps_constants, ONLY: R_dry, g_wmo
55! USE ropp_apps_utils, not_this => ropp_apps_calc_geop
56
57 IMPLICIT NONE
58
59! In/out
60 TYPE(ROprof), INTENT(IN) :: ro_data ! RO profile containing lev2b data
61 REAL(wp), DIMENSION(:), INTENT(INOUT) :: geop
62
63! Local variables
64 REAL(wp), DIMENSION(:), ALLOCATABLE :: p_hlv, geop_hlv
65 REAL(wp), DIMENSION(:), ALLOCATABLE :: Tvflv, ln_prflv
66 REAL(wp), DIMENSION(:), ALLOCATABLE :: del_p, del_geop, alpha
67
68 INTEGER :: n_hlv, n_flv, lvl
69
70 CHARACTER(len=256) :: level_type
71 CHARACTER(len=256) :: routine
72
73! 1.1 Initialisation
74! ------------------
75
76 CALL message_get_routine(routine)
77
78 CALL message_set_routine('ropp_apps_calc_geop')
79
80! 1.2 Set up
81! ----------
82
83 level_type = ro_data%Lev2d%level_type
84 IF ( INDEX(level_type,'UNKNOWN') > 0) level_type = ro_data%bg%source
85
86 CALL To_Upper(level_type)
87
88 IF ( (INDEX(level_type,'HYBRID') > 0) .OR. &
89 (INDEX(level_type,'ECMWF') > 0) ) THEN ! ECMWF background
90
91 n_hlv = SIZE(ro_data%lev2d%level_coeff_a)
92 n_flv = SIZE(ro_data%lev2b%press)
93
94 IF ( n_hlv /= n_flv+1 ) &
95 CALL message(msg_fatal, "Inconsistent number of pressure levels " // &
96 "and Ak coefficients \n")
97
98 ALLOCATE (p_hlv(n_hlv), geop_hlv(n_hlv))
99 ALLOCATE (Tvflv(n_flv), ln_prflv(n_flv))
100 ALLOCATE (del_p(n_flv), del_geop(n_flv), alpha(n_flv))
101
102 p_hlv = ro_data%lev2d%level_coeff_a + &
103 ro_data%lev2d%level_coeff_b * ro_data%lev2c%press_sfc
104
105 WHERE (ro_data%lev2b%shum > ropp_MDTV) ! very occasionally they're not
106 Tvflv = (1.0_wp + 0.61e-3_wp * ro_data%lev2b%shum) * ro_data%lev2b%temp !NB: shum in g/kg
107 ELSEWHERE
108 Tvflv = (1.0_wp + 0.61e-3_wp * 0.0_wp) * ro_data%lev2b%temp !NB: shum in g/kg
109 END WHERE
110
111 del_p(1:n_flv) = p_hlv(1:n_hlv-1) - p_hlv(2:n_hlv)
112
113 ln_prflv(1:n_flv-1) = LOG(p_hlv(1:n_hlv-2)/p_hlv(2:n_hlv-1))
114 ln_prflv(n_flv) = 0.0_wp ! For completeness (never used)
115
116 alpha(1:n_flv-1) = 1.0_wp - p_hlv(2:n_hlv-1)/del_p(1:n_flv-1) * ln_prflv(1:n_flv-1)
117 alpha(n_flv) = LOG(2.0_wp)
118
119 del_geop = R_dry * Tvflv * ln_prflv / g_wmo
120
121! 1.3 Integration
122! ---------------
123
124 geop_hlv(1) = ro_data%lev2c%geop_sfc
125 DO lvl = 2, n_hlv-1
126 geop_hlv(lvl) = geop_hlv(lvl-1) + del_geop(lvl-1)
127 END DO
128
129! 1.4 Interpolation
130! -----------------
131
132 geop = geop_hlv(1:n_hlv-1) + alpha * R_dry * Tvflv / g_wmo
133
134! 1.5 Cleaning up
135! ---------------
136
137 DEALLOCATE (del_p, del_geop, alpha)
138 DEALLOCATE (Tvflv, ln_prflv)
139 DEALLOCATE (p_hlv, geop_hlv)
140
141 ELSE
142
143! 1.6 Computer says no
144! --------------------
145
146 CALL message(msg_info, "Background profile not in ECMWF format ... " // &
147 "cannot calculate geopotentials \n")
148 geop = ropp_MDFV
149
150 END IF
151
152! 1.7 Clear up
153! ------------
154
155 CALL message_set_routine(routine)
156
157 END SUBROUTINE ropp_apps_calc_geop
158
159
160! ==============================================================================
161
162!****s* APPSUtils/ropp_apps_calc_tdry *
163!
164! NAME
165! ropp_apps_calc_tdry - calculate dry temperature from refractivity profile
166!
167! SYNOPSIS
168! CALL ropp_apps_calc_tdry(ro_data, tdry)
169!
170!****
171
172 SUBROUTINE ropp_apps_calc_tdry(ro_data, tdry)
173!
174! Calculate geopotential dry temperature from refractivity profile.
175! We integrate dp/dz = (-g(z)/k_1 R) N using a very crude climatological
176! estimate of dT/dz at the top, based on z(top).
177! Euler forward, so 1st order accurate. Probably good enough here - we just
178! need something with the same gradients as Tdry in the lower troposphere.
179
180 USE typesizes, ONLY: wp => EightByteReal
181 USE ropp_utils
182 USE ropp_io, ONLY: ropp_io_ascend
183 USE ropp_io_types, ONLY: ROprof
184 USE ropp_apps_constants, ONLY: R_dry, g_wmo, kappa1
185
186 IMPLICIT NONE
187
188! In/out
189 TYPE(ROprof), INTENT(IN) :: ro_data ! RO profile containing lev2a data
190 REAL(wp), DIMENSION(:), INTENT(INOUT) :: tdry
191
192! Local variables
193 REAL(wp), DIMENSION(:), ALLOCATABLE :: z, ref, t, lnp
194 REAL(wp) :: const ! g_wmo/R/kappa1
195 REAL(wp) :: p_top, z_top, dlognbydz_top, dtbydz_top
196
197! 'Climatological' lapse rate parameters (actually dT/dz = - usual lapse rate)
198 REAL(wp), PARAMETER :: z_surf=0.0_wp ! gpm
199 REAL(wp), PARAMETER :: z_trop_lower=11.0E3_wp ! gpm
200 REAL(wp), PARAMETER :: z_trop_midpt=12.0E3_wp ! gpm
201 REAL(wp), PARAMETER :: z_trop_upper=13.0E3_wp ! gpm
202 REAL(wp), PARAMETER :: z_strat_lower=47.5E3_wp ! gpm
203 REAL(wp), PARAMETER :: z_strat_midpt=50.0E3_wp ! gpm
204 REAL(wp), PARAMETER :: z_strat_upper=52.5E3_wp ! gpm
205 REAL(wp), PARAMETER :: z_mes_lower=87.5E3_wp ! gpm
206 REAL(wp), PARAMETER :: z_mes_midpt=90.0E3_wp ! gpm
207 REAL(wp), PARAMETER :: z_mes_upper=92.5E3_wp ! gpm
208
209 REAL(wp), PARAMETER :: lrt_trop=-8.5E-3_wp ! K/m
210 REAL(wp), PARAMETER :: lrt_strat=1.7E-3_wp ! K/m
211 REAL(wp), PARAMETER :: lrt_mes=-2.5E-3_wp ! K/m
212 REAL(wp), PARAMETER :: lrt_therm=4.8E-3_wp ! K/m
213
214 INTEGER :: i, nlev
215
216 CHARACTER(len=256) :: routine
217 CHARACTER(len=10) :: s_z_top, s_dtbydz_top
218
219 LOGICAL, DIMENSION(:), ALLOCATABLE :: lvalid
220
221! 1.1 Initialisation
222! ------------------
223
224 CALL message_get_routine(routine)
225
226 CALL message_set_routine('ropp_apps_calc_tdry')
227
228! 1.2 Set up
229! ----------
230
231 ALLOCATE ( lvalid(ro_data%lev2a%npoints) )
232
233 lvalid = .FALSE.
234
235 WHERE ( (ro_data%lev2a%refrac > ropp_ZERO) .AND. &
236 (ro_data%lev2a%alt_refrac > ropp_MDTV) ) lvalid = .TRUE.
237
238 nlev = COUNT(lvalid)
239
240 IF ( nlev < 3 ) THEN
241
242 CALL message(msg_info, "Not enough valid refractivities and heights ... " // &
243 "cannot calculate dry temperatures \n")
244
245 tdry = ropp_MDFV
246
247 CALL message_set_routine(routine)
248
249 RETURN
250
251 END IF
252
253! 1.3 Make sure heights are increasing
254! ------------------------------------
255
256 CALL ropp_io_ascend(ro_data)
257
258! 1.4 Define constants and variables
259! ----------------------------------
260
261 ALLOCATE ( z(nlev), ref(nlev), t(nlev), lnp(nlev) )
262
263 z = PACK(ro_data%lev2a%alt_refrac, lvalid)
264
265 z = geometric2geopotential(ro_data%georef%lat, z) ! Because we're assuming a constant g_wmo
266
267 ref = PACK(ro_data%lev2a%refrac, lvalid)
268
269 const = g_wmo / r_dry / kappa1
270
271! 1.5 Boundary condition on p
272! ---------------------------
273
274! Estimate dlogN/dz(top=(nlev-1)).
275
276 dlognbydz_top = LOG(ref(nlev)/ref(nlev-2)) / (z(nlev) - z(nlev-2))
277
278 IF ( dlognbydz_top > -ropp_ZDTV ) THEN
279
280 CALL message(msg_warn, "dlogN/dz non-negative at top of profile ... " // &
281 "cannot calculate dry temperature.\n")
282
283 tdry = ropp_MDFV
284
285 CALL message_set_routine(routine)
286
287 RETURN
288
289 END IF
290
291! Set dT/dz(top). Do not rearrange the order of these statements!
292
293 z_top = z(nlev-1)
294
295 IF (z_top > z_surf) dtbydz_top = lrt_trop
296 IF (z_top > z_trop_lower) dtbydz_top = 0.0_wp
297 IF (z_top > z_trop_upper) dtbydz_top = lrt_strat
298 IF (z_top > z_strat_lower) dtbydz_top = 0.0_wp
299 IF (z_top > z_strat_upper) dtbydz_top = lrt_mes
300 IF (z_top > z_mes_lower) dtbydz_top = 0.0_wp
301 IF (z_top > z_mes_upper) dtbydz_top = lrt_therm
302
303 p_top = -ref(nlev-1) * ((g_wmo/R_dry) + dtbydz_top) / &
304 kappa1 / dlognbydz_top
305
306 WRITE (s_z_top, FMT='(F10.5)') z_top / 1.E3_wp
307 WRITE (s_dtbydz_top, FMT='(F10.5)') dtbydz_top / 1.E-3_wp
308
309 CALL message(msg_diag, "Setting dT/dz = " // s_dtbydz_top // " K/km," // &
310 " at z_top = " // s_z_top // " km.\n")
311
312 lnp(nlev-1) = LOG(p_top)
313
314! 1.6 Integrate downwards from z(nlev-1) to z(1)
315! ----------------------------------------------
316
317 DO i=nlev-1,2,-1
318
319! lnp(i-1) = lnp(i) - const * ref(i) * (z(i-1) - z(i)) * EXP(-lnp(i))
320
321 lnp(i-1) = lnp(i) - 0.5_wp * const * ref(i) * (z(i-1) - z(i)) * EXP(-lnp(i)) ! Estimate of lnp(i-1/2)
322
323 lnp(i-1) = lnp(i) - const * SQRT(ref(i)*ref(i-1)) * (z(i-1) - z(i)) * EXP(-lnp(i-1)) ! Using previous estimate of lnp(i-1/2) on RHS
324
325 END DO
326
327! 1.7 Recalculate lnp at top
328! --------------------------
329
330 lnp(nlev) = lnp(nlev-2) - const * ref(nlev-1) * (z(nlev) - z(nlev-2)) * EXP(-lnp(nlev-1))
331
332! 1.8 Convert p to t
333! ------------------
334
335 t = kappa1 * EXP(lnp) / ref
336
337! 1.9 Unpack t to return tdry
338! ---------------------------
339
340 tdry = UNPACK(t, lvalid, ropp_MDFV)
341
342! 2.0 Clear up
343! ------------
344
345 DEALLOCATE (lnp, t, ref, z, lvalid)
346
347 CALL message_set_routine(routine)
348
349 END SUBROUTINE ropp_apps_calc_tdry
350
351! ==============================================================================
352
353!****s* APPSUtils/ropp_apps_impact2geom *
354!
355! NAME
356! ropp_apps_impact2geom - calculate geometric height from impact parameter
357!
358! SYNOPSIS
359! CALL ropp_apps_impact2geom(ro_data, geom, pblh_qc_flag)
360!
361!****
362
363 SUBROUTINE ropp_apps_impact2geom(ro_data, geom, pblh_qc_flag)
364!
365! Calculate geometric height h corresponding to a given impact parameter a
366! by solving (iteratively) a = r n(r), where
367! n = 1 + refrac*10^-6 and
368! r = h + RoC + und
369! Refractivity on (lev1b) r values is found by interpolating (lev2a)
370! logs of the refracs.
371
372 USE typesizes, ONLY: wp => EightByteReal
373 USE ropp_utils
374 USE ropp_io
375 USE ropp_io_types, ONLY: ROprof
376 USE ropp_apps_types
377! USE ropp_apps_utils, not_this => ropp_apps_impact2geom
378
379 IMPLICIT NONE
380
381! In/out
382 TYPE(ROprof), INTENT(IN) :: ro_data ! RO profile containing lev2b data
383 REAL(wp), DIMENSION(:), INTENT(INOUT) :: geom
384 INTEGER, INTENT(INOUT) :: pblh_qc_flag
385
386! Local parameters
387 REAL(wp), PARAMETER :: refrac_surf_default=300.0_wp
388 REAL(wp), PARAMETER :: tol=1.0e-3_wp ! m
389 INTEGER, PARAMETER :: iter_max=250
390
391! Local variables
392 REAL(wp) :: roc, und
393 REAL(wp) :: refrac_surf
394 REAL(wp), DIMENSION(MAX(1, ro_data%lev1b%npoints)) &
395 :: old_rad, new_rad, refrac
396
397 INTEGER, DIMENSION(:), POINTER :: idx => NULL() ! Array indices
398 INTEGER :: nidx
399
400 INTEGER :: iter
401 CHARACTER(LEN=3) :: siter_max
402 CHARACTER(len=256) :: routine
403
404! 1.1 Initialisation
405! ------------------
406
407 CALL message_get_routine(routine)
408
409 CALL message_set_routine('ropp_apps_impact2geom')
410
411! 1.2 Definition
412! --------------
413
414 IF ( ro_data%lev1b%npoints < 1 ) THEN
415 CALL message(msg_fatal, "No level 1b data in profile \n")
416 CALL message_set_routine(routine)
417 RETURN
418 END IF
419
420 IF ( ro_data%GEOref%roc < ropp_MDTV ) THEN
421 CALL message(msg_error, "Missing radius of curvature ... " // &
422 "cannot calculate geometric height \n")
423 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
424 CALL message_set_routine(routine)
425 RETURN
426 ELSE
427 roc = ro_data%GEOref%roc
428 END IF
429
430 IF ( ro_data%GEOref%undulation < ropp_MDTV ) THEN
431 CALL message(msg_warn, "Missing undulation ... assuming it to be zero \n")
432 und = 0.0_wp
433 ELSE
434 und = ro_data%GEOref%undulation
435 END IF
436
437 idx => WHERE (ro_data%lev2a%refrac > ropp_ZERO, nidx)
438 IF (nidx == 0) THEN
439 CALL message(msg_error, "No positive refractivities ... " // &
440 "cannot calculate heights from impact parameters \n")
441 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
442 CALL message_set_routine(routine)
443 RETURN
444 END IF
445
446 IF ( (ro_data%lev2a%refrac(idx(1)) < 10.0_wp) .OR. &
447 (ro_data%lev2a%refrac(idx(1)) > 1000.0_wp) ) THEN
448 CALL message(msg_info, "First refractivity < 10 N-units or > 1000 N-units ... " // &
449 "using 300 N-units on first iteration \n")
450 refrac_surf = refrac_surf_default
451 ELSE
452 refrac_surf = ro_data%lev2a%refrac(idx(1))
453 END IF
454
455! 1.3 Iteration
456! -------------
457
458 iter = 0
459
460 old_rad = ro_data%lev1b%impact / (1.0_wp + 1.0e-6_wp*refrac_surf)
461
462 DO
463
464! Check monotonicity
465
466 IF ( ANY( (old_rad(2:ro_data%lev1b%npoints ) - &
467 old_rad(1:ro_data%lev1b%npoints-1)) < 0) ) THEN
468 CALL message(msg_error, 'Calculation generated non-monotonic refractivity altitudes')
469 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
470 EXIT
471 END IF
472
473 refrac = EXP( ropp_apps_interpol_1d(ro_data%lev2a%alt_refrac(idx), &
474 LOG(ro_data%lev2a%refrac(idx)), &
475 old_rad - roc - und, .TRUE.) )
476
477 new_rad = ro_data%lev1b%impact / (1.0_wp + 1.0e-6_wp*refrac)
478
479 IF ( MAXVAL( ABS( old_rad - new_rad ) ) <= tol ) EXIT
480
481 iter = iter + 1
482
483 IF ( iter > iter_max ) THEN
484 WRITE ( siter_max, '(i3)' ) iter_max
485 CALL message(msg_error, 'Iteration failed to converge within ' // &
486 siter_max // ' iterations')
487 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
488 EXIT
489 END IF
490
491 old_rad = new_rad
492
493 END DO
494
495! 1.4 Calculation
496! ---------------
497
498 IF ( .NOT. BTEST(pblh_qc_flag, PBLH_QC_data_invalid) ) &
499 geom = new_rad - roc - und
500
501! 1.5 Clear up
502! ------------
503
504 CALL message_set_routine(routine)
505
506
507 END SUBROUTINE ropp_apps_impact2geom
508
509! ==============================================================================
510
511!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_locate *
512!
513! NAME
514! ropp_apps_pblh_locate
515!
516! SYNOPSIS
517! Locate planetary boundary layer height diagnostic by interpolating dlapse/dz
518!
519! CALL ropp_apps_pblh_locate(pblh_index, var_x, var_y, var_grad, pblh, pblx)
520!
521! DESCRIPTION
522! Hone PBLH estimate by find zero of lapse rate in vicinity of max/min.
523!
524! INPUTS
525! REAL(wp), DIMENSION(:), INTENT(IN) :: var_x ! Variable
526! REAL(wp), DIMENSION(:), INTENT(IN) :: var_y ! Height
527! REAL(wp), DIMENSION(:), INTENT(IN) :: var_grad ! Gradient
528! INTEGER, INTENT(IN) :: pblh_index ! Location of max/min
529!
530! OUTPUTS
531! REAL(wp), INTENT(OUT) :: pblh ! Height of PBL
532! REAL(wp), INTENT(OUT) :: pblx ! Variable at PBLH
533! INTEGER, INTENT(INOUT) :: pblh_qc_flag ! QC flag
534!
535! AUTHOR
536! Met Office, Exeter, UK.
537! Any comments on this software should be given via the ROM SAF
538! Helpdesk at http://www.romsaf.org
539!
540! COPYRIGHT
541! (c) EUMETSAT. All rights reserved.
542! For further details please refer to the file COPYRIGHT
543! which you should have received as part of this distribution.
544!
545!****
546
547 SUBROUTINE ropp_apps_pblh_locate(pblh_index, var_x, var_y, var_grad, pblh, pblx, pblh_qc_flag)
548
549 USE typesizes, ONLY: wp => EightByteReal
550 USE ropp_utils
551 USE ropp_io
552 USE ropp_apps_types
553
554 IMPLICIT NONE
555
556! I/O
557 REAL(wp), INTENT(OUT) :: pblh ! Height of PBL
558 REAL(wp), INTENT(OUT) :: pblx ! Variable at PBLH
559
560 REAL(wp), DIMENSION(:), INTENT(IN) :: var_x
561 REAL(wp), DIMENSION(:), INTENT(IN) :: var_y
562 REAL(wp), DIMENSION(:), INTENT(IN) :: var_grad
563
564 INTEGER, INTENT(IN) :: pblh_index
565 INTEGER, INTENT(INOUT) :: pblh_qc_flag
566
567! Local
568 REAL(wp) :: dgamma_by_dz_m,dgamma_by_dz_p ! Lapse rates
569 REAL(wp) :: alpha, beta, dy_p, dy_m
570 INTEGER :: n_points1
571 CHARACTER(LEN=256) :: routine
572
573! 1.1 Initialisation
574! ------------------
575
576 CALL message_get_routine(routine)
577
578 CALL message_set_routine('ropp_apps_pblh_locate')
579
580! 1.2 Estimate location of PBLH by interpolating dlapse/dz
581! --------------------------------------------------------
582
583 n_points1 = SIZE(var_x)
584
585 IF (.NOT. BTEST(pblh_qc_flag, PBLH_QC_data_invalid)) THEN
586
587 IF ( (pblh_index <= 2) .OR. (pblh_index >= (n_points1-1)) ) THEN
588
589 CALL message(msg_diag, "Maximum gradient at top or bottom of profile \n")
590
591 pblh = (var_y(MAX(pblh_index, 1)) + var_y(MIN(pblh_index+1, n_points1))) / 2.0_wp
592
593 pblx = var_x(MAX(pblh_index, 1)) + & ! using var_grad between pblh_index and pblh_index+1
594 var_grad(MAX(pblh_index, 1)) * (pblh - var_y(MAX(pblh_index, 1)))
595
596 ELSE
597
598 CALL message(msg_diag, "Calculating PBLH by interpolating lapse rate gradient")
599
600 dy_m = 0.5_wp * (var_y(pblh_index+1) - var_y(pblh_index-1)) ! at i=pblh_index
601 dy_p = 0.5_wp * (var_y(pblh_index+2) - var_y(pblh_index )) ! at i=pblh_index+1
602
603 dgamma_by_dz_m = (var_grad(pblh_index ) - var_grad(pblh_index-1)) / dy_m ! at i=pblh_index
604
605 dgamma_by_dz_p = (var_grad(pblh_index+1) - var_grad(pblh_index )) / dy_p ! at i=pblh_index+1
606
607 IF ( (dgamma_by_dz_m*dgamma_by_dz_p > ropp_ZERO) .OR. &
608 (ABS(dgamma_by_dz_m-dgamma_by_dz_p) < ropp_ZDTV) ) THEN
609
610 CALL message(msg_info, "Could not find PBLH by interpolating ... " // &
611 "leaving as missing data \n")
612 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_data_invalid)
613
614 ELSE
615
616! Fit gamma = gamma(i*) + alpha*(y-y(i*)) + beta*(y-y(i*))^2 through (i*-1, i*, i*+1) where i*=pblh_index
617
618 alpha = (dgamma_by_dz_p*dy_m + dgamma_by_dz_m*dy_p) / (dy_m + dy_p)
619
620 beta = (dgamma_by_dz_p - dgamma_by_dz_m ) / (dy_m + dy_p)
621
622 pblh = 0.5_wp * (var_y(pblh_index) + var_y(pblh_index+1)) - &
623 (0.5_wp*alpha/beta)
624
625 pblx = 0.5_wp * (var_x(pblh_index) + var_x(pblh_index+1)) + &
626 (0.5_wp*alpha/beta)*(-var_grad(pblh_index) + (alpha**2/6.0_wp/beta))
627
628 END IF
629
630 END IF
631
632 END IF
633
634! 1.3 Clear up
635! ------------
636
637 CALL message_set_routine(routine)
638
639 END SUBROUTINE ropp_apps_pblh_locate
640
641! ==============================================================================
642
643!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_check_min_height *
644!
645! NAME
646! ropp_apps_pblh_check_min_height
647!
648! SYNOPSIS
649! Check that the derived PBLH(s) is/are not less than pblh_min.
650!
651! CALL ropp_apps_pblh_check_min_height(n_pblh, pblh_min, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)
652!
653! DESCRIPTION
654! Check that (up to) 2 PBLHs are >= pblh_min
655!
656! INPUTS
657! REAL(wp), INTENT(IN) :: pblh_min ! Minimum allowed PBLH
658! INTEGER, INTENT(IN) :: n_pblh ! No. of PBLHs
659!
660! REAL(wp), INTENT(INOUT) :: pblh1, pblh2 ! Heights of PBL
661! REAL(wp), INTENT(INOUT) :: pblx1, pblx2 ! Variable at PBLH
662! INTEGER, INTENT(INOUT) :: pblh_qc_flag ! PBLH QC flag
663!
664! OUTPUTS
665! REAL(wp), INTENT(INOUT) :: pblh1, pblh2 ! Heights of PBL
666! REAL(wp), INTENT(INOUT) :: pblx1, pblx2 ! Variable at PBLH
667! INTEGER, INTENT(INOUT) :: pblh_qc_flag ! PBLH QC flag
668!
669! AUTHOR
670! Met Office, Exeter, UK.
671! Any comments on this software should be given via the ROM SAF
672! Helpdesk at http://www.romsaf.org
673!
674! COPYRIGHT
675! (c) EUMETSAT. All rights reserved.
676! For further details please refer to the file COPYRIGHT
677! which you should have received as part of this distribution.
678!
679!****
680
681 SUBROUTINE ropp_apps_pblh_check_min_height(n_pblh, pblh_min, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)
682
683 USE typesizes, ONLY: wp => EightByteReal
684 USE ropp_utils
685 USE ropp_io
686 USE ropp_apps_types
687
688 IMPLICIT NONE
689
690! I/O
691 REAL(wp), INTENT(IN) :: pblh_min ! Minimum allowed PBLH
692 INTEGER, INTENT(IN) :: n_pblh ! No. of PBLHs
693
694 REAL(wp), INTENT(INOUT) :: pblh1, pblh2 ! Heights of PBL
695 REAL(wp), INTENT(INOUT) :: pblx1, pblx2 ! Variable at PBLH
696 INTEGER, INTENT(INOUT) :: pblh_qc_flag ! PBLH QC flag
697
698! Local
699 INTEGER :: n_pblh_temp ! No. of PBLHs
700 CHARACTER(LEN=7) :: str_pblh,str_pblh_min
701 CHARACTER(LEN=256) :: routine
702
703! 1.1 Initialisation
704! ------------------
705
706 CALL message_get_routine(routine)
707
708 CALL message_set_routine('ropp_apps_pblh_check_min_height')
709
710 CALL message(msg_diag, "Checking that the PBLH(s) is/are not too low")
711
712! 1.2 Check PBLH(s) >= pblh_min
713! -----------------------------
714
715 n_pblh_temp = n_pblh
716
717 IF (n_pblh > 1) THEN
718
719 IF ((pblh1 < pblh_min) .AND. (pblh1 > ropp_MDTV)) THEN
720 WRITE(str_pblh , FMT='(F7.3)') 1.0e-3_wp*pblh1
721 WRITE(str_pblh_min, FMT='(F7.3)') 1.0e-3_wp*pblh_min
722 CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
723 " is below minimum of " // str_pblh_min // " km")
724 pblh1 = ropp_MDFV
725 pblx1 = ropp_MDFV
726 n_pblh_temp = n_pblh_temp - 1
727 END IF
728
729 IF ((pblh2 < pblh_min) .AND. (pblh2 > ropp_MDTV)) THEN
730 WRITE(str_pblh , FMT='(F7.3)') 1.0e-3_wp*pblh2
731 WRITE(str_pblh_min, FMT='(F7.3)') 1.0e-3_wp*pblh_min
732 CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
733 " is below minimum of " // str_pblh_min // " km")
734 pblh2 = ropp_MDFV
735 pblx2 = ropp_MDFV
736 n_pblh_temp = n_pblh_temp - 1
737 END IF
738
739 IF ((pblh1 < ropp_MDTV) .AND. (pblh2 < ropp_MDTV)) THEN
740 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_low)
741 ELSE
742 IF (pblh1 < ropp_MDTV) THEN ! swap them around
743 pblh1 = pblh2 ; pblh2 = ropp_MDFV
744 pblx1 = pblx2 ; pblx2 = ropp_MDFV
745 END IF
746 END IF
747
748 ELSE ! n_pblh = 1
749
750 IF ((pblh1 < pblh_min) .AND. (pblh1 > ropp_MDTV)) THEN
751 WRITE(str_pblh , FMT='(F7.3)') 1.0e-3_wp*pblh1
752 WRITE(str_pblh_min, FMT='(F7.3)') 1.0e-3_wp*pblh_min
753 CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
754 " is below minimum of " // str_pblh_min // " km")
755 pblh1 = ropp_MDFV
756 pblx1 = ropp_MDFV
757 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_low)
758 n_pblh_temp = n_pblh_temp - 1
759 END IF
760
761 END IF
762
763! 1.3 Reset PBLH_QC_flag if necessary
764! -----------------------------------
765
766 SELECT CASE (n_pblh_temp)
767
768 CASE(:1)
769 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
770 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)
771
772 CASE(2)
773 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_double_pblh)
774 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)
775
776 CASE(3:)
777 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
778 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_multiple_pblh)
779
780 END SELECT
781
782! 1.4 Clear up
783! ------------
784
785 CALL message_set_routine(routine)
786
787 END SUBROUTINE ropp_apps_pblh_check_min_height
788
789! ==============================================================================
790
791!****s* PlanetaryBoundaryLayerHeight/ropp_apps_pblh_check_max_height *
792!
793! NAME
794! ropp_apps_pblh_check_max_height
795!
796! SYNOPSIS
797! Check that the derived PBLH(s) is/are not greater than pblh_max.
798!
799! CALL ropp_apps_pblh_check_max_height(n_pblh, pblh_max, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)
800!
801! DESCRIPTION
802! Check that (up to) 2 PBLHs are <= pblh_max
803!
804! INPUTS
805! REAL(wp), INTENT(IN) :: pblh_max ! Maximum allowed PBLH
806! INTEGER, INTENT(IN) :: n_pblh ! No. of PBLHs
807!
808! REAL(wp), INTENT(INOUT) :: pblh1, pblh2 ! Heights of PBL
809! REAL(wp), INTENT(INOUT) :: pblx1, pblx2 ! Variable at PBLH
810! INTEGER, INTENT(INOUT) :: pblh_qc_flag ! PBLH QC flag
811!
812! OUTPUTS
813! REAL(wp), INTENT(INOUT) :: pblh1, pblh2 ! Heights of PBL
814! REAL(wp), INTENT(INOUT) :: pblx1, pblx2 ! Variable at PBLH
815! INTEGER, INTENT(INOUT) :: pblh_qc_flag ! PBLH QC flag
816!
817! AUTHOR
818! Met Office, Exeter, UK.
819! Any comments on this software should be given via the ROM SAF
820! Helpdesk at http://www.romsaf.org
821!
822! COPYRIGHT
823! (c) EUMETSAT. All rights reserved.
824! For further details please refer to the file COPYRIGHT
825! which you should have received as part of this distribution.
826!
827!****
828
829 SUBROUTINE ropp_apps_pblh_check_max_height(n_pblh, pblh_max, pblh1, pblx1, pblh2, pblx2, pblh_qc_flag)
830
831 USE typesizes, ONLY: wp => EightByteReal
832 USE ropp_utils
833 USE ropp_io
834 USE ropp_apps_types
835
836 IMPLICIT NONE
837
838! I/O
839 REAL(wp), INTENT(IN) :: pblh_max ! Maximum allowed PBLH
840 INTEGER, INTENT(IN) :: n_pblh ! No. of PBLHs
841
842 REAL(wp), INTENT(INOUT) :: pblh1, pblh2 ! Heights of PBL
843 REAL(wp), INTENT(INOUT) :: pblx1, pblx2 ! Variable at PBLH
844 INTEGER, INTENT(INOUT) :: pblh_qc_flag ! PBLH QC flag
845
846! Local
847 INTEGER :: n_pblh_temp ! No. of PBLHs
848 CHARACTER(LEN=7) :: str_pblh,str_pblh_max
849 CHARACTER(LEN=256) :: routine
850
851! 1.1 Initialisation
852! ------------------
853
854 CALL message_get_routine(routine)
855
856 CALL message_set_routine('ropp_apps_pblh_check_max_height')
857
858 CALL message(msg_diag, "Checking that the PBLH(s) is/are not too high")
859
860! 1.2 Check PBLH(s) <= pblh_max
861! -----------------------------
862
863 n_pblh_temp = n_pblh
864
865 IF (n_pblh > 1) THEN
866
867 IF ((pblh1 > pblh_max) .AND. (pblh1 > ropp_MDTV)) THEN
868 WRITE(str_pblh , FMT='(F7.3)') 1.0e-3_wp*pblh1
869 WRITE(str_pblh_max, FMT='(F7.3)') 1.0e-3_wp*pblh_max
870 CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
871 " is above maximum of " // str_pblh_max // " km")
872 pblh1 = ropp_MDFV
873 pblx1 = ropp_MDFV
874 n_pblh_temp = n_pblh_temp - 1
875 END IF
876
877 IF ((pblh2 > pblh_max) .AND. (pblh2 > ropp_MDTV)) THEN
878 WRITE(str_pblh , FMT='(F7.3)') 1.0e-3_wp*pblh2
879 WRITE(str_pblh_max, FMT='(F7.3)') 1.0e-3_wp*pblh_max
880 CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
881 " is above maximum of " // str_pblh_max // " km")
882 pblh2 = ropp_MDFV
883 pblx2 = ropp_MDFV
884 n_pblh_temp = n_pblh_temp - 1
885 END IF
886
887 IF ((pblh1 < ropp_MDTV) .AND. (pblh2 < ropp_MDTV)) THEN
888 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_high)
889 ELSE
890 IF (pblh1 < ropp_MDTV) THEN ! swap them around
891 pblh1 = pblh2 ; pblh2 = ropp_MDFV
892 pblx1 = pblx2 ; pblx2 = ropp_MDFV
893 END IF
894 END IF
895
896 ELSE ! n_pblh = 1
897
898 IF ((pblh1 > pblh_max) .AND. (pblh1 > ropp_MDTV)) THEN
899 WRITE(str_pblh , FMT='(F7.3)') 1.0e-3_wp*pblh1
900 WRITE(str_pblh_max, FMT='(F7.3)') 1.0e-3_wp*pblh_max
901 CALL message(msg_info, "Derived PBLH of " // str_pblh // " km" // &
902 " is above maximum of " // str_pblh_max // " km")
903 pblh1 = ropp_MDFV
904 pblx1 = ropp_MDFV
905 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_too_high)
906 n_pblh_temp = n_pblh_temp - 1
907 END IF
908
909 END IF
910
911! 1.3 Reset PBLH_QC_flag if necessary
912! -----------------------------------
913
914 SELECT CASE (n_pblh_temp)
915
916 CASE(:1)
917 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
918 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)
919
920 CASE(2)
921 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_double_pblh)
922 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_multiple_pblh)
923
924 CASE(3:)
925 pblh_qc_flag = IBCLR(pblh_qc_flag, PBLH_QC_double_pblh)
926 pblh_qc_flag = IBSET(pblh_qc_flag, PBLH_QC_multiple_pblh)
927
928 END SELECT
929
930! 1.4 Clear up
931! ------------
932
933 CALL message_set_routine(routine)
934
935 END SUBROUTINE ropp_apps_pblh_check_max_height
936
937! ==============================================================================
938
939!****f* APPSUtils/ropp_apps_interpol_1d *
940!
941! NAME
942! ropp_apps_interpol_1d - interpolate in 1D
943!
944! SYNOPSIS
945! y_out = ropp_apps_interpol_1d(x_in, y_in, x_out, flat)
946!
947!****
948 FUNCTION ropp_apps_interpol_1d(x_in, y_in, x_out, flat) RESULT (y_out)
949
950! Linearly interpolates 1d arrays, with constant or linear extrapolation.
951! Assumes x_in and x_out are monotonically increasing.
952
953 USE typesizes, ONLY: wp => EightByteReal
954! USE ropp_apps_utils, not_this => ropp_apps_interpol
955
956 REAL(wp), DIMENSION(:), INTENT(IN) :: x_in, y_in, x_out
957 REAL(wp), DIMENSION(SIZE(x_out)) :: y_out
958 LOGICAL, OPTIONAL :: flat ! If you want flat extrapolation
959
960 REAL(wp) :: grad
961 INTEGER :: i_in, i_out, n_in, n_out
962 LOGICAL :: l_flat
963
964 n_in = SIZE(x_in)
965 n_out = SIZE(x_out)
966
967 l_flat = .FALSE.
968 IF ( PRESENT(flat) ) l_flat = flat
969
970 DO i_out=1,n_out
971
972 grad = 0.0_wp
973
974 IF ( x_out(i_out) <= x_in(1) ) THEN
975 i_in = 1
976 IF ( .NOT. l_flat ) grad = (y_in(i_in+1) - y_in(i_in)) / (x_in(i_in+1) - x_in(i_in))
977 ELSE IF ( x_out(i_out) >= x_in(n_in) ) THEN
978 i_in = n_in
979 IF ( .NOT. l_flat ) grad = (y_in(i_in) - y_in(i_in-1)) / (x_in(i_in) - x_in(i_in-1))
980 ELSE
981 i_in = SUM ( MINLOC ( x_out(i_out)-x_in, MASK=(x_out(i_out)-x_in >= 0) ) )
982 grad = (y_in(i_in+1) - y_in(i_in)) / (x_in(i_in+1) - x_in(i_in))
983 END IF
984
985 y_out(i_out) = y_in(i_in) + grad * (x_out(i_out) - x_in(i_in)) ! Linear (inter/extra)polation
986
987 END DO
988
989 END FUNCTION ropp_apps_interpol_1d
990
991
992END MODULE ropp_apps_utils