Line data Source code
1 : ! This module is used to diagnose the location of the tropopause. Multiple
2 : ! algorithms are provided, some of which may not be able to identify a
3 : ! tropopause in all situations. To handle these cases, an analytic
4 : ! definition and a climatology are provided that can be used to fill in
5 : ! when the original algorithm fails. The tropopause temperature and
6 : ! pressure are determined and can be output to the history file.
7 : !
8 : ! Author: Charles Bardeen
9 : ! Created: April, 2009
10 : !
11 : ! CCPP-ized: Haipeng Lin, August 2024
12 : module tropopause_find
13 : !---------------------------------------------------------------
14 : ! ... variables for the tropopause module
15 : !---------------------------------------------------------------
16 :
17 : use ccpp_kinds, only : kind_phys
18 :
19 : implicit none
20 :
21 : private
22 :
23 : ! CCPP-compliant subroutines
24 : public :: tropopause_find_init
25 : public :: tropopause_find_run
26 :
27 : ! "Wrapped" routine for use in old CAM for backward compatibility.
28 : ! Also called by tropopause_find_run driver routine.
29 : public :: tropopause_findWithBackup
30 :
31 : ! Switches for tropopause method
32 : public :: TROP_ALG_NONE, TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
33 : public :: TROP_ALG_STOBIE, TROP_ALG_HYBSTOB, TROP_ALG_TWMO, TROP_ALG_WMO
34 : public :: TROP_ALG_CPP
35 : public :: NOTFOUND
36 :
37 : save
38 :
39 : ! These parameters define an enumeration to be used to define the primary
40 : ! and backup algorithms to be used with the tropopause_find() method. The
41 : ! backup algorithm is meant to provide a solution when the primary algorithm
42 : ! fails. The algorithms that can't fail (i.e., always find a tropopause) are:
43 : ! TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE, and TROP_ALG_STOBIE.
44 : integer, parameter :: TROP_ALG_NONE = 1 ! Don't evaluate
45 : integer, parameter :: TROP_ALG_ANALYTIC = 2 ! Analytic Expression
46 : integer, parameter :: TROP_ALG_CLIMATE = 3 ! Climatology
47 : integer, parameter :: TROP_ALG_STOBIE = 4 ! Stobie Algorithm
48 : integer, parameter :: TROP_ALG_TWMO = 5 ! WMO Definition, Reichler et al. [2003]
49 : integer, parameter :: TROP_ALG_WMO = 6 ! WMO Definition
50 : integer, parameter :: TROP_ALG_HYBSTOB = 7 ! Hybrid Stobie Algorithm
51 : integer, parameter :: TROP_ALG_CPP = 8 ! Cold Point Parabolic
52 : integer, parameter :: TROP_ALG_CHEMTROP = 9 ! Chemical tropopause
53 :
54 : integer, parameter :: default_primary = TROP_ALG_TWMO ! default primary algorithm
55 : integer, parameter :: default_backup = TROP_ALG_CLIMATE ! default backup algorithm
56 :
57 : integer, parameter :: NOTFOUND = -1
58 :
59 : real(kind_phys), parameter :: ALPHA = 0.03_kind_phys
60 :
61 : ! physical constants
62 : ! These constants are set in module variables rather than as parameters
63 : ! to support the aquaplanet mode in which the constants have values determined
64 : ! by the experiment protocol
65 : real(kind_phys) :: cnst_kap ! = cappa
66 : real(kind_phys) :: cnst_faktor ! = -gravit/rair
67 : real(kind_phys) :: cnst_rga ! = 1/gravit
68 : real(kind_phys) :: cnst_ka1 ! = cnst_kap - 1._kind_phys
69 : real(kind_phys) :: cnst_rad2deg ! = 180/pi
70 :
71 : !================================================================================================
72 : contains
73 : !================================================================================================
74 :
75 : !> \section arg_table_tropopause_find_init Argument Table
76 : !! \htmlinclude tropopause_find_init.html
77 1536 : subroutine tropopause_find_init(cappa, rair, gravit, pi, errmsg, errflg)
78 :
79 : real(kind_phys), intent(in) :: cappa ! R/Cp
80 : real(kind_phys), intent(in) :: rair ! Dry air gas constant (J K-1 kg-1)
81 : real(kind_phys), intent(in) :: gravit ! Gravitational acceleration (m s-2)
82 : real(kind_phys), intent(in) :: pi ! Pi
83 :
84 : character(len=512), intent(out) :: errmsg
85 : integer, intent(out) :: errflg
86 :
87 1536 : errmsg = ' '
88 1536 : errflg = 0
89 :
90 : ! define physical constants
91 1536 : cnst_kap = cappa
92 1536 : cnst_faktor = -gravit/rair
93 1536 : cnst_rga = 1._kind_phys/gravit ! Reciprocal of gravit (s2 m-1)
94 1536 : cnst_ka1 = cnst_kap - 1._kind_phys
95 1536 : cnst_rad2deg = 180._kind_phys/pi ! radians to degrees conversion factor
96 :
97 1536 : end subroutine tropopause_find_init
98 :
99 : ! "Driver" routine for tropopause_find. Identifies the tropopause using several methods
100 : ! and populates them into the model state.
101 : ! Most methods use climatological tropopause as a backup, and as such is guaranteed to
102 : ! find a tropopause in all columns. Others are explicitly single-method and used by
103 : ! other parameterizations as-is with NOTFOUND values being intentional.
104 : !> \section arg_table_tropopause_find_run Argument Table
105 : !! \htmlinclude tropopause_find_run.html
106 0 : subroutine tropopause_find_run(ncol, pver, fillvalue, lat, pint, pmid, t, zi, zm, phis, &
107 0 : calday, tropp_p_loc, tropp_days, &
108 0 : tropLev, tropP, tropT, tropZ, & ! Default primary+backup (twmo+climate)
109 0 : tropLev_twmo, tropP_twmo, tropT_twmo, tropZ_twmo, & ! Primary only (twmo)
110 0 : tropLev_clim, tropP_clim, tropT_clim, tropZ_clim, & ! Climate-only
111 0 : tropLev_hybstob, tropP_hybstob, tropT_hybstob, tropZ_hybstob, & ! Hybridstobie + climate backup
112 0 : tropLev_cpp, tropP_cpp, tropT_cpp, tropZ_cpp, & ! Cold point only
113 0 : tropLev_chem, tropP_chem, tropT_chem, tropZ_chem, & ! Chemical tropopause only
114 0 : hstobie_trop, hstobie_linoz, hstobie_tropop, & ! Hybridstobie only for chemistry diagnostics
115 0 : scheme_name, errmsg, errflg)
116 :
117 : integer, intent(in) :: ncol ! Number of atmospheric columns
118 : integer, intent(in) :: pver ! Number of vertical levels
119 : real(kind_phys), intent(in) :: fillvalue ! Fill value for diagnostic outputs
120 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
121 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
122 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
123 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
124 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
125 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
126 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
127 :
128 : real(kind_phys), intent(in) :: calday ! Day of year including fraction of day
129 :
130 : ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
131 : real(kind_phys), intent(in) :: tropp_p_loc(:,:)
132 : real(kind_phys), intent(in) :: tropp_days(:) ! Day-of-year for climo data, 12
133 :
134 : integer, intent(out) :: tropLev(:) ! tropopause level index
135 : real(kind_phys), intent(out) :: tropP(:) ! tropopause pressure (Pa)
136 : real(kind_phys), intent(out) :: tropT(:) ! tropopause temperature (K)
137 : real(kind_phys), intent(out) :: tropZ(:) ! tropopause height (m)
138 :
139 : integer, intent(out) :: tropLev_twmo(:) ! lapse-rate tropopause level index
140 : real(kind_phys), intent(out) :: tropP_twmo(:) ! lapse-rate tropopause pressure (Pa)
141 : real(kind_phys), intent(out) :: tropT_twmo(:) ! lapse-rate tropopause temperature (K)
142 : real(kind_phys), intent(out) :: tropZ_twmo(:) ! lapse-rate tropopause height (m)
143 :
144 : integer, intent(out) :: tropLev_clim(:) ! climatology-backed tropopause level index
145 : real(kind_phys), intent(out) :: tropP_clim(:) ! climatology-backed tropopause pressure (Pa)
146 : real(kind_phys), intent(out) :: tropT_clim(:) ! climatology-backed tropopause temperature (K)
147 : real(kind_phys), intent(out) :: tropZ_clim(:) ! climatology-backed tropopause height (m)
148 :
149 : integer, intent(out) :: tropLev_hybstob(:) ! hybridstobie climatology-backed tropopause level index
150 : real(kind_phys), intent(out) :: tropP_hybstob(:) ! hybridstobie climatology-backed tropopause pressure (Pa)
151 : real(kind_phys), intent(out) :: tropT_hybstob(:) ! hybridstobie climatology-backed tropopause temperature (K)
152 : real(kind_phys), intent(out) :: tropZ_hybstob(:) ! hybridstobie climatology-backed tropopause height (m)
153 :
154 : integer, intent(out) :: tropLev_cpp(:) ! cold point tropopause level index
155 : real(kind_phys), intent(out) :: tropP_cpp(:) ! cold point tropopause pressure (Pa)
156 : real(kind_phys), intent(out) :: tropT_cpp(:) ! cold point tropopause temperature (K)
157 : real(kind_phys), intent(out) :: tropZ_cpp(:) ! cold point tropopause height (m)
158 :
159 : integer, intent(out) :: tropLev_chem(:) ! chemical tropopause level index
160 : real(kind_phys), intent(out) :: tropP_chem(:) ! chemical tropopause pressure (Pa)
161 : real(kind_phys), intent(out) :: tropT_chem(:) ! chemical tropopause temperature (K)
162 : real(kind_phys), intent(out) :: tropZ_chem(:) ! chemical tropopause height (m)
163 :
164 : ! Optional output arguments for hybridstobie with chemistry
165 : real(kind_phys), intent(out) :: hstobie_trop(:,:) ! Lowest level with strat. chem
166 : real(kind_phys), intent(out) :: hstobie_linoz(:,:) ! Lowest possible Linoz level
167 : real(kind_phys), intent(out) :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
168 :
169 : character(len=64), intent(out) :: scheme_name
170 : character(len=512), intent(out) :: errmsg
171 : integer, intent(out) :: errflg
172 :
173 0 : scheme_name = 'tropopause_find'
174 0 : errmsg = ' '
175 0 : errflg = 0
176 :
177 : ! Obtain the primary output, which is TWMO + climate
178 : call tropopause_findWithBackup( &
179 : ncol = ncol, &
180 : pver = pver, &
181 : fillvalue = fillvalue, &
182 : lat = lat, &
183 : pint = pint, &
184 : pmid = pmid, &
185 : t = t, &
186 : zi = zi, &
187 : zm = zm, &
188 : phis = phis, &
189 : calday = calday, &
190 : tropp_p_loc = tropp_p_loc, &
191 : tropp_days = tropp_days, &
192 : tropLev = tropLev, &
193 : tropP = tropP, &
194 : tropT = tropT, &
195 : tropZ = tropZ, &
196 : primary = default_primary, &
197 : backup = default_backup, &
198 : errmsg = errmsg, &
199 : errflg = errflg &
200 0 : )
201 :
202 : ! Any other intended outputs
203 : ! Primary (TWMO) only
204 : call tropopause_findWithBackup( &
205 : ncol = ncol, &
206 : pver = pver, &
207 : fillvalue = fillvalue, &
208 : lat = lat, &
209 : pint = pint, &
210 : pmid = pmid, &
211 : t = t, &
212 : zi = zi, &
213 : zm = zm, &
214 : phis = phis, &
215 : calday = calday, &
216 : tropp_p_loc = tropp_p_loc, &
217 : tropp_days = tropp_days, &
218 : tropLev = tropLev_twmo, &
219 : tropP = tropP_twmo, &
220 : tropT = tropT_twmo, &
221 : tropZ = tropZ_twmo, &
222 : primary = default_primary, &
223 : backup = TROP_ALG_NONE, &
224 : errmsg = errmsg, &
225 : errflg = errflg &
226 0 : )
227 :
228 : ! Climatology only
229 : call tropopause_findWithBackup( &
230 : ncol = ncol, &
231 : pver = pver, &
232 : fillvalue = fillvalue, &
233 : lat = lat, &
234 : pint = pint, &
235 : pmid = pmid, &
236 : t = t, &
237 : zi = zi, &
238 : zm = zm, &
239 : phis = phis, &
240 : calday = calday, &
241 : tropp_p_loc = tropp_p_loc, &
242 : tropp_days = tropp_days, &
243 : tropLev = tropLev_clim, &
244 : tropP = tropP_clim, &
245 : tropT = tropT_clim, &
246 : tropZ = tropZ_clim, &
247 : primary = TROP_ALG_CLIMATE, &
248 : backup = TROP_ALG_NONE, &
249 : errmsg = errmsg, &
250 : errflg = errflg &
251 0 : )
252 :
253 : ! Cold point (CPP) only
254 : call tropopause_findWithBackup( &
255 : ncol = ncol, &
256 : pver = pver, &
257 : fillvalue = fillvalue, &
258 : lat = lat, &
259 : pint = pint, &
260 : pmid = pmid, &
261 : t = t, &
262 : zi = zi, &
263 : zm = zm, &
264 : phis = phis, &
265 : calday = calday, &
266 : tropp_p_loc = tropp_p_loc, &
267 : tropp_days = tropp_days, &
268 : tropLev = tropLev_cpp, &
269 : tropP = tropP_cpp, &
270 : tropT = tropT_cpp, &
271 : tropZ = tropZ_cpp, &
272 : primary = TROP_ALG_CPP, &
273 : backup = TROP_ALG_NONE, &
274 : errmsg = errmsg, &
275 : errflg = errflg &
276 0 : )
277 :
278 : ! Hybridstobie with climatology-backed
279 : call tropopause_findWithBackup( &
280 : ncol = ncol, &
281 : pver = pver, &
282 : fillvalue = fillvalue, &
283 : lat = lat, &
284 : pint = pint, &
285 : pmid = pmid, &
286 : t = t, &
287 : zi = zi, &
288 : zm = zm, &
289 : phis = phis, &
290 : calday = calday, &
291 : tropp_p_loc = tropp_p_loc, &
292 : tropp_days = tropp_days, &
293 : tropLev = tropLev_hybstob, &
294 : tropP = tropP_hybstob, &
295 : tropT = tropT_hybstob, &
296 : tropZ = tropZ_hybstob, &
297 : primary = TROP_ALG_HYBSTOB, &
298 : backup = TROP_ALG_CLIMATE, &
299 : hstobie_trop = hstobie_trop, & ! Only used if TROP_ALG_HYBSTOB
300 : hstobie_linoz = hstobie_linoz, & ! Only used if TROP_ALG_HYBSTOB
301 : hstobie_tropop = hstobie_tropop, & ! Only used if TROP_ALG_HYBSTOB
302 : errmsg = errmsg, &
303 : errflg = errflg &
304 0 : )
305 :
306 : ! Chemical tropopause (used for chemistry)
307 : call tropopause_findWithBackup( &
308 : ncol = ncol, &
309 : pver = pver, &
310 : fillvalue = fillvalue, &
311 : lat = lat, &
312 : pint = pint, &
313 : pmid = pmid, &
314 : t = t, &
315 : zi = zi, &
316 : zm = zm, &
317 : phis = phis, &
318 : calday = calday, &
319 : tropp_p_loc = tropp_p_loc, &
320 : tropp_days = tropp_days, &
321 : tropLev = tropLev_chem, &
322 : tropP = tropP_chem, &
323 : tropT = tropT_chem, &
324 : tropZ = tropZ_chem, &
325 : primary = TROP_ALG_CHEMTROP, &
326 : backup = TROP_ALG_CLIMATE, &
327 : errmsg = errmsg, &
328 : errflg = errflg &
329 0 : )
330 :
331 0 : end subroutine tropopause_find_run
332 :
333 : ! Searches all the columns and attempts to identify the tropopause.
334 : ! Two routines can be specifed, a primary routine which is tried first and a
335 : ! backup routine which will be tried only if the first routine fails. If the
336 : ! tropopause can not be identified by either routine, then a NOTFOUND is returned
337 : ! for the tropopause level, temperature and pressure.
338 5235336 : subroutine tropopause_findWithBackup(ncol, pver, fillvalue, lat, pint, pmid, t, zi, zm, phis, &
339 5235336 : calday, tropp_p_loc, tropp_days, &
340 5235336 : tropLev, tropP, tropT, tropZ, &
341 5235336 : hstobie_trop, hstobie_linoz, hstobie_tropop, &
342 : primary, backup, &
343 5235336 : errmsg, errflg)
344 :
345 : integer, intent(in) :: ncol ! Number of atmospheric columns
346 : integer, intent(in) :: pver ! Number of vertical levels
347 : real(kind_phys), intent(in) :: fillvalue ! Fill value for diagnostic outputs
348 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
349 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
350 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
351 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
352 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
353 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
354 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
355 :
356 : real(kind_phys), intent(in) :: calday ! Day of year including fraction of day
357 :
358 : ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
359 : real(kind_phys), intent(in) :: tropp_p_loc(:,:)
360 : real(kind_phys), intent(in) :: tropp_days(:) ! Day-of-year for climo data, 12
361 :
362 : integer, intent(out) :: tropLev(:) ! tropopause level index
363 : real(kind_phys), optional, intent(out) :: tropP(:) ! tropopause pressure (Pa)
364 : real(kind_phys), optional, intent(out) :: tropT(:) ! tropopause temperature (K)
365 : real(kind_phys), optional, intent(out) :: tropZ(:) ! tropopause height (m)
366 :
367 : ! Optional output arguments for hybridstobie with chemistry
368 : real(kind_phys), optional, intent(inout) :: hstobie_trop(:,:) ! Lowest level with strat. chem
369 : real(kind_phys), optional, intent(inout) :: hstobie_linoz(:,:) ! Lowest possible Linoz level
370 : real(kind_phys), optional, intent(inout) :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
371 :
372 : ! primary and backup are no longer optional arguments for CCPP-compliance.
373 : ! specify defaults when calling (TWMO, CLIMO)
374 : integer, intent(in) :: primary ! primary detection algorithm
375 : integer, intent(in) :: backup ! backup detection algorithm
376 :
377 : character(len=512), intent(out) :: errmsg
378 : integer, intent(out) :: errflg
379 :
380 5235336 : errmsg = ' '
381 5235336 : errflg = 0
382 :
383 : ! Initialize the results to a missing value, so that the algorithms will
384 : ! attempt to find the tropopause for all of them.
385 87417936 : tropLev(:) = NOTFOUND
386 77013000 : if(present(tropP)) tropP(:) = fillvalue
387 77013000 : if(present(tropT)) tropT(:) = fillvalue
388 77013000 : if(present(tropZ)) tropZ(:) = fillvalue
389 :
390 : ! Try to find the tropopause using the primary algorithm.
391 5235336 : if (primary /= TROP_ALG_NONE) then
392 : call tropopause_findUsing(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
393 : calday, tropp_p_loc, tropp_days, &
394 : primary, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, &
395 : hstobie_trop=hstobie_trop, hstobie_linoz=hstobie_linoz, hstobie_tropop=hstobie_tropop, & ! only for HYBSTOB
396 7483032 : errmsg=errmsg, errflg=errflg)
397 : end if
398 :
399 87115226 : if ((backup /= TROP_ALG_NONE) .and. any(tropLev(:) == NOTFOUND)) then
400 : call tropopause_findUsing(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
401 : calday, tropp_p_loc, tropp_days, &
402 : backup, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, &
403 30349 : errmsg=errmsg, errflg=errflg)
404 : end if
405 :
406 5235336 : end subroutine tropopause_findWithBackup
407 :
408 : ! Call the appropriate tropopause detection routine based upon the algorithm
409 : ! specifed.
410 : !
411 : ! NOTE: It is assumed that the output fields have been initialized by the
412 : ! caller, and only output values set to fillvalue will be detected.
413 5250478 : subroutine tropopause_findUsing(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
414 5250478 : calday, tropp_p_loc, tropp_days, &
415 5250478 : algorithm, tropLev, tropP, tropT, tropZ, &
416 5250478 : hstobie_trop, hstobie_linoz, hstobie_tropop, &
417 5250478 : errmsg, errflg)
418 :
419 : integer, intent(in) :: ncol ! Number of atmospheric columns
420 : integer, intent(in) :: pver ! Number of vertical levels
421 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
422 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
423 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
424 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
425 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
426 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
427 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
428 :
429 : real(kind_phys), intent(in) :: calday ! Day of year including fraction of day
430 :
431 : ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
432 : real(kind_phys), intent(in) :: tropp_p_loc(:,:)
433 : real(kind_phys), intent(in) :: tropp_days(:) ! Day-of-year for climo data, 12
434 :
435 : integer, intent(in) :: algorithm ! detection algorithm
436 : integer, intent(inout) :: tropLev(:) ! tropopause level index
437 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
438 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
439 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
440 :
441 : ! Optional output arguments for hybridstobie with chemistry
442 : real(kind_phys), optional, intent(inout) :: hstobie_trop(:,:) ! Lowest level with strat. chem
443 : real(kind_phys), optional, intent(inout) :: hstobie_linoz(:,:) ! Lowest possible Linoz level
444 : real(kind_phys), optional, intent(inout) :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
445 :
446 : character(len=512), intent(out) :: errmsg
447 : integer, intent(out) :: errflg
448 :
449 5250478 : errmsg = ' '
450 5250478 : errflg = 0
451 :
452 : ! Dispatch the request to the appropriate routine.
453 5250478 : select case(algorithm)
454 : case(TROP_ALG_ANALYTIC)
455 0 : call tropopause_analytic(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
456 :
457 : case(TROP_ALG_CLIMATE)
458 : call tropopause_climate(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
459 : calday, tropp_p_loc, tropp_days, &
460 30349 : tropLev, tropP, tropT, tropZ)
461 :
462 : case(TROP_ALG_STOBIE)
463 0 : call tropopause_stobie(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
464 :
465 : case(TROP_ALG_HYBSTOB)
466 0 : if(present(hstobie_trop) .and. present(hstobie_linoz) .and. present(hstobie_tropop)) then
467 : call tropopause_hybridstobie(ncol, pver, pmid, t, zm, &
468 : tropLev, tropP, tropT, tropZ, &
469 0 : hstobie_trop, hstobie_linoz, hstobie_tropop)
470 : else
471 : call tropopause_hybridstobie(ncol, pver, pmid, t, zm, &
472 0 : tropLev, tropP, tropT, tropZ)
473 : endif
474 :
475 : case(TROP_ALG_TWMO)
476 5987664 : call tropopause_twmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
477 :
478 : case(TROP_ALG_WMO)
479 0 : call tropopause_wmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
480 :
481 : case(TROP_ALG_CPP)
482 1495368 : call tropopause_cpp(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
483 :
484 : case(TROP_ALG_CHEMTROP)
485 : ! hplin: needs climatological arguments as calling tropopause_climate from within findChemTrop
486 : call tropopause_findChemTrop(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
487 : calday, tropp_p_loc, tropp_days, &
488 : tropLev, tropP, tropT, tropZ, &
489 0 : errmsg, errflg)
490 :
491 : case default
492 0 : errflg = 1
493 5250478 : write(errmsg,*) 'tropopause: Invalid detection algorithm (', algorithm, ') specified.'
494 : end select
495 :
496 5250478 : end subroutine tropopause_findUsing
497 :
498 : ! This analytic expression closely matches the mean tropopause determined
499 : ! by the NCEP reanalysis and has been used by the radiation code.
500 0 : subroutine tropopause_analytic(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
501 0 : tropLev, tropP, tropT, tropZ)
502 :
503 : integer, intent(in) :: ncol ! Number of atmospheric columns
504 : integer, intent(in) :: pver ! Number of vertical levels
505 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
506 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
507 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
508 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
509 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
510 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
511 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
512 : integer, intent(inout) :: tropLev(:) ! tropopause level index
513 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
514 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
515 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
516 :
517 : ! Local Variables
518 : integer :: i
519 : integer :: k
520 : real(kind_phys) :: tP ! tropopause pressure (Pa)
521 :
522 : ! Iterate over all of the columns.
523 0 : do i = 1, ncol
524 :
525 : ! Skip column in which the tropopause has already been found.
526 0 : if (tropLev(i) == NOTFOUND) then
527 :
528 : ! Calculate the pressure of the tropopause.
529 0 : tP = (25000.0_kind_phys - 15000.0_kind_phys * (cos(lat(i)))**2)
530 :
531 : ! Find the level that contains the tropopause.
532 0 : do k = pver, 2, -1
533 0 : if (tP >= pint(i, k)) then
534 0 : tropLev(i) = k
535 0 : exit
536 : end if
537 : end do
538 :
539 : ! Return the optional outputs
540 0 : if (present(tropP)) tropP(i) = tP
541 :
542 0 : if (present(tropT)) then
543 0 : tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
544 : end if
545 :
546 0 : if (present(tropZ)) then
547 0 : tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
548 : end if
549 : end if
550 : end do
551 :
552 0 : end subroutine tropopause_analytic
553 :
554 : ! Determine the tropopause pressure from a climatology,
555 : ! interpolated to the current day of year and latitude.
556 30284 : subroutine tropopause_climate(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
557 30284 : calday, tropp_p_loc, tropp_days, tropLev, tropP, tropT, tropZ)
558 :
559 : integer, intent(in) :: ncol ! Number of atmospheric columns
560 : integer, intent(in) :: pver ! Number of vertical levels
561 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
562 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa), pverp
563 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
564 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
565 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
566 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
567 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
568 :
569 : real(kind_phys), intent(in) :: calday ! Day of year including fraction of day
570 :
571 : ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
572 : real(kind_phys), intent(in) :: tropp_p_loc(:,:)
573 : real(kind_phys), intent(in) :: tropp_days(:) ! Day-of-year for climo data, 12
574 :
575 : integer, intent(inout) :: tropLev(:) ! tropopause level index
576 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
577 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
578 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
579 :
580 : ! Local Variables
581 : integer :: i
582 : integer :: k
583 : integer :: m
584 : real(kind_phys) :: tP ! tropopause pressure (Pa)
585 : real(kind_phys) :: dels
586 : integer :: last
587 : integer :: next
588 :
589 : ! If any columns remain to be indentified, then get the current
590 : ! day from the calendar.
591 :
592 71839 : if (any(tropLev == NOTFOUND)) then
593 :
594 : !--------------------------------------------------------
595 : ! ... setup the time interpolation
596 : !--------------------------------------------------------
597 15142 : if( calday < tropp_days(1) ) then
598 0 : next = 1
599 0 : last = 12
600 0 : dels = (365._kind_phys + calday - tropp_days(12)) / (365._kind_phys + tropp_days(1) - tropp_days(12))
601 15142 : else if( calday >= tropp_days(12) ) then
602 0 : next = 1
603 0 : last = 12
604 0 : dels = (calday - tropp_days(12)) / (365._kind_phys + tropp_days(1) - tropp_days(12))
605 : else
606 136278 : do m = 11,1,-1
607 136278 : if( calday >= tropp_days(m) ) then
608 : exit
609 : end if
610 : end do
611 15142 : last = m
612 15142 : next = m + 1
613 15142 : dels = (calday - tropp_days(m)) / (tropp_days(m+1) - tropp_days(m))
614 : end if
615 :
616 15142 : dels = max( min( 1._kind_phys,dels ),0._kind_phys )
617 :
618 :
619 : ! Iterate over all of the columns.
620 253322 : do i = 1, ncol
621 :
622 : ! Skip column in which the tropopause has already been found.
623 253322 : if (tropLev(i) == NOTFOUND) then
624 :
625 : !--------------------------------------------------------
626 : ! ... get tropopause level from climatology
627 : !--------------------------------------------------------
628 : ! Interpolate the tropopause pressure.
629 91206 : tP = tropp_p_loc(i,last) &
630 182412 : + dels * (tropp_p_loc(i,next) - tropp_p_loc(i,last))
631 :
632 : ! Find the associated level.
633 1641214 : do k = pver, 2, -1
634 1641214 : if (tP >= pint(i, k)) then
635 91206 : tropLev(i) = k
636 91206 : exit
637 : end if
638 : end do
639 :
640 : ! Return the optional outputs
641 91206 : if (present(tropP)) tropP(i) = tP
642 :
643 91206 : if (present(tropT)) then
644 60614 : tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
645 : end if
646 :
647 91206 : if (present(tropZ)) then
648 60614 : tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
649 : end if
650 : end if
651 : end do
652 : end if
653 :
654 15142 : end subroutine tropopause_climate
655 :
656 : !-----------------------------------------------------------------------
657 : !-----------------------------------------------------------------------
658 0 : subroutine tropopause_hybridstobie(ncol, pver, pmid, t, zm, &
659 0 : tropLev, tropP, tropT, tropZ, &
660 0 : hstobie_trop, hstobie_linoz, hstobie_tropop)
661 :
662 : !-----------------------------------------------------------------------
663 : ! Originally written by Philip Cameron-Smith, LLNL
664 : !
665 : ! Stobie-Linoz hybrid: the highest altitude of
666 : ! a) Stobie algorithm, or
667 : ! b) minimum Linoz pressure.
668 : !
669 : ! NOTE: the ltrop(i) gridbox itself is assumed to be a STRATOSPHERIC gridbox.
670 : !-----------------------------------------------------------------------
671 : !-----------------------------------------------------------------------
672 : ! ... Local variables
673 : !-----------------------------------------------------------------------
674 : integer, intent(in) :: ncol ! Number of atmospheric columns
675 : integer, intent(in) :: pver ! Number of vertical levelserp
676 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
677 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
678 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
679 :
680 : integer, intent(inout) :: tropLev(:) ! tropopause level index
681 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
682 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
683 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
684 :
685 : ! Optional output arguments for hybridstobie with chemistry
686 : real(kind_phys), optional, intent(inout) :: hstobie_trop(:,:) ! Lowest level with strat. chem
687 : real(kind_phys), optional, intent(inout) :: hstobie_linoz(:,:) ! Lowest possible Linoz level
688 : real(kind_phys), optional, intent(inout) :: hstobie_tropop(:,:) ! Troposphere boundary calculated in chem.
689 :
690 : real(kind_phys),parameter :: min_Stobie_Pressure= 40.E2_kind_phys !For case 2 & 4. [Pa]
691 : real(kind_phys),parameter :: max_Linoz_Pressure =208.E2_kind_phys !For case 4. [Pa]
692 :
693 : integer :: i, k
694 : real(kind_phys) :: stobie_min, shybrid_temp !temporary variable for case 2 & 3.
695 0 : integer :: ltrop_linoz(ncol) !Lowest possible Linoz vertical level
696 0 : integer :: ltrop_trop(ncol) !Tropopause level for hybrid case.
697 : logical :: ltrop_linoz_set !Flag that lowest linoz level already found.
698 0 : real(kind_phys) :: trop_output(ncol,pver) !For output purposes only.
699 0 : real(kind_phys) :: trop_linoz_output(ncol,pver) !For output purposes only.
700 0 : real(kind_phys) :: trop_trop_output(ncol,pver) !For output purposes only.
701 :
702 0 : ltrop_linoz(:) = 1 ! Initialize to default value.
703 0 : ltrop_trop(:) = 1 ! Initialize to default value.
704 :
705 0 : LOOP_COL4: do i=1,ncol
706 :
707 : ! Skip column in which the tropopause has already been found.
708 0 : not_found: if (tropLev(i) == NOTFOUND) then
709 :
710 : stobie_min = 1.e10_kind_phys ! An impossibly large number
711 : ltrop_linoz_set = .FALSE.
712 0 : LOOP_LEV: do k=pver,1,-1
713 0 : IF (pmid(i,k) < min_stobie_pressure) cycle
714 0 : shybrid_temp = ALPHA * t(i,k) - Log10(pmid(i,k))
715 : !PJC_NOTE: the units of pmid won't matter, because it is just an additive offset.
716 0 : IF (shybrid_temp<stobie_min) then
717 0 : ltrop_trop(i)=k
718 0 : stobie_min = shybrid_temp
719 : ENDIF
720 0 : IF (pmid(i,k) < max_Linoz_pressure .AND. .NOT. ltrop_linoz_set) THEN
721 0 : ltrop_linoz(i) = k
722 0 : ltrop_linoz_set = .TRUE.
723 : ENDIF
724 : enddo LOOP_LEV
725 :
726 0 : tropLev(i) = MIN(ltrop_trop(i),ltrop_linoz(i))
727 :
728 0 : if (present(tropP)) then
729 0 : tropP(i) = pmid(i,tropLev(i))
730 : endif
731 0 : if (present(tropT)) then
732 0 : tropT(i) = t(i,tropLev(i))
733 : endif
734 0 : if (present(tropZ)) then
735 0 : tropZ(i) = zm(i,tropLev(i))
736 : endif
737 :
738 : endif not_found
739 :
740 : enddo LOOP_COL4
741 :
742 0 : trop_output(:,:)=0._kind_phys
743 0 : trop_linoz_output(:,:)=0._kind_phys
744 0 : trop_trop_output(:,:)=0._kind_phys
745 0 : do i=1,ncol
746 0 : if (tropLev(i)>0) then
747 0 : trop_output(i,tropLev(i))=1._kind_phys
748 0 : trop_linoz_output(i,ltrop_linoz(i))=1._kind_phys
749 0 : trop_trop_output(i,ltrop_trop(i))=1._kind_phys
750 : endif
751 : enddo
752 :
753 0 : if(present(hstobie_trop)) then
754 0 : hstobie_trop(:,:) = trop_output(:,:)
755 : endif
756 :
757 0 : if(present(hstobie_linoz)) then
758 0 : hstobie_linoz(:,:) = trop_linoz_output(:,:)
759 : endif
760 :
761 0 : if(present(hstobie_tropop)) then
762 0 : hstobie_tropop(:,:) = trop_trop_output(:,:)
763 : endif
764 :
765 0 : end subroutine tropopause_hybridstobie
766 :
767 : ! This routine originates with Stobie at NASA Goddard, but does not have a
768 : ! known reference. It was supplied by Philip Cameron-Smith of LLNL.
769 : !
770 0 : subroutine tropopause_stobie(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
771 0 : tropLev, tropP, tropT, tropZ)
772 :
773 : integer, intent(in) :: ncol ! Number of atmospheric columns
774 : integer, intent(in) :: pver ! Number of vertical levels
775 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
776 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
777 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
778 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
779 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
780 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
781 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
782 :
783 : integer, intent(inout) :: tropLev(:) ! tropopause level index
784 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
785 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
786 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
787 :
788 : ! Local Variables
789 : integer :: i
790 : integer :: k
791 : integer :: tLev ! tropopause level
792 : real(kind_phys) :: tP ! tropopause pressure (Pa)
793 0 : real(kind_phys) :: stobie(pver) ! stobie weighted temperature
794 : real(kind_phys) :: sTrop ! stobie value at the tropopause
795 :
796 : ! Iterate over all of the columns.
797 0 : do i = 1, ncol
798 :
799 : ! Skip column in which the tropopause has already been found.
800 0 : if (tropLev(i) == NOTFOUND) then
801 :
802 : ! Caclulate a pressure weighted temperature.
803 0 : stobie(:) = ALPHA * t(i,:) - log10(pmid(i, :))
804 :
805 : ! Search from the bottom up, looking for the first minimum.
806 0 : tLev = -1
807 :
808 0 : do k = pver-1, 1, -1
809 :
810 0 : if (pmid(i, k) <= 4000._kind_phys) then
811 : exit
812 : end if
813 :
814 0 : if (pmid(i, k) >= 55000._kind_phys) then
815 : cycle
816 : end if
817 :
818 0 : if ((tLev == -1) .or. (stobie(k) < sTrop)) then
819 0 : tLev = k
820 0 : tP = pmid(i, k)
821 0 : sTrop = stobie(k)
822 : end if
823 : end do
824 :
825 0 : if (tLev /= -1) then
826 0 : tropLev(i) = tLev
827 :
828 : ! Return the optional outputs
829 0 : if (present(tropP)) tropP(i) = tP
830 :
831 0 : if (present(tropT)) then
832 0 : tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
833 : end if
834 :
835 0 : if (present(tropZ)) then
836 0 : tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
837 : end if
838 : end if
839 : end if
840 : end do
841 :
842 0 : end subroutine tropopause_stobie
843 :
844 :
845 : ! This routine is an implementation of Reichler et al. [2003] done by
846 : ! Reichler and downloaded from his web site. Minimal modifications were
847 : ! made to have the routine work within the CAM framework (i.e. using
848 : ! CAM constants and types).
849 : !
850 : ! NOTE: I am not a big fan of the goto's and multiple returns in this
851 : ! code, but for the moment I have left them to preserve as much of the
852 : ! original and presumably well tested code as possible.
853 : ! UPDATE: The most "obvious" substitutions have been made to replace
854 : ! goto/return statements with cycle/exit. The structure is still
855 : ! somewhat tangled.
856 : ! UPDATE 2: "gamma" renamed to "gam" in order to avoid confusion
857 : ! with the Fortran 2008 intrinsic. "level" argument removed because
858 : ! a physics column is not contiguous, so using explicit dimensions
859 : ! will cause the data to be needlessly copied.
860 : !
861 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
862 : !
863 : ! determination of tropopause height from gridded temperature data
864 : !
865 : ! Reichler, T., M. Dameris, and R. Sausen (2003),
866 : ! Determining the tropopause height from gridded data,
867 : ! Geophys. Res. Lett., 30, 2042, doi:10.1029/2003GL018240, 20.
868 : !
869 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
870 58708800 : subroutine twmo(t, p, plimu, pliml, gam, trp)
871 :
872 : real(kind_phys), intent(in), dimension(:) :: t, p
873 : real(kind_phys), intent(in) :: plimu, pliml, gam
874 : real(kind_phys), intent(out) :: trp
875 :
876 : real(kind_phys), parameter :: deltaz = 2000.0_kind_phys
877 :
878 : real(kind_phys) :: pmk, pm, a, b, tm, dtdp, dtdz
879 : real(kind_phys) :: ag, bg, ptph
880 : real(kind_phys) :: pm0, pmk0, dtdz0
881 : real(kind_phys) :: p2km, asum, aquer
882 : real(kind_phys) :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2
883 : integer :: level
884 : integer :: icount, jj
885 : integer :: j
886 :
887 :
888 58708800 : trp=-99.0_kind_phys ! negative means not valid
889 :
890 : ! initialize start level
891 : ! dt/dz
892 58708800 : level = size(t)
893 58708800 : pmk= .5_kind_phys * (p(level-1)**cnst_kap+p(level)**cnst_kap)
894 58708800 : pm = pmk**(1/cnst_kap)
895 58708800 : a = (t(level-1)-t(level))/(p(level-1)**cnst_kap-p(level)**cnst_kap)
896 58708800 : b = t(level)-(a*p(level)**cnst_kap)
897 58708800 : tm = a * pmk + b
898 58708800 : dtdp = a * cnst_kap * (pm**cnst_ka1)
899 58708800 : dtdz = cnst_faktor*dtdp*pm/tm
900 :
901 860489835 : main_loop: do j=level-1,2,-1
902 : pm0 = pm
903 860338015 : pmk0 = pmk
904 860338015 : dtdz0 = dtdz
905 :
906 : ! dt/dz
907 860338015 : pmk= .5_kind_phys * (p(j-1)**cnst_kap+p(j)**cnst_kap)
908 860338015 : pm = pmk**(1/cnst_kap)
909 860338015 : a = (t(j-1)-t(j))/(p(j-1)**cnst_kap-p(j)**cnst_kap)
910 860338015 : b = t(j)-(a*p(j)**cnst_kap)
911 860338015 : tm = a * pmk + b
912 860338015 : dtdp = a * cnst_kap * (pm**cnst_ka1)
913 860338015 : dtdz = cnst_faktor*dtdp*pm/tm
914 : ! dt/dz valid?
915 860338015 : if (dtdz.le.gam) cycle main_loop ! no, dt/dz < -2 K/km
916 82117522 : if (pm.gt.plimu) cycle main_loop ! no, too low
917 :
918 : ! dtdz is valid, calculate tropopause pressure
919 62963492 : if (dtdz0.lt.gam) then
920 61947215 : ag = (dtdz-dtdz0) / (pmk-pmk0)
921 61947215 : bg = dtdz0 - (ag * pmk0)
922 61947215 : ptph = exp(log((gam-bg)/ag)/cnst_kap)
923 : else
924 : ptph = pm
925 : endif
926 :
927 62963492 : if (ptph.lt.pliml) cycle main_loop
928 62052572 : if (ptph.gt.plimu) cycle main_loop
929 :
930 : ! 2nd test: dtdz above 2 km must not exceed gam
931 62051289 : p2km = ptph + deltaz*(pm/tm)*cnst_faktor ! p at ptph + 2km
932 62051289 : asum = 0.0_kind_phys ! dtdz above
933 62051289 : icount = 0 ! number of levels above
934 :
935 : ! test until apm < p2km
936 190153617 : in_loop: do jj=j,2,-1
937 :
938 190153617 : pmk2 = .5_kind_phys * (p(jj-1)**cnst_kap+p(jj)**cnst_kap) ! p mean ^kappa
939 190153617 : pm2 = pmk2**(1/cnst_kap) ! p mean
940 190153617 : if(pm2.gt.ptph) cycle in_loop ! doesn't happen
941 190153617 : if(pm2.lt.p2km) exit in_loop ! ptropo is valid
942 :
943 131596637 : a2 = (t(jj-1)-t(jj)) ! a
944 131596637 : a2 = a2/(p(jj-1)**cnst_kap-p(jj)**cnst_kap)
945 131596637 : b2 = t(jj)-(a2*p(jj)**cnst_kap) ! b
946 131596637 : tm2 = a2 * pmk2 + b2 ! T mean
947 131596637 : dtdp2 = a2 * cnst_kap * (pm2**(cnst_kap-1)) ! dt/dp
948 131596637 : dtdz2 = cnst_faktor*dtdp2*pm2/tm2
949 131596637 : asum = asum+dtdz2
950 131596637 : icount = icount+1
951 131596637 : aquer = asum/float(icount) ! dt/dz mean
952 :
953 : ! discard ptropo ?
954 190153617 : if (aquer.le.gam) cycle main_loop ! dt/dz above < gam
955 :
956 : enddo in_loop ! test next level
957 :
958 58556980 : trp = ptph
959 860489835 : exit main_loop
960 : enddo main_loop
961 :
962 58708800 : end subroutine twmo
963 :
964 :
965 : ! This routine uses an implementation of Reichler et al. [2003] done by
966 : ! Reichler and downloaded from his web site. This is similar to the WMO
967 : ! routines, but is designed for GCMs with a coarse vertical grid.
968 7479936 : subroutine tropopause_twmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
969 7479936 : tropLev, tropP, tropT, tropZ)
970 :
971 : integer, intent(in) :: ncol ! Number of atmospheric columns
972 : integer, intent(in) :: pver ! Number of vertical levels
973 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
974 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
975 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
976 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
977 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
978 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
979 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
980 : integer, intent(inout) :: tropLev(:) ! tropopause level index
981 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
982 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
983 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
984 :
985 : ! Local Variables
986 : real(kind_phys), parameter :: gam = -0.002_kind_phys ! K/m
987 : real(kind_phys), parameter :: plimu = 45000._kind_phys ! Pa
988 : real(kind_phys), parameter :: pliml = 7500._kind_phys ! Pa
989 :
990 : integer :: i
991 : integer :: k
992 : real(kind_phys) :: tP ! tropopause pressure (Pa)
993 :
994 : ! Iterate over all of the columns.
995 62448768 : do i = 1, ncol
996 :
997 : ! Skip column in which the tropopause has already been found.
998 62448768 : if (tropLev(i) == NOTFOUND) then
999 :
1000 : ! Use the routine from Reichler.
1001 58708800 : call twmo(t(i, :), pmid(i, :), plimu, pliml, gam, tP)
1002 :
1003 : ! if successful, store of the results and find the level and temperature.
1004 58708800 : if (tP > 0) then
1005 :
1006 : ! Find the associated level.
1007 914070955 : do k = pver, 2, -1
1008 914070955 : if (tP >= pint(i, k)) then
1009 58556980 : tropLev(i) = k
1010 58556980 : exit
1011 : end if
1012 : end do
1013 :
1014 : ! Return the optional outputs
1015 58556980 : if (present(tropP)) tropP(i) = tP
1016 :
1017 58556980 : if (present(tropT)) then
1018 46826372 : tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
1019 : end if
1020 :
1021 58556980 : if (present(tropZ)) then
1022 46826372 : tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
1023 : end if
1024 : end if
1025 : end if
1026 : end do
1027 :
1028 3739968 : end subroutine tropopause_twmo
1029 :
1030 : ! This routine implements the WMO definition of the tropopause (WMO, 1957; Seidel and Randel, 2006).
1031 : ! Seidel, D. J., and W. J. Randel (2006),
1032 : ! Variability and trends in the global tropopause estimated from radiosonde data,
1033 : ! J. Geophys. Res., 111, D21101, doi:10.1029/2006JD007363.
1034 : !
1035 : ! This requires that the lapse rate be less than 2 K/km for an altitude range
1036 : ! of 2 km. The search starts at the surface and stops the first time this
1037 : ! criteria is met.
1038 0 : subroutine tropopause_wmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
1039 0 : tropLev, tropP, tropT, tropZ)
1040 :
1041 : integer, intent(in) :: ncol ! Number of atmospheric columns
1042 : integer, intent(in) :: pver ! Number of vertical levels
1043 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
1044 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
1045 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
1046 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
1047 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
1048 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
1049 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
1050 :
1051 : integer, intent(inout) :: tropLev(:) ! tropopause level index
1052 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
1053 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
1054 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
1055 :
1056 : ! Local Variables
1057 : real(kind_phys), parameter :: ztrop_low = 5000._kind_phys ! lowest tropopause level allowed (m)
1058 : real(kind_phys), parameter :: ztrop_high = 20000._kind_phys ! highest tropopause level allowed (m)
1059 : real(kind_phys), parameter :: max_dtdz = 0.002_kind_phys ! max dt/dz for tropopause level (K/m)
1060 : real(kind_phys), parameter :: min_trop_dz = 2000._kind_phys ! min tropopause thickness (m)
1061 :
1062 : integer :: i
1063 : integer :: k
1064 : integer :: k2
1065 : real(kind_phys) :: tP ! tropopause pressure (Pa)
1066 : real(kind_phys) :: dt
1067 :
1068 : ! Iterate over all of the columns.
1069 0 : do i = 1, ncol
1070 :
1071 : ! Skip column in which the tropopause has already been found.
1072 0 : if (tropLev(i) == NOTFOUND) then
1073 :
1074 0 : kloop: do k = pver-1, 2, -1
1075 :
1076 : ! Skip levels below the minimum and stop if nothing is found
1077 : ! before the maximum.
1078 0 : if (zm(i, k) < ztrop_low) then
1079 : cycle kloop
1080 0 : else if (zm(i, k) > ztrop_high) then
1081 : exit kloop
1082 : end if
1083 :
1084 : ! Compare the actual lapse rate to the threshold
1085 0 : dt = t(i, k) - t(i, k-1)
1086 :
1087 0 : if (dt <= (max_dtdz * (zm(i, k-1) - zm(i, k)))) then
1088 :
1089 : ! Make sure that the lapse rate stays below the threshold for the
1090 : ! specified range.
1091 0 : k2loop: do k2 = k-1, 2, -1
1092 0 : if ((zm(i, k2) - zm(i, k)) >= min_trop_dz) then
1093 0 : tP = pmid(i, k)
1094 0 : tropLev(i) = k
1095 0 : exit k2loop
1096 : end if
1097 :
1098 0 : dt = t(i, k) - t(i, k2)
1099 0 : if (dt > (max_dtdz * (zm(i, k2) - zm(i, k)))) then
1100 : exit k2loop
1101 : end if
1102 : end do k2loop
1103 :
1104 0 : if (tropLev(i) == NOTFOUND) then
1105 : cycle kloop
1106 : else
1107 :
1108 : ! Return the optional outputs
1109 0 : if (present(tropP)) tropP(i) = tP
1110 :
1111 0 : if (present(tropT)) then
1112 0 : tropT(i) = tropopause_interpolateT(pver, pmid, t, i, tropLev(i), tP)
1113 : end if
1114 :
1115 0 : if (present(tropZ)) then
1116 0 : tropZ(i) = tropopause_interpolateZ(pint, pmid, zi, zm, phis, i, tropLev(i), tP)
1117 : end if
1118 :
1119 : exit kloop
1120 : end if
1121 : end if
1122 : end do kloop
1123 : end if
1124 : end do
1125 :
1126 0 : end subroutine tropopause_wmo
1127 :
1128 :
1129 : ! This routine searches for the cold point tropopause, and uses a parabolic
1130 : ! fit of the coldest point and two adjacent points to interpolate the cold point
1131 : ! between model levels.
1132 2990736 : subroutine tropopause_cpp(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
1133 2990736 : tropLev, tropP, tropT, tropZ)
1134 :
1135 : integer, intent(in) :: ncol ! Number of atmospheric columns
1136 : integer, intent(in) :: pver ! Number of vertical levels
1137 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
1138 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
1139 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
1140 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
1141 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
1142 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
1143 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
1144 : integer, intent(inout) :: tropLev(:) ! tropopause level index
1145 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
1146 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
1147 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
1148 :
1149 : ! Local Variables
1150 : real(kind_phys), parameter :: ztrop_low = 5000._kind_phys ! lowest tropopause level allowed (m)
1151 : real(kind_phys), parameter :: ztrop_high = 25000._kind_phys ! highest tropopause level allowed (m)
1152 :
1153 : integer :: i
1154 : integer :: k, firstk, lastk
1155 : integer :: k2
1156 : real(kind_phys) :: tZ ! tropopause height (m)
1157 : real(kind_phys) :: tmin
1158 : real(kind_phys) :: f0, f1, f2
1159 : real(kind_phys) :: x0, x1, x2
1160 : real(kind_phys) :: c0, c1, c2
1161 : real(kind_phys) :: a, b, c
1162 :
1163 : ! Iterate over all of the columns.
1164 24969168 : do i = 1, ncol
1165 :
1166 23473800 : firstk = 0
1167 23473800 : lastk = pver+1
1168 :
1169 : ! Skip column in which the tropopause has already been found.
1170 24969168 : if (tropLev(i) == NOTFOUND) then
1171 23473800 : tmin = 1e6_kind_phys
1172 :
1173 524406067 : kloop: do k = pver-1, 2, -1
1174 :
1175 : ! Skip levels below the minimum and stop if nothing is found
1176 : ! before the maximum.
1177 524406067 : if (zm(i, k) < ztrop_low) then
1178 : firstk = k
1179 : cycle kloop
1180 381514883 : else if (zm(i, k) > ztrop_high) then
1181 : lastk = k
1182 : exit kloop
1183 : end if
1184 :
1185 : ! Find the coldest point
1186 358041083 : if (t(i, k) < tmin) then
1187 222384290 : tropLev(i) = k
1188 222384290 : tmin = t(i,k)
1189 : end if
1190 : end do kloop
1191 :
1192 :
1193 : ! If the minimum is at the edge of the search range, then don't
1194 : ! consider this to be a minima
1195 23473800 : if ((tropLev(i) >= (firstk-1)) .or. (tropLev(i) <= (lastk+1))) then
1196 86 : tropLev(i) = NOTFOUND
1197 : else
1198 :
1199 : ! If returning P, Z, or T, then do a parabolic fit using the
1200 : ! cold point and it its 2 surrounding points to interpolate
1201 : ! between model levels.
1202 23473714 : if (present(tropP) .or. present(tropZ) .or. present(tropT)) then
1203 23473714 : f0 = t(i, tropLev(i)-1)
1204 23473714 : f1 = t(i, tropLev(i))
1205 23473714 : f2 = t(i, tropLev(i)+1)
1206 :
1207 23473714 : x0 = zm(i, tropLev(i)-1)
1208 23473714 : x1 = zm(i, tropLev(i))
1209 23473714 : x2 = zm(i, tropLev(i)+1)
1210 :
1211 23473714 : c0 = (x0-x1)*(x0-x2)
1212 23473714 : c1 = (x1-x0)*(x1-x2)
1213 23473714 : c2 = (x2-x0)*(x2-x1)
1214 :
1215 : ! Determine the quadratic coefficients of:
1216 : ! T = a * z^2 - b*z + c
1217 23473714 : a = (f0/c0 + f1/c1 + f2/c2)
1218 23473714 : b = (f0/c0*(x1+x2) + f1/c1*(x0+x2) + f2/c2*(x0+x1))
1219 23473714 : c = f0/c0*x1*x2 + f1/c1*x0*x2 + f2/c2*x0*x1
1220 :
1221 : ! Find the altitude of the minimum temperature
1222 23473714 : tZ = 0.5_kind_phys * b / a
1223 :
1224 : ! The fit should be between the upper and lower points,
1225 : ! so skip the point if the fit fails.
1226 23473714 : if ((tZ >= x0) .or. (tZ <= x2)) then
1227 0 : tropLev(i) = NOTFOUND
1228 : else
1229 :
1230 : ! Return the optional outputs
1231 23473714 : if (present(tropP)) then
1232 23473714 : tropP(i) = tropopause_interpolateP(pver, pmid, zm, i, tropLev(i), tZ)
1233 : end if
1234 :
1235 23473714 : if (present(tropT)) then
1236 23473714 : tropT(i) = a * tZ*tZ - b*tZ + c
1237 : end if
1238 :
1239 23473714 : if (present(tropZ)) then
1240 23473714 : tropZ(i) = tZ
1241 : end if
1242 : end if
1243 : end if
1244 : end if
1245 : end if
1246 : end do
1247 :
1248 1495368 : end subroutine tropopause_cpp
1249 :
1250 : ! Searches all the columns and attempts to identify the "chemical"
1251 : ! tropopause. This is the lapse rate tropopause, backed up by the climatology
1252 : ! if the lapse rate fails to find the tropopause at pressures higher than a certain
1253 : ! threshold. This pressure threshold depends on latitude. Between 50S and 50N,
1254 : ! the climatology is used if the lapse rate tropopause is not found at P > 75 hPa.
1255 : ! At high latitude (poleward of 50), the threshold is increased to 125 hPa to
1256 : ! eliminate false events that are sometimes detected in the cold polar stratosphere.
1257 : !
1258 : ! NOTE: This routine was adapted from code in chemistry.F90 and mo_gasphase_chemdr.F90.
1259 : ! During the CCPP-ization, findChemTrop is now called from tropopause_find_run using method CHEMTROP
1260 : ! and now also returns the standard tropLev, tropP, tropT, tropZ outputs (optional).
1261 : ! The "backup" option is dropped as it is not used anywhere in current CAM.
1262 0 : subroutine tropopause_findChemTrop(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
1263 0 : calday, tropp_p_loc, tropp_days, &
1264 0 : tropLev, tropP, tropT, tropZ, errmsg, errflg)
1265 :
1266 : integer, intent(in) :: ncol ! Number of atmospheric columns
1267 : integer, intent(in) :: pver ! Number of vertical levels
1268 : real(kind_phys), intent(in) :: lat(:) ! Latitudes (radians)
1269 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
1270 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
1271 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
1272 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
1273 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
1274 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
1275 :
1276 : real(kind_phys), intent(in) :: calday ! Day of year including fraction of day
1277 :
1278 : ! Climatological tropopause pressures (Pa), (ncol,ntimes=12).
1279 : real(kind_phys), intent(in) :: tropp_p_loc(:,:)
1280 : real(kind_phys), intent(in) :: tropp_days(:) ! Day-of-year for climo data, 12
1281 :
1282 : integer, intent(out) :: tropLev(:) ! tropopause level index
1283 : real(kind_phys), optional, intent(inout) :: tropP(:) ! tropopause pressure (Pa)
1284 : real(kind_phys), optional, intent(inout) :: tropT(:) ! tropopause temperature (K)
1285 : real(kind_phys), optional, intent(inout) :: tropZ(:) ! tropopause height (m)
1286 :
1287 : character(len=512), intent(out) :: errmsg
1288 : integer, intent(out) :: errflg
1289 :
1290 : ! Local Variable
1291 0 : real(kind_phys) :: dlats(ncol)
1292 : integer :: i
1293 :
1294 0 : errmsg = ' '
1295 0 : errflg = 0
1296 :
1297 : ! First use the lapse rate tropopause.
1298 0 : call tropopause_twmo(ncol, pver, lat, pint, pmid, t, zi, zm, phis, tropLev, tropP, tropT, tropZ)
1299 :
1300 : ! Now check high latitudes (poleward of 50) and set the level to the
1301 : ! climatology if the level was not found or is at P <= 125 hPa.
1302 0 : dlats(:ncol) = lat(:ncol) * cnst_rad2deg ! convert to degrees
1303 :
1304 0 : do i = 1, ncol
1305 0 : if (abs(dlats(i)) > 50._kind_phys) then
1306 0 : if (tropLev(i) .ne. NOTFOUND) then
1307 0 : if (pmid(i, tropLev(i)) <= 12500._kind_phys) then
1308 0 : tropLev(i) = NOTFOUND
1309 : end if
1310 : end if
1311 : end if
1312 : end do
1313 :
1314 : ! Now use the climatology backup
1315 0 : if (any(tropLev(:) == NOTFOUND)) then
1316 : call tropopause_climate(ncol, pver, lat, pint, pmid, t, zi, zm, phis, &
1317 : calday, tropp_p_loc, tropp_days, &
1318 0 : tropLev, tropP, tropT, tropZ)
1319 : end if
1320 :
1321 0 : end subroutine tropopause_findChemTrop
1322 :
1323 : ! This routine interpolates the pressures in the physics state to
1324 : ! find the pressure at the specified tropopause altitude.
1325 23473714 : function tropopause_interpolateP(pver, pmid, zm, icol, tropLev, tropZ)
1326 :
1327 : implicit none
1328 :
1329 : integer, intent(in) :: pver ! Number of vertical levels
1330 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
1331 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
1332 : integer, intent(in) :: icol ! column being processed
1333 : integer, intent(in) :: tropLev ! tropopause level index
1334 : real(kind_phys), optional, intent(in) :: tropZ ! tropopause pressure (m)
1335 : real(kind_phys) :: tropopause_interpolateP
1336 :
1337 : ! Local Variables
1338 : real(kind_phys) :: tropP ! tropopause pressure (Pa)
1339 : real(kind_phys) :: dlogPdZ ! dlog(p)/dZ
1340 :
1341 : ! Interpolate the temperature linearly against log(P)
1342 :
1343 : ! Is the tropopause at the midpoint?
1344 23473714 : if (tropZ == zm(icol, tropLev)) then
1345 0 : tropP = pmid(icol, tropLev)
1346 :
1347 23473714 : else if (tropZ > zm(icol, tropLev)) then
1348 :
1349 : ! It is above the midpoint? Make sure we aren't at the top.
1350 13207819 : if (tropLev > 1) then
1351 26415638 : dlogPdZ = (log(pmid(icol, tropLev)) - log(pmid(icol, tropLev - 1))) / &
1352 26415638 : (zm(icol, tropLev) - zm(icol, tropLev - 1))
1353 13207819 : tropP = pmid(icol, tropLev) + exp((tropZ - zm(icol, tropLev)) * dlogPdZ)
1354 : end if
1355 : else
1356 :
1357 : ! It is below the midpoint. Make sure we aren't at the bottom.
1358 10265895 : if (tropLev < pver) then
1359 10265895 : dlogPdZ = (log(pmid(icol, tropLev + 1)) - log(pmid(icol, tropLev))) / &
1360 20531790 : (zm(icol, tropLev + 1) - zm(icol, tropLev))
1361 10265895 : tropP = pmid(icol, tropLev) + exp((tropZ - zm(icol, tropLev)) * dlogPdZ)
1362 : end if
1363 : end if
1364 :
1365 23473714 : tropopause_interpolateP = tropP
1366 23473714 : end function tropopause_interpolateP
1367 :
1368 :
1369 : ! This routine interpolates the temperatures in the physics state to
1370 : ! find the temperature at the specified tropopause pressure.
1371 46886986 : function tropopause_interpolateT(pver, pmid, t, icol, tropLev, tropP)
1372 :
1373 : implicit none
1374 :
1375 : integer, intent(in) :: pver ! Number of vertical levels
1376 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
1377 : real(kind_phys), intent(in) :: t(:,:) ! Temperature (K)
1378 : integer, intent(in) :: icol ! column being processed
1379 : integer, intent(in) :: tropLev ! tropopause level index
1380 : real(kind_phys), optional, intent(in) :: tropP ! tropopause pressure (Pa)
1381 : real(kind_phys) :: tropopause_interpolateT
1382 :
1383 : ! Local Variables
1384 : real(kind_phys) :: tropT ! tropopause temperature (K)
1385 : real(kind_phys) :: dTdlogP ! dT/dlog(P)
1386 :
1387 : ! Interpolate the temperature linearly against log(P)
1388 :
1389 : ! Is the tropopause at the midpoint?
1390 46886986 : if (tropP == pmid(icol, tropLev)) then
1391 0 : tropT = t(icol, tropLev)
1392 :
1393 46886986 : else if (tropP < pmid(icol, tropLev)) then
1394 :
1395 : ! It is above the midpoint? Make sure we aren't at the top.
1396 23471692 : if (tropLev > 1) then
1397 46943384 : dTdlogP = (t(icol, tropLev) - t(icol, tropLev - 1)) / &
1398 46943384 : (log(pmid(icol, tropLev)) - log(pmid(icol, tropLev - 1)))
1399 23471692 : tropT = t(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dTdlogP
1400 : end if
1401 : else
1402 :
1403 : ! It is below the midpoint. Make sure we aren't at the bottom.
1404 23415294 : if (tropLev < pver) then
1405 23415294 : dTdlogP = (t(icol, tropLev + 1) - t(icol, tropLev)) / &
1406 46830588 : (log(pmid(icol, tropLev + 1)) - log(pmid(icol, tropLev)))
1407 23415294 : tropT = t(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dTdlogP
1408 : end if
1409 : end if
1410 :
1411 46886986 : tropopause_interpolateT = tropT
1412 46886986 : end function tropopause_interpolateT
1413 :
1414 :
1415 : ! This routine interpolates the geopotential height in the physics state to
1416 : ! find the geopotential height at the specified tropopause pressure.
1417 46886986 : function tropopause_interpolateZ(pint, pmid, zi, zm, phis, icol, tropLev, tropP)
1418 :
1419 : real(kind_phys), intent(in) :: pint(:,:) ! Interface pressures (Pa)
1420 : real(kind_phys), intent(in) :: pmid(:,:) ! Midpoint pressures (Pa)
1421 : real(kind_phys), intent(in) :: zi(:,:) ! Geopotential height above surface at interfaces (m)
1422 : real(kind_phys), intent(in) :: zm(:,:) ! Geopotential height above surface at midpoints (m)
1423 : real(kind_phys), intent(in) :: phis(:) ! Surface geopotential (m2 s-2)
1424 : integer, intent(in) :: icol ! column being processed
1425 : integer, intent(in) :: tropLev ! tropopause level index
1426 : real(kind_phys), optional, intent(in) :: tropP ! tropopause pressure (Pa)
1427 : real(kind_phys) :: tropopause_interpolateZ
1428 :
1429 : ! Local Variables
1430 : real(kind_phys) :: tropZ ! tropopause geopotential height (m)
1431 : real(kind_phys) :: dZdlogP ! dZ/dlog(P)
1432 :
1433 : ! Interpolate the geopotential height linearly against log(P)
1434 :
1435 : ! Is the tropopause at the midpoint?
1436 46886986 : if (tropP == pmid(icol, tropLev)) then
1437 0 : tropZ = zm(icol, tropLev)
1438 :
1439 46886986 : else if (tropP < pmid(icol, tropLev)) then
1440 :
1441 : ! It is above the midpoint? Make sure we aren't at the top.
1442 23471692 : dZdlogP = (zm(icol, tropLev) - zi(icol, tropLev)) / &
1443 23471692 : (log(pmid(icol, tropLev)) - log(pint(icol, tropLev)))
1444 23471692 : tropZ = zm(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dZdlogP
1445 : else
1446 :
1447 : ! It is below the midpoint. Make sure we aren't at the bottom.
1448 46830588 : dZdlogP = (zm(icol, tropLev) - zi(icol, tropLev+1)) / &
1449 46830588 : (log(pmid(icol, tropLev)) - log(pint(icol, tropLev+1)))
1450 23415294 : tropZ = zm(icol, tropLev) + (log(tropP) - log(pmid(icol, tropLev))) * dZdlogP
1451 : end if
1452 :
1453 46886986 : tropopause_interpolateZ = tropZ + phis(icol)*cnst_rga
1454 46886986 : end function tropopause_interpolateZ
1455 : end module tropopause_find
|
