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 : ! These routines are based upon code in the WACCM chemistry module
9 : ! including mo_tropoause.F90 and llnl_set_chem_trop.F90. The code
10 : ! for the Reichler et al. [2003] algorithm is from:
11 : !
12 : ! http://www.gfdl.noaa.gov/~tjr/TROPO/tropocode.htm
13 : !
14 : ! Author: Charles Bardeen
15 : ! Created: April, 2009
16 :
17 : module tropopause
18 : !---------------------------------------------------------------
19 : ! ... variables for the tropopause module
20 : !---------------------------------------------------------------
21 :
22 : use shr_kind_mod, only : r8 => shr_kind_r8
23 : use shr_const_mod, only : pi => shr_const_pi
24 : use ppgrid, only : pcols, pver, begchunk, endchunk
25 : use cam_abortutils, only : endrun
26 : use cam_logfile, only : iulog
27 : use cam_history_support, only : fillvalue
28 : use physics_types, only : physics_state
29 : use physconst, only : cappa, rair, gravit
30 : use spmd_utils, only : masterproc
31 :
32 : implicit none
33 :
34 : private
35 :
36 : public :: tropopause_readnl, tropopause_init, tropopause_find, tropopause_output
37 : public :: tropopause_findChemTrop
38 : public :: TROP_ALG_NONE, TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
39 : public :: TROP_ALG_STOBIE, TROP_ALG_HYBSTOB, TROP_ALG_TWMO, TROP_ALG_WMO
40 : public :: TROP_ALG_CPP
41 : public :: NOTFOUND
42 :
43 : save
44 :
45 : ! These parameters define and enumeration to be used to define the primary
46 : ! and backup algorithms to be used with the tropopause_find() method. The
47 : ! backup algorithm is meant to provide a solution when the primary algorithm
48 : ! fail. The algorithms that can't fail are: TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
49 : ! and TROP_ALG_STOBIE.
50 : integer, parameter :: TROP_ALG_NONE = 1 ! Don't evaluate
51 : integer, parameter :: TROP_ALG_ANALYTIC = 2 ! Analytic Expression
52 : integer, parameter :: TROP_ALG_CLIMATE = 3 ! Climatology
53 : integer, parameter :: TROP_ALG_STOBIE = 4 ! Stobie Algorithm
54 : integer, parameter :: TROP_ALG_TWMO = 5 ! WMO Definition, Reichler et al. [2003]
55 : integer, parameter :: TROP_ALG_WMO = 6 ! WMO Definition
56 : integer, parameter :: TROP_ALG_HYBSTOB = 7 ! Hybrid Stobie Algorithm
57 : integer, parameter :: TROP_ALG_CPP = 8 ! Cold Point Parabolic
58 :
59 : integer, parameter :: TROP_NALG = 8 ! Number of Algorithms
60 : character,parameter :: TROP_LETTER(TROP_NALG) = (/ ' ', 'A', 'C', 'S', 'T', 'W', 'H', 'F' /)
61 : ! unique identifier for output, don't use P
62 :
63 : ! These variables should probably be controlled by namelist entries.
64 : logical ,parameter :: output_all = .False. ! output tropopause info from all algorithms
65 : integer ,parameter :: default_primary = TROP_ALG_TWMO ! default primary algorithm
66 : integer ,parameter :: default_backup = TROP_ALG_CLIMATE ! default backup algorithm
67 :
68 : ! Namelist variables
69 : character(len=256) :: tropopause_climo_file = 'trop_climo' ! absolute filepath of climatology file
70 :
71 : ! These variables are used to store the climatology data.
72 : real(r8) :: days(12) ! days in the climatology
73 : real(r8), pointer :: tropp_p_loc(:,:,:) ! climatological tropopause pressures
74 :
75 : integer, parameter :: NOTFOUND = -1
76 :
77 : real(r8),parameter :: ALPHA = 0.03_r8
78 :
79 : ! physical constants
80 : ! These constants are set in module variables rather than as parameters
81 : ! to support the aquaplanet mode in which the constants have values determined
82 : ! by the experiment protocol
83 : real(r8) :: cnst_kap ! = cappa
84 : real(r8) :: cnst_faktor ! = -gravit/rair
85 : real(r8) :: cnst_ka1 ! = cnst_kap - 1._r8
86 :
87 : !================================================================================================
88 : contains
89 : !================================================================================================
90 :
91 : ! Read namelist variables.
92 60360 : subroutine tropopause_readnl(nlfile)
93 :
94 : use namelist_utils, only: find_group_name
95 : use units, only: getunit, freeunit
96 : use mpishorthand
97 :
98 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
99 :
100 : ! Local variables
101 : integer :: unitn, ierr
102 : character(len=*), parameter :: subname = 'tropopause_readnl'
103 :
104 : namelist /tropopause_nl/ tropopause_climo_file
105 : !-----------------------------------------------------------------------------
106 :
107 1536 : if (masterproc) then
108 2 : unitn = getunit()
109 2 : open( unitn, file=trim(nlfile), status='old' )
110 2 : call find_group_name(unitn, 'tropopause_nl', status=ierr)
111 2 : if (ierr == 0) then
112 2 : read(unitn, tropopause_nl, iostat=ierr)
113 2 : if (ierr /= 0) then
114 0 : call endrun(subname // ':: ERROR reading namelist')
115 : end if
116 : end if
117 2 : close(unitn)
118 2 : call freeunit(unitn)
119 : end if
120 :
121 : #ifdef SPMD
122 : ! Broadcast namelist variables
123 1536 : call mpibcast(tropopause_climo_file, len(tropopause_climo_file), mpichar, 0, mpicom)
124 : #endif
125 :
126 1536 : end subroutine tropopause_readnl
127 :
128 :
129 : ! This routine is called during intialization and must be called before the
130 : ! other methods in this module can be used. Its main tasks are to read in the
131 : ! climatology from a file and to define the output fields. Much of this code
132 : ! is taken from mo_tropopause.
133 1536 : subroutine tropopause_init()
134 :
135 : use cam_history, only: addfld, horiz_only
136 :
137 :
138 : implicit none
139 :
140 : ! define physical constants
141 1536 : cnst_kap = cappa
142 1536 : cnst_faktor = -gravit/rair
143 1536 : cnst_ka1 = cnst_kap - 1._r8
144 :
145 : ! Define the output fields.
146 1536 : call addfld('TROP_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure', flag_xyfill=.True.)
147 1536 : call addfld('TROP_T', horiz_only, 'A', 'K', 'Tropopause Temperature', flag_xyfill=.True.)
148 1536 : call addfld('TROP_Z', horiz_only, 'A', 'm', 'Tropopause Height', flag_xyfill=.True.)
149 3072 : call addfld('TROP_DZ', (/ 'lev' /), 'A', 'm', 'Relative Tropopause Height')
150 3072 : call addfld('TROP_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Probabilty')
151 1536 : call addfld('TROP_FD', horiz_only, 'A', 'probability', 'Tropopause Found')
152 :
153 1536 : call addfld('TROPP_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (primary)', flag_xyfill=.True.)
154 1536 : call addfld('TROPP_T', horiz_only, 'A', 'K', 'Tropopause Temperature (primary)', flag_xyfill=.True.)
155 1536 : call addfld('TROPP_Z', horiz_only, 'A', 'm', 'Tropopause Height (primary)', flag_xyfill=.True.)
156 3072 : call addfld('TROPP_DZ', (/ 'lev' /), 'A', 'm', 'Relative Tropopause Height (primary)')
157 3072 : call addfld('TROPP_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (primary)')
158 1536 : call addfld('TROPP_FD', horiz_only, 'A', 'probability', 'Tropopause Found (primary)')
159 :
160 1536 : call addfld('TROPF_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (cold point)', flag_xyfill=.True.)
161 1536 : call addfld('TROPF_T', horiz_only, 'A', 'K', 'Tropopause Temperature (cold point)', flag_xyfill=.True.)
162 1536 : call addfld('TROPF_Z', horiz_only, 'A', 'm', 'Tropopause Height (cold point)', flag_xyfill=.True.)
163 3072 : call addfld('TROPF_DZ', (/ 'lev' /), 'A', 'm', 'Relative Tropopause Height (cold point)', flag_xyfill=.True.)
164 3072 : call addfld('TROPF_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (cold point)')
165 1536 : call addfld('TROPF_FD', horiz_only, 'A', 'probability', 'Tropopause Found (cold point)')
166 :
167 3072 : call addfld( 'hstobie_trop', (/ 'lev' /), 'I', 'fraction of model time', 'Lowest level with stratospheric chemsitry')
168 3072 : call addfld( 'hstobie_linoz', (/ 'lev' /), 'I', 'fraction of model time', 'Lowest possible Linoz level')
169 : call addfld( 'hstobie_tropop', (/ 'lev' /), 'I', 'fraction of model time', &
170 3072 : 'Troposphere boundary calculated in chemistry' )
171 :
172 : ! If requested, be prepared to output results from all of the methods.
173 : if (output_all) then
174 : call addfld('TROPA_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (analytic)', flag_xyfill=.True.)
175 : call addfld('TROPA_T', horiz_only, 'A', 'K', 'Tropopause Temperature (analytic)', flag_xyfill=.True.)
176 : call addfld('TROPA_Z', horiz_only, 'A', 'm', 'Tropopause Height (analytic)', flag_xyfill=.True.)
177 : call addfld('TROPA_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (analytic)')
178 : call addfld('TROPA_FD', horiz_only, 'A', 'probability', 'Tropopause Found (analytic)')
179 :
180 : call addfld('TROPC_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (climatology)', flag_xyfill=.True.)
181 : call addfld('TROPC_T', horiz_only, 'A', 'K', 'Tropopause Temperature (climatology)', flag_xyfill=.True.)
182 : call addfld('TROPC_Z', horiz_only, 'A', 'm', 'Tropopause Height (climatology)', flag_xyfill=.True.)
183 : call addfld('TROPC_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (climatology)')
184 : call addfld('TROPC_FD', horiz_only, 'A', 'probability', 'Tropopause Found (climatology)')
185 :
186 : call addfld('TROPS_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (stobie)', flag_xyfill=.True.)
187 : call addfld('TROPS_T', horiz_only, 'A', 'K', 'Tropopause Temperature (stobie)', flag_xyfill=.True.)
188 : call addfld('TROPS_Z', horiz_only, 'A', 'm', 'Tropopause Height (stobie)', flag_xyfill=.True.)
189 : call addfld('TROPS_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (stobie)')
190 : call addfld('TROPS_FD', horiz_only, 'A', 'probability', 'Tropopause Found (stobie)')
191 :
192 : call addfld('TROPT_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (twmo)', flag_xyfill=.True.)
193 : call addfld('TROPT_T', horiz_only, 'A', 'K', 'Tropopause Temperature (twmo)', flag_xyfill=.True.)
194 : call addfld('TROPT_Z', horiz_only, 'A', 'm', 'Tropopause Height (twmo)', flag_xyfill=.True.)
195 : call addfld('TROPT_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (twmo)')
196 : call addfld('TROPT_FD', horiz_only, 'A', 'probability', 'Tropopause Found (twmo)')
197 :
198 : call addfld('TROPW_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (WMO)', flag_xyfill=.True.)
199 : call addfld('TROPW_T', horiz_only, 'A', 'K', 'Tropopause Temperature (WMO)', flag_xyfill=.True.)
200 : call addfld('TROPW_Z', horiz_only, 'A', 'm', 'Tropopause Height (WMO)', flag_xyfill=.True.)
201 : call addfld('TROPW_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (WMO)')
202 : call addfld('TROPW_FD', horiz_only, 'A', 'probability', 'Tropopause Found (WMO)')
203 :
204 : call addfld('TROPH_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (Hybrid Stobie)', flag_xyfill=.True.)
205 : call addfld('TROPH_T', horiz_only, 'A', 'K', 'Tropopause Temperature (Hybrid Stobie)', flag_xyfill=.True.)
206 : call addfld('TROPH_Z', horiz_only, 'A', 'm', 'Tropopause Height (Hybrid Stobie)', flag_xyfill=.True.)
207 : call addfld('TROPH_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (Hybrid Stobie)')
208 : call addfld('TROPH_FD', horiz_only, 'A', 'probability', 'Tropopause Found (Hybrid Stobie)')
209 : end if
210 :
211 :
212 1536 : call tropopause_read_file()
213 :
214 :
215 1536 : end subroutine tropopause_init
216 :
217 :
218 1536 : subroutine tropopause_read_file
219 : !------------------------------------------------------------------
220 : ! ... initialize upper boundary values
221 : !------------------------------------------------------------------
222 1536 : use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
223 : use dyn_grid, only : get_dyn_grid_parm
224 : use phys_grid, only : get_ncols_p, get_rlat_all_p, get_rlon_all_p
225 : use ioFileMod, only : getfil
226 : use time_manager, only : get_calday
227 : use physconst, only : pi
228 : use cam_pio_utils, only: cam_pio_openfile
229 : use pio, only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, &
230 : pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite
231 :
232 : !------------------------------------------------------------------
233 : ! ... local variables
234 : !------------------------------------------------------------------
235 : integer :: i, j, n
236 : integer :: ierr
237 : type(file_desc_t) :: pio_id
238 : integer :: dimid
239 : type(var_desc_t) :: vid
240 : integer :: nlon, nlat, ntimes
241 : integer :: start(3)
242 : integer :: count(3)
243 : integer, parameter :: dates(12) = (/ 116, 214, 316, 415, 516, 615, &
244 : 716, 816, 915, 1016, 1115, 1216 /)
245 : integer :: plon, plat
246 : type(interp_type) :: lon_wgts, lat_wgts
247 1536 : real(r8), allocatable :: tropp_p_in(:,:,:)
248 1536 : real(r8), allocatable :: lat(:)
249 1536 : real(r8), allocatable :: lon(:)
250 : real(r8) :: to_lats(pcols), to_lons(pcols)
251 : real(r8), parameter :: d2r=pi/180._r8, zero=0._r8, twopi=pi*2._r8
252 : character(len=256) :: locfn
253 : integer :: c, ncols
254 :
255 :
256 : plon = get_dyn_grid_parm('plon')
257 : plat = get_dyn_grid_parm('plat')
258 :
259 :
260 : !-----------------------------------------------------------------------
261 : ! ... open netcdf file
262 : !-----------------------------------------------------------------------
263 1536 : call getfil (tropopause_climo_file, locfn, 0)
264 1536 : call cam_pio_openfile(pio_id, trim(locfn), PIO_NOWRITE)
265 :
266 : !-----------------------------------------------------------------------
267 : ! ... get time dimension
268 : !-----------------------------------------------------------------------
269 1536 : ierr = pio_inq_dimid( pio_id, 'time', dimid )
270 1536 : ierr = pio_inq_dimlen( pio_id, dimid, ntimes )
271 1536 : if( ntimes /= 12 )then
272 0 : write(iulog,*) 'tropopause_init: number of months = ',ntimes,'; expecting 12'
273 0 : call endrun
274 : end if
275 : !-----------------------------------------------------------------------
276 : ! ... get latitudes
277 : !-----------------------------------------------------------------------
278 1536 : ierr = pio_inq_dimid( pio_id, 'lat', dimid )
279 1536 : ierr = pio_inq_dimlen( pio_id, dimid, nlat )
280 4608 : allocate( lat(nlat), stat=ierr )
281 1536 : if( ierr /= 0 ) then
282 0 : write(iulog,*) 'tropopause_init: lat allocation error = ',ierr
283 0 : call endrun
284 : end if
285 1536 : ierr = pio_inq_varid( pio_id, 'lat', vid )
286 1536 : ierr = pio_get_var( pio_id, vid, lat )
287 115200 : lat(:nlat) = lat(:nlat) * d2r
288 : !-----------------------------------------------------------------------
289 : ! ... get longitudes
290 : !-----------------------------------------------------------------------
291 1536 : ierr = pio_inq_dimid( pio_id, 'lon', dimid )
292 1536 : ierr = pio_inq_dimlen( pio_id, dimid, nlon )
293 4608 : allocate( lon(nlon), stat=ierr )
294 1536 : if( ierr /= 0 ) then
295 0 : write(iulog,*) 'tropopause_init: lon allocation error = ',ierr
296 0 : call endrun
297 : end if
298 1536 : ierr = pio_inq_varid( pio_id, 'lon', vid )
299 1536 : ierr = pio_get_var( pio_id, vid, lon )
300 224256 : lon(:nlon) = lon(:nlon) * d2r
301 :
302 : !------------------------------------------------------------------
303 : ! ... allocate arrays
304 : !------------------------------------------------------------------
305 7680 : allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr )
306 1536 : if( ierr /= 0 ) then
307 0 : write(iulog,*) 'tropopause_init: tropp_p_in allocation error = ',ierr
308 0 : call endrun
309 : end if
310 : !------------------------------------------------------------------
311 : ! ... read in the tropopause pressure
312 : !------------------------------------------------------------------
313 1536 : ierr = pio_inq_varid( pio_id, 'trop_p', vid )
314 1536 : start = (/ 1, 1, 1 /)
315 6144 : count = (/ nlon, nlat, ntimes /)
316 1536 : ierr = pio_get_var( pio_id, vid, start, count, tropp_p_in )
317 :
318 : !------------------------------------------------------------------
319 : ! ... close the netcdf file
320 : !------------------------------------------------------------------
321 1536 : call pio_closefile( pio_id )
322 :
323 : !--------------------------------------------------------------------
324 : ! ... regrid
325 : !--------------------------------------------------------------------
326 :
327 6144 : allocate( tropp_p_loc(pcols,begchunk:endchunk,ntimes), stat=ierr )
328 :
329 1536 : if( ierr /= 0 ) then
330 0 : write(iulog,*) 'tropopause_init: tropp_p_loc allocation error = ',ierr
331 0 : call endrun
332 : end if
333 :
334 7728 : do c=begchunk,endchunk
335 6192 : ncols = get_ncols_p(c)
336 6192 : call get_rlat_all_p(c, pcols, to_lats)
337 6192 : call get_rlon_all_p(c, pcols, to_lons)
338 6192 : call lininterp_init(lon, nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
339 6192 : call lininterp_init(lat, nlat, to_lats, ncols, 1, lat_wgts)
340 80496 : do n=1,ntimes
341 80496 : call lininterp(tropp_p_in(:,:,n), nlon, nlat, tropp_p_loc(1:ncols,c,n), ncols, lon_wgts, lat_wgts)
342 : end do
343 6192 : call lininterp_finish(lon_wgts)
344 7728 : call lininterp_finish(lat_wgts)
345 : end do
346 1536 : deallocate(lon)
347 1536 : deallocate(lat)
348 1536 : deallocate(tropp_p_in)
349 :
350 : !--------------------------------------------------------
351 : ! ... initialize the monthly day of year times
352 : !--------------------------------------------------------
353 :
354 19968 : do n = 1,12
355 19968 : days(n) = get_calday( dates(n), 0 )
356 : end do
357 1536 : if (masterproc) then
358 2 : write(iulog,*) 'tropopause_init : days'
359 2 : write(iulog,'(1p,5g15.8)') days(:)
360 : endif
361 :
362 3072 : end subroutine tropopause_read_file
363 :
364 :
365 : ! This analytic expression closely matches the mean tropopause determined
366 : ! by the NCEP reanalysis and has been used by the radiation code.
367 0 : subroutine tropopause_analytic(pstate, tropLev, tropP, tropT, tropZ)
368 :
369 : implicit none
370 :
371 : type(physics_state), intent(in) :: pstate
372 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
373 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
374 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
375 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
376 :
377 : ! Local Variables
378 : integer :: i
379 : integer :: k
380 : integer :: ncol ! number of columns in the chunk
381 : integer :: lchnk ! chunk identifier
382 : real(r8) :: tP ! tropopause pressure (Pa)
383 :
384 : ! Information about the chunk.
385 0 : lchnk = pstate%lchnk
386 0 : ncol = pstate%ncol
387 :
388 : ! Iterate over all of the columns.
389 0 : do i = 1, ncol
390 :
391 : ! Skip column in which the tropopause has already been found.
392 0 : if (tropLev(i) == NOTFOUND) then
393 :
394 : ! Calculate the pressure of the tropopause.
395 0 : tP = (25000.0_r8 - 15000.0_r8 * (cos(pstate%lat(i)))**2)
396 :
397 : ! Find the level that contains the tropopause.
398 0 : do k = pver, 2, -1
399 0 : if (tP >= pstate%pint(i, k)) then
400 0 : tropLev(i) = k
401 0 : exit
402 : end if
403 : end do
404 :
405 : ! Return the optional outputs
406 0 : if (present(tropP)) tropP(i) = tP
407 :
408 0 : if (present(tropT)) then
409 0 : tropT(i) = tropopause_interpolateT(pstate, i, tropLev(i), tP)
410 : end if
411 :
412 0 : if (present(tropZ)) then
413 0 : tropZ(i) = tropopause_interpolateZ(pstate, i, tropLev(i), tP)
414 : end if
415 : end if
416 : end do
417 1536 : end subroutine tropopause_analytic
418 :
419 :
420 : ! Read the tropopause pressure in from a file containging a climatology. The
421 : ! data is interpolated to the current dat of year and latitude.
422 : !
423 : ! NOTE: The data is read in during tropopause_init and stored in the module
424 : ! variable trop
425 171983 : subroutine tropopause_climate(pstate, tropLev, tropP, tropT, tropZ)
426 : use time_manager, only : get_curr_calday
427 :
428 : implicit none
429 :
430 : type(physics_state), intent(in) :: pstate
431 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
432 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
433 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
434 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
435 :
436 : ! Local Variables
437 : integer :: i
438 : integer :: k
439 : integer :: m
440 : integer :: ncol ! number of columns in the chunk
441 : integer :: lchnk ! chunk identifier
442 : real(r8) :: tP ! tropopause pressure (Pa)
443 : real(r8) :: calday ! day of year including fraction
444 : real(r8) :: dels
445 : integer :: last
446 : integer :: next
447 :
448 : ! Information about the chunk.
449 171983 : lchnk = pstate%lchnk
450 171983 : ncol = pstate%ncol
451 :
452 : ! If any columns remain to be indentified, the nget the current
453 : ! day from the calendar.
454 :
455 2598827 : if (any(tropLev == NOTFOUND)) then
456 :
457 : ! Determine the calendar day.
458 171983 : calday = get_curr_calday()
459 :
460 : !--------------------------------------------------------
461 : ! ... setup the time interpolation
462 : !--------------------------------------------------------
463 171983 : if( calday < days(1) ) then
464 171983 : next = 1
465 171983 : last = 12
466 171983 : dels = (365._r8 + calday - days(12)) / (365._r8 + days(1) - days(12))
467 0 : else if( calday >= days(12) ) then
468 0 : next = 1
469 0 : last = 12
470 0 : dels = (calday - days(12)) / (365._r8 + days(1) - days(12))
471 : else
472 0 : do m = 11,1,-1
473 0 : if( calday >= days(m) ) then
474 : exit
475 : end if
476 : end do
477 0 : last = m
478 0 : next = m + 1
479 0 : dels = (calday - days(m)) / (days(m+1) - days(m))
480 : end if
481 :
482 171983 : dels = max( min( 1._r8,dels ),0._r8 )
483 :
484 :
485 : ! Iterate over all of the columns.
486 2725279 : do i = 1, ncol
487 :
488 : ! Skip column in which the tropopause has already been found.
489 2725279 : if (tropLev(i) == NOTFOUND) then
490 :
491 : !--------------------------------------------------------
492 : ! ... get tropopause level from climatology
493 : !--------------------------------------------------------
494 : ! Interpolate the tropopause pressure.
495 0 : tP = tropp_p_loc(i,lchnk,last) &
496 56617 : + dels * (tropp_p_loc(i,lchnk,next) - tropp_p_loc(i,lchnk,last))
497 :
498 : ! Find the associated level.
499 2668505 : do k = pver, 2, -1
500 2668505 : if (tP >= pstate%pint(i, k)) then
501 56617 : tropLev(i) = k
502 56617 : exit
503 : end if
504 : end do
505 :
506 : ! Return the optional outputs
507 56617 : if (present(tropP)) tropP(i) = tP
508 :
509 56617 : if (present(tropT)) then
510 5049 : tropT(i) = tropopause_interpolateT(pstate, i, tropLev(i), tP)
511 : end if
512 :
513 56617 : if (present(tropZ)) then
514 5049 : tropZ(i) = tropopause_interpolateZ(pstate, i, tropLev(i), tP)
515 : end if
516 : end if
517 : end do
518 : end if
519 :
520 171983 : return
521 171983 : end subroutine tropopause_climate
522 :
523 : !-----------------------------------------------------------------------
524 : !-----------------------------------------------------------------------
525 0 : subroutine tropopause_hybridstobie(pstate, tropLev, tropP, tropT, tropZ)
526 171983 : use cam_history, only : outfld
527 :
528 : !-----------------------------------------------------------------------
529 : ! Originally written by Philip Cameron-Smith, LLNL
530 : !
531 : ! Stobie-Linoz hybrid: the highest altitude of
532 : ! a) Stobie algorithm, or
533 : ! b) minimum Linoz pressure.
534 : !
535 : ! NOTE: the ltrop(i) gridbox itself is assumed to be a STRATOSPHERIC gridbox.
536 : !-----------------------------------------------------------------------
537 : !-----------------------------------------------------------------------
538 : ! ... Local variables
539 : !-----------------------------------------------------------------------
540 :
541 : implicit none
542 :
543 : type(physics_state), intent(in) :: pstate
544 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
545 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
546 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
547 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
548 :
549 : real(r8),parameter :: min_Stobie_Pressure= 40.E2_r8 !For case 2 & 4. [Pa]
550 : real(r8),parameter :: max_Linoz_Pressure =208.E2_r8 !For case 4. [Pa]
551 :
552 : integer :: i, k, ncol
553 : real(r8) :: stobie_min, shybrid_temp !temporary variable for case 2 & 3.
554 : integer :: ltrop_linoz(pcols) !Lowest possible Linoz vertical level
555 : integer :: ltrop_trop(pcols) !Tropopause level for hybrid case.
556 : logical :: ltrop_linoz_set !Flag that lowest linoz level already found.
557 : real(r8) :: trop_output(pcols,pver) !For output purposes only.
558 : real(r8) :: trop_linoz_output(pcols,pver) !For output purposes only.
559 : real(r8) :: trop_trop_output(pcols,pver) !For output purposes only.
560 :
561 : ! write(iulog,*) 'In set_chem_trop, o3_ndx =',o3_ndx
562 0 : ltrop_linoz(:) = 1 ! Initialize to default value.
563 0 : ltrop_trop(:) = 1 ! Initialize to default value.
564 0 : ncol = pstate%ncol
565 :
566 0 : LOOP_COL4: do i=1,ncol
567 :
568 : ! Skip column in which the tropopause has already been found.
569 0 : not_found: if (tropLev(i) == NOTFOUND) then
570 :
571 : stobie_min = 1.e10_r8 ! An impossibly large number
572 : ltrop_linoz_set = .FALSE.
573 0 : LOOP_LEV: do k=pver,1,-1
574 0 : IF (pstate%pmid(i,k) < min_stobie_pressure) cycle
575 0 : shybrid_temp = ALPHA * pstate%t(i,k) - Log10(pstate%pmid(i,k))
576 : !PJC_NOTE: the units of pmid won't matter, because it is just an additive offset.
577 0 : IF (shybrid_temp<stobie_min) then
578 0 : ltrop_trop(i)=k
579 0 : stobie_min = shybrid_temp
580 : ENDIF
581 0 : IF (pstate%pmid(i,k) < max_Linoz_pressure .AND. .NOT. ltrop_linoz_set) THEN
582 0 : ltrop_linoz(i) = k
583 0 : ltrop_linoz_set = .TRUE.
584 : ENDIF
585 : enddo LOOP_LEV
586 :
587 0 : tropLev(i) = MIN(ltrop_trop(i),ltrop_linoz(i))
588 :
589 0 : if (present(tropP)) then
590 0 : tropP(i) = pstate%pmid(i,tropLev(i))
591 : endif
592 0 : if (present(tropT)) then
593 0 : tropT(i) = pstate% t(i,tropLev(i))
594 : endif
595 0 : if (present(tropZ)) then
596 0 : tropZ(i) = pstate% zm(i,tropLev(i))
597 : endif
598 :
599 : endif not_found
600 :
601 : enddo LOOP_COL4
602 :
603 0 : trop_output(:,:)=0._r8
604 0 : trop_linoz_output(:,:)=0._r8
605 0 : trop_trop_output(:,:)=0._r8
606 0 : do i=1,ncol
607 0 : if (tropLev(i)>0) then
608 0 : trop_output(i,tropLev(i))=1._r8
609 0 : trop_linoz_output(i,ltrop_linoz(i))=1._r8
610 0 : trop_trop_output(i,ltrop_trop(i))=1._r8
611 : endif
612 : enddo
613 :
614 0 : call outfld( 'hstobie_trop', trop_output(:ncol,:), ncol, pstate%lchnk )
615 0 : call outfld( 'hstobie_linoz', trop_linoz_output(:ncol,:), ncol, pstate%lchnk )
616 0 : call outfld( 'hstobie_tropop', trop_trop_output(:ncol,:), ncol, pstate%lchnk )
617 :
618 0 : endsubroutine tropopause_hybridstobie
619 :
620 : ! This routine originates with Stobie at NASA Goddard, but does not have a
621 : ! known reference. It was supplied by Philip Cameron-Smith of LLNL.
622 : !
623 0 : subroutine tropopause_stobie(pstate, tropLev, tropP, tropT, tropZ)
624 :
625 : implicit none
626 :
627 : type(physics_state), intent(in) :: pstate
628 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
629 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
630 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
631 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
632 :
633 : ! Local Variables
634 : integer :: i
635 : integer :: k
636 : integer :: ncol ! number of columns in the chunk
637 : integer :: lchnk ! chunk identifier
638 : integer :: tLev ! tropopause level
639 : real(r8) :: tP ! tropopause pressure (Pa)
640 : real(r8) :: stobie(pver) ! stobie weighted temperature
641 : real(r8) :: sTrop ! stobie value at the tropopause
642 :
643 : ! Information about the chunk.
644 0 : lchnk = pstate%lchnk
645 0 : ncol = pstate%ncol
646 :
647 : ! Iterate over all of the columns.
648 0 : do i = 1, ncol
649 :
650 : ! Skip column in which the tropopause has already been found.
651 0 : if (tropLev(i) == NOTFOUND) then
652 :
653 : ! Caclulate a pressure weighted temperature.
654 0 : stobie(:) = ALPHA * pstate%t(i,:) - log10(pstate%pmid(i, :))
655 :
656 : ! Search from the bottom up, looking for the first minimum.
657 : tLev = -1
658 :
659 0 : do k = pver-1, 1, -1
660 :
661 0 : if (pstate%pmid(i, k) <= 4000._r8) then
662 : exit
663 : end if
664 :
665 0 : if (pstate%pmid(i, k) >= 55000._r8) then
666 : cycle
667 : end if
668 :
669 0 : if ((tLev == -1) .or. (stobie(k) < sTrop)) then
670 0 : tLev = k
671 0 : tP = pstate%pmid(i, k)
672 : sTrop = stobie(k)
673 : end if
674 : end do
675 :
676 0 : if (tLev /= -1) then
677 0 : tropLev(i) = tLev
678 :
679 : ! Return the optional outputs
680 0 : if (present(tropP)) tropP(i) = tP
681 :
682 0 : if (present(tropT)) then
683 0 : tropT(i) = tropopause_interpolateT(pstate, i, tropLev(i), tP)
684 : end if
685 :
686 0 : if (present(tropZ)) then
687 0 : tropZ(i) = tropopause_interpolateZ(pstate, i, tropLev(i), tP)
688 : end if
689 : end if
690 : end if
691 : end do
692 :
693 0 : return
694 0 : end subroutine tropopause_stobie
695 :
696 :
697 : ! This routine is an implementation of Reichler et al. [2003] done by
698 : ! Reichler and downloaded from his web site. Minimal modifications were
699 : ! made to have the routine work within the CAM framework (i.e. using
700 : ! CAM constants and types).
701 : !
702 : ! NOTE: I am not a big fan of the goto's and multiple returns in this
703 : ! code, but for the moment I have left them to preserve as much of the
704 : ! original and presumably well tested code as possible.
705 : ! UPDATE: The most "obvious" substitutions have been made to replace
706 : ! goto/return statements with cycle/exit. The structure is still
707 : ! somewhat tangled.
708 : ! UPDATE 2: "gamma" renamed to "gam" in order to avoid confusion
709 : ! with the Fortran 2008 intrinsic. "level" argument removed because
710 : ! a physics column is not contiguous, so using explicit dimensions
711 : ! will cause the data to be needlessly copied.
712 : !
713 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
714 : !
715 : ! determination of tropopause height from gridded temperature data
716 : !
717 : ! reference: Reichler, T., M. Dameris, and R. Sausen (2003)
718 : !
719 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
720 11226600 : subroutine twmo(t, p, plimu, pliml, gam, trp)
721 :
722 : real(r8), intent(in), dimension(:) :: t, p
723 : real(r8), intent(in) :: plimu, pliml, gam
724 : real(r8), intent(out) :: trp
725 :
726 : real(r8), parameter :: deltaz = 2000.0_r8
727 :
728 : real(r8) :: pmk, pm, a, b, tm, dtdp, dtdz
729 : real(r8) :: ag, bg, ptph
730 : real(r8) :: pm0, pmk0, dtdz0
731 : real(r8) :: p2km, asum, aquer
732 : real(r8) :: pmk2, pm2, a2, b2, tm2, dtdp2, dtdz2
733 : integer :: level
734 : integer :: icount, jj
735 : integer :: j
736 :
737 :
738 11226600 : trp=-99.0_r8 ! negative means not valid
739 :
740 : ! initialize start level
741 : ! dt/dz
742 11226600 : level = size(t)
743 11226600 : pmk= .5_r8 * (p(level-1)**cnst_kap+p(level)**cnst_kap)
744 11226600 : pm = pmk**(1/cnst_kap)
745 11226600 : a = (t(level-1)-t(level))/(p(level-1)**cnst_kap-p(level)**cnst_kap)
746 11226600 : b = t(level)-(a*p(level)**cnst_kap)
747 11226600 : tm = a * pmk + b
748 11226600 : dtdp = a * cnst_kap * (pm**cnst_ka1)
749 11226600 : dtdz = cnst_faktor*dtdp*pm/tm
750 :
751 462415622 : main_loop: do j=level-1,2,-1
752 : pm0 = pm
753 462353956 : pmk0 = pmk
754 462353956 : dtdz0 = dtdz
755 :
756 : ! dt/dz
757 462353956 : pmk= .5_r8 * (p(j-1)**cnst_kap+p(j)**cnst_kap)
758 462353956 : pm = pmk**(1/cnst_kap)
759 462353956 : a = (t(j-1)-t(j))/(p(j-1)**cnst_kap-p(j)**cnst_kap)
760 462353956 : b = t(j)-(a*p(j)**cnst_kap)
761 462353956 : tm = a * pmk + b
762 462353956 : dtdp = a * cnst_kap * (pm**cnst_ka1)
763 462353956 : dtdz = cnst_faktor*dtdp*pm/tm
764 : ! dt/dz valid?
765 462353956 : if (dtdz.le.gam) cycle main_loop ! no, dt/dz < -2 K/km
766 51920171 : if (pm.gt.plimu) cycle main_loop ! no, too low
767 :
768 : ! dtdz is valid, calculate tropopause pressure
769 17023674 : if (dtdz0.lt.gam) then
770 13733505 : ag = (dtdz-dtdz0) / (pmk-pmk0)
771 13733505 : bg = dtdz0 - (ag * pmk0)
772 13733505 : ptph = exp(log((gam-bg)/ag)/cnst_kap)
773 : else
774 : ptph = pm
775 : endif
776 :
777 17023674 : if (ptph.lt.pliml) cycle main_loop
778 14661857 : if (ptph.gt.plimu) cycle main_loop
779 :
780 : ! 2nd test: dtdz above 2 km must not exceed gam
781 14656366 : p2km = ptph + deltaz*(pm/tm)*cnst_faktor ! p at ptph + 2km
782 14656366 : asum = 0.0_r8 ! dtdz above
783 14656366 : icount = 0 ! number of levels above
784 :
785 : ! test until apm < p2km
786 82324859 : in_loop: do jj=j,2,-1
787 :
788 82324859 : pmk2 = .5_r8 * (p(jj-1)**cnst_kap+p(jj)**cnst_kap) ! p mean ^kappa
789 82324859 : pm2 = pmk2**(1/cnst_kap) ! p mean
790 82324859 : if(pm2.gt.ptph) cycle in_loop ! doesn't happen
791 82324859 : if(pm2.lt.p2km) exit in_loop ! ptropo is valid
792 :
793 71159925 : a2 = (t(jj-1)-t(jj)) ! a
794 71159925 : a2 = a2/(p(jj-1)**cnst_kap-p(jj)**cnst_kap)
795 71159925 : b2 = t(jj)-(a2*p(jj)**cnst_kap) ! b
796 71159925 : tm2 = a2 * pmk2 + b2 ! T mean
797 71159925 : dtdp2 = a2 * cnst_kap * (pm2**(cnst_kap-1)) ! dt/dp
798 71159925 : dtdz2 = cnst_faktor*dtdp2*pm2/tm2
799 71159925 : asum = asum+dtdz2
800 71159925 : icount = icount+1
801 71159925 : aquer = asum/float(icount) ! dt/dz mean
802 :
803 : ! discard ptropo ?
804 82324859 : if (aquer.le.gam) cycle main_loop ! dt/dz above < gam
805 :
806 : enddo in_loop ! test next level
807 :
808 11164934 : trp = ptph
809 462415622 : exit main_loop
810 : enddo main_loop
811 :
812 11226600 : end subroutine twmo
813 :
814 :
815 : ! This routine uses an implementation of Reichler et al. [2003] done by
816 : ! Reichler and downloaded from his web site. This is similar to the WMO
817 : ! routines, but is designed for GCMs with a coarse vertical grid.
818 715176 : subroutine tropopause_twmo(pstate, tropLev, tropP, tropT, tropZ)
819 :
820 : implicit none
821 :
822 : type(physics_state), intent(in) :: pstate
823 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
824 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
825 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
826 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
827 :
828 : ! Local Variables
829 : real(r8), parameter :: gam = -0.002_r8 ! K/m
830 : real(r8), parameter :: plimu = 45000._r8 ! Pa
831 : real(r8), parameter :: pliml = 7500._r8 ! Pa
832 :
833 : integer :: i
834 : integer :: k
835 : integer :: ncol ! number of columns in the chunk
836 : integer :: lchnk ! chunk identifier
837 : real(r8) :: tP ! tropopause pressure (Pa)
838 :
839 : ! Information about the chunk.
840 715176 : lchnk = pstate%lchnk
841 715176 : ncol = pstate%ncol
842 :
843 : ! Iterate over all of the columns.
844 11941776 : do i = 1, ncol
845 :
846 : ! Skip column in which the tropopause has already been found.
847 11941776 : if (tropLev(i) == NOTFOUND) then
848 :
849 : ! Use the routine from Reichler.
850 11226600 : call twmo(pstate%t(i, :), pstate%pmid(i, :), plimu, pliml, gam, tP)
851 :
852 : ! if successful, store of the results and find the level and temperature.
853 11226600 : if (tP > 0) then
854 :
855 : ! Find the associated level.
856 467788385 : do k = pver, 2, -1
857 467788385 : if (tP >= pstate%pint(i, k)) then
858 11164934 : tropLev(i) = k
859 11164934 : exit
860 : end if
861 : end do
862 :
863 : ! Return the optional outputs
864 11164934 : if (present(tropP)) tropP(i) = tP
865 :
866 11164934 : if (present(tropT)) then
867 1836702 : tropT(i) = tropopause_interpolateT(pstate, i, tropLev(i), tP)
868 : end if
869 :
870 11164934 : if (present(tropZ)) then
871 1836702 : tropZ(i) = tropopause_interpolateZ(pstate, i, tropLev(i), tP)
872 : end if
873 : end if
874 : end if
875 : end do
876 :
877 715176 : return
878 : end subroutine tropopause_twmo
879 :
880 : ! This routine implements the WMO definition of the tropopause (WMO, 1957; Seidel and Randel, 2006).
881 : ! This requires that the lapse rate be less than 2 K/km for an altitude range
882 : ! of 2 km. The search starts at the surface and stops the first time this
883 : ! criteria is met.
884 : !
885 : ! NOTE: This code was modeled after the code in mo_tropopause; however, the
886 : ! requirement that dt be greater than 0 was removed and the check to make
887 : ! sure that the lapse rate is maintained for 2 km was added.
888 0 : subroutine tropopause_wmo(pstate, tropLev, tropP, tropT, tropZ)
889 :
890 : implicit none
891 :
892 : type(physics_state), intent(in) :: pstate
893 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
894 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
895 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
896 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
897 :
898 : ! Local Variables
899 : real(r8), parameter :: ztrop_low = 5000._r8 ! lowest tropopause level allowed (m)
900 : real(r8), parameter :: ztrop_high = 20000._r8 ! highest tropopause level allowed (m)
901 : real(r8), parameter :: max_dtdz = 0.002_r8 ! max dt/dz for tropopause level (K/m)
902 : real(r8), parameter :: min_trop_dz = 2000._r8 ! min tropopause thickness (m)
903 :
904 : integer :: i
905 : integer :: k
906 : integer :: k2
907 : integer :: ncol ! number of columns in the chunk
908 : integer :: lchnk ! chunk identifier
909 : real(r8) :: tP ! tropopause pressure (Pa)
910 : real(r8) :: dt
911 :
912 : ! Information about the chunk.
913 0 : lchnk = pstate%lchnk
914 0 : ncol = pstate%ncol
915 :
916 : ! Iterate over all of the columns.
917 0 : do i = 1, ncol
918 :
919 : ! Skip column in which the tropopause has already been found.
920 0 : if (tropLev(i) == NOTFOUND) then
921 :
922 0 : kloop: do k = pver-1, 2, -1
923 :
924 : ! Skip levels below the minimum and stop if nothing is found
925 : ! before the maximum.
926 0 : if (pstate%zm(i, k) < ztrop_low) then
927 : cycle kloop
928 0 : else if (pstate%zm(i, k) > ztrop_high) then
929 : exit kloop
930 : end if
931 :
932 : ! Compare the actual lapse rate to the threshold
933 0 : dt = pstate%t(i, k) - pstate%t(i, k-1)
934 :
935 0 : if (dt <= (max_dtdz * (pstate%zm(i, k-1) - pstate%zm(i, k)))) then
936 :
937 : ! Make sure that the lapse rate stays below the threshold for the
938 : ! specified range.
939 0 : k2loop: do k2 = k-1, 2, -1
940 0 : if ((pstate%zm(i, k2) - pstate%zm(i, k)) >= min_trop_dz) then
941 0 : tP = pstate%pmid(i, k)
942 0 : tropLev(i) = k
943 0 : exit k2loop
944 : end if
945 :
946 0 : dt = pstate%t(i, k) - pstate%t(i, k2)
947 0 : if (dt > (max_dtdz * (pstate%zm(i, k2) - pstate%zm(i, k)))) then
948 : exit k2loop
949 : end if
950 : end do k2loop
951 :
952 0 : if (tropLev(i) == NOTFOUND) then
953 : cycle kloop
954 : else
955 :
956 : ! Return the optional outputs
957 0 : if (present(tropP)) tropP(i) = tP
958 :
959 0 : if (present(tropT)) then
960 0 : tropT(i) = tropopause_interpolateT(pstate, i, tropLev(i), tP)
961 : end if
962 :
963 0 : if (present(tropZ)) then
964 0 : tropZ(i) = tropopause_interpolateZ(pstate, i, tropLev(i), tP)
965 : end if
966 :
967 : exit kloop
968 : end if
969 : end if
970 : end do kloop
971 : end if
972 : end do
973 :
974 0 : return
975 : end subroutine tropopause_wmo
976 :
977 :
978 : ! This routine searches for the cold point tropopause, and uses a parabolic
979 : ! fit of the coldest point and two adjacent points to interpolate the cold point
980 : ! between model levels.
981 58824 : subroutine tropopause_cpp(pstate, tropLev, tropP, tropT, tropZ)
982 :
983 : implicit none
984 :
985 : type(physics_state), intent(in) :: pstate
986 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
987 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
988 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
989 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
990 :
991 : ! Local Variables
992 : real(r8), parameter :: ztrop_low = 5000._r8 ! lowest tropopause level allowed (m)
993 : real(r8), parameter :: ztrop_high = 25000._r8 ! highest tropopause level allowed (m)
994 :
995 : integer :: i
996 : integer :: k, firstk, lastk
997 : integer :: k2
998 : integer :: ncol ! number of columns in the chunk
999 : integer :: lchnk ! chunk identifier
1000 : real(r8) :: tZ ! tropopause height (m)
1001 : real(r8) :: tmin
1002 : real(r8) :: f0, f1, f2
1003 : real(r8) :: x0, x1, x2
1004 : real(r8) :: c0, c1, c2
1005 : real(r8) :: a, b, c
1006 :
1007 : ! Information about the chunk.
1008 58824 : lchnk = pstate%lchnk
1009 58824 : ncol = pstate%ncol
1010 :
1011 : ! Iterate over all of the columns.
1012 982224 : do i = 1, ncol
1013 :
1014 923400 : firstk = 0
1015 923400 : lastk = pver+1
1016 :
1017 : ! Skip column in which the tropopause has already been found.
1018 982224 : if (tropLev(i) == NOTFOUND) then
1019 : tmin = 1e6_r8
1020 :
1021 62137872 : kloop: do k = pver-1, 2, -1
1022 :
1023 : ! Skip levels below the minimum and stop if nothing is found
1024 : ! before the maximum.
1025 62137872 : if (pstate%zm(i, k) < ztrop_low) then
1026 : firstk = k
1027 : cycle kloop
1028 41419936 : else if (pstate%zm(i, k) > ztrop_high) then
1029 : lastk = k
1030 : exit kloop
1031 : end if
1032 :
1033 : ! Find the coldest point
1034 40496536 : if (pstate%t(i, k) < tmin) then
1035 20927322 : tropLev(i) = k
1036 20927322 : tmin = pstate%t(i,k)
1037 : end if
1038 : end do kloop
1039 :
1040 :
1041 : ! If the minimum is at the edge of the search range, then don't
1042 : ! consider this to be a minima
1043 923400 : if ((tropLev(i) >= (firstk-1)) .or. (tropLev(i) <= (lastk+1))) then
1044 7766 : tropLev(i) = NOTFOUND
1045 : else
1046 :
1047 : ! If returning P, Z, or T, then do a parabolic fit using the
1048 : ! cold point and it its 2 surrounding points to interpolate
1049 : ! between model levels.
1050 915634 : if (present(tropP) .or. present(tropZ) .or. present(tropT)) then
1051 915634 : f0 = pstate%t(i, tropLev(i)-1)
1052 915634 : f1 = pstate%t(i, tropLev(i))
1053 915634 : f2 = pstate%t(i, tropLev(i)+1)
1054 :
1055 915634 : x0 = pstate%zm(i, tropLev(i)-1)
1056 915634 : x1 = pstate%zm(i, tropLev(i))
1057 915634 : x2 = pstate%zm(i, tropLev(i)+1)
1058 :
1059 915634 : c0 = (x0-x1)*(x0-x2)
1060 915634 : c1 = (x1-x0)*(x1-x2)
1061 915634 : c2 = (x2-x0)*(x2-x1)
1062 :
1063 : ! Determine the quadratic coefficients of:
1064 : ! T = a * z^2 - b*z + c
1065 915634 : a = (f0/c0 + f1/c1 + f2/c2)
1066 915634 : b = (f0/c0*(x1+x2) + f1/c1*(x0+x2) + f2/c2*(x0+x1))
1067 915634 : c = f0/c0*x1*x2 + f1/c1*x0*x2 + f2/c2*x0*x1
1068 :
1069 : ! Find the altitude of the minimum temperature
1070 915634 : tZ = 0.5_r8 * b / a
1071 :
1072 : ! The fit should be between the upper and lower points,
1073 : ! so skip the point if the fit fails.
1074 915634 : if ((tZ >= x0) .or. (tZ <= x2)) then
1075 0 : tropLev(i) = NOTFOUND
1076 : else
1077 :
1078 : ! Return the optional outputs
1079 915634 : if (present(tropP)) then
1080 915634 : tropP(i) = tropopause_interpolateP(pstate, i, tropLev(i), tZ)
1081 : end if
1082 :
1083 915634 : if (present(tropT)) then
1084 915634 : tropT(i) = a * tZ*tZ - b*tZ + c
1085 : end if
1086 :
1087 915634 : if (present(tropZ)) then
1088 915634 : tropZ(i) = tZ
1089 : end if
1090 : end if
1091 : end if
1092 : end if
1093 : end if
1094 : end do
1095 :
1096 58824 : return
1097 : end subroutine tropopause_cpp
1098 :
1099 :
1100 : ! Searches all the columns in the chunk and attempts to identify the tropopause.
1101 : ! Two routines can be specifed, a primary routine which is tried first and a
1102 : ! backup routine which will be tried only if the first routine fails. If the
1103 : ! tropopause can not be identified by either routine, then a NOTFOUND is returned
1104 : ! for the tropopause level, temperature and pressure.
1105 774000 : subroutine tropopause_find(pstate, tropLev, tropP, tropT, tropZ, primary, backup)
1106 :
1107 : implicit none
1108 :
1109 : type(physics_state), intent(in) :: pstate
1110 : integer, optional, intent(in) :: primary ! primary detection algorithm
1111 : integer, optional, intent(in) :: backup ! backup detection algorithm
1112 : integer, intent(out) :: tropLev(pcols) ! tropopause level index
1113 : real(r8), optional, intent(out) :: tropP(pcols) ! tropopause pressure (Pa)
1114 : real(r8), optional, intent(out) :: tropT(pcols) ! tropopause temperature (K)
1115 : real(r8), optional, intent(out) :: tropZ(pcols) ! tropopause height (m)
1116 :
1117 : ! Local Variable
1118 : integer :: primAlg ! Primary algorithm
1119 : integer :: backAlg ! Backup algorithm
1120 :
1121 : ! Initialize the results to a missing value, so that the algorithms will
1122 : ! attempt to find the tropopause for all of them.
1123 13158000 : tropLev(:) = NOTFOUND
1124 3774024 : if (present(tropP)) tropP(:) = fillvalue
1125 3774024 : if (present(tropT)) tropT(:) = fillvalue
1126 3774024 : if (present(tropZ)) tropZ(:) = fillvalue
1127 :
1128 : ! Set the algorithms to be used, either the ones provided or the defaults.
1129 774000 : if (present(primary)) then
1130 58824 : primAlg = primary
1131 : else
1132 715176 : primAlg = default_primary
1133 : end if
1134 :
1135 774000 : if (present(backup)) then
1136 625392 : backAlg = backup
1137 : else
1138 148608 : backAlg = default_backup
1139 : end if
1140 :
1141 : ! Try to find the tropopause using the primary algorithm.
1142 774000 : if (primAlg /= TROP_ALG_NONE) then
1143 774000 : call tropopause_findUsing(pstate, primAlg, tropLev, tropP, tropT, tropZ)
1144 : end if
1145 :
1146 12772805 : if ((backAlg /= TROP_ALG_NONE) .and. any(tropLev(:) == NOTFOUND)) then
1147 38922 : call tropopause_findUsing(pstate, backAlg, tropLev, tropP, tropT, tropZ)
1148 : end if
1149 :
1150 774000 : return
1151 : end subroutine tropopause_find
1152 :
1153 : ! Searches all the columns in the chunk and attempts to identify the "chemical"
1154 : ! tropopause. This is the lapse rate tropopause, backed up by the climatology
1155 : ! if the lapse rate fails to find the tropopause at pressures higher than a certain
1156 : ! threshold. This pressure threshold depends on latitude. Between 50S and 50N,
1157 : ! the climatology is used if the lapse rate tropopause is not found at P > 75 hPa.
1158 : ! At high latitude (poleward of 50), the threshold is increased to 125 hPa to
1159 : ! eliminate false events that are sometimes detected in the cold polar stratosphere.
1160 : !
1161 : ! NOTE: This routine was adapted from code in chemistry.F90 and mo_gasphase_chemdr.F90.
1162 507744 : subroutine tropopause_findChemTrop(pstate, tropLev, primary, backup)
1163 :
1164 : implicit none
1165 :
1166 : type(physics_state), intent(in) :: pstate
1167 : integer, optional, intent(in) :: primary ! primary detection algorithm
1168 : integer, optional, intent(in) :: backup ! backup detection algorithm
1169 : integer, intent(out) :: tropLev(pcols) ! tropopause level index
1170 :
1171 : ! Local Variable
1172 : real(r8), parameter :: rad2deg = 180._r8/pi ! radians to degrees conversion factor
1173 : real(r8) :: dlats(pcols)
1174 : integer :: i
1175 : integer :: ncol
1176 : integer :: backAlg
1177 :
1178 : ! First use the lapse rate tropopause.
1179 507744 : ncol = pstate%ncol
1180 507744 : call tropopause_find(pstate, tropLev, primary=primary, backup=TROP_ALG_NONE)
1181 :
1182 : ! Now check high latitudes (poleward of 50) and set the level to the
1183 : ! climatology if the level was not found or is at P <= 125 hPa.
1184 8478144 : dlats(:ncol) = pstate%lat(:ncol) * rad2deg ! convert to degrees
1185 :
1186 507744 : if (present(backup)) then
1187 0 : backAlg = backup
1188 : else
1189 507744 : backAlg = default_backup
1190 : end if
1191 :
1192 8478144 : do i = 1, ncol
1193 8478144 : if (abs(dlats(i)) > 50._r8) then
1194 1788256 : if (tropLev(i) .ne. NOTFOUND) then
1195 1788256 : if (pstate%pmid(i, tropLev(i)) <= 12500._r8) then
1196 0 : tropLev(i) = NOTFOUND
1197 : end if
1198 : end if
1199 : end if
1200 : end do
1201 :
1202 : ! Now use the backup algorithm
1203 8380036 : if ((backAlg /= TROP_ALG_NONE) .and. any(tropLev(:) == NOTFOUND)) then
1204 133061 : call tropopause_findUsing(pstate, backAlg, tropLev)
1205 : end if
1206 :
1207 507744 : return
1208 : end subroutine tropopause_findChemTrop
1209 :
1210 :
1211 : ! Call the appropriate tropopause detection routine based upon the algorithm
1212 : ! specifed.
1213 : !
1214 : ! NOTE: It is assumed that the output fields have been initialized by the
1215 : ! caller, and only output values set to fillvalue will be detected.
1216 945983 : subroutine tropopause_findUsing(pstate, algorithm, tropLev, tropP, tropT, tropZ)
1217 :
1218 : implicit none
1219 :
1220 : type(physics_state), intent(in) :: pstate
1221 : integer, intent(in) :: algorithm ! detection algorithm
1222 : integer, intent(inout) :: tropLev(pcols) ! tropopause level index
1223 : real(r8), optional, intent(inout) :: tropP(pcols) ! tropopause pressure (Pa)
1224 : real(r8), optional, intent(inout) :: tropT(pcols) ! tropopause temperature (K)
1225 : real(r8), optional, intent(inout) :: tropZ(pcols) ! tropopause height (m)
1226 :
1227 : ! Dispatch the request to the appropriate routine.
1228 945983 : select case(algorithm)
1229 : case(TROP_ALG_ANALYTIC)
1230 0 : call tropopause_analytic(pstate, tropLev, tropP, tropT, tropZ)
1231 :
1232 : case(TROP_ALG_CLIMATE)
1233 171983 : call tropopause_climate(pstate, tropLev, tropP, tropT, tropZ)
1234 :
1235 : case(TROP_ALG_STOBIE)
1236 0 : call tropopause_stobie(pstate, tropLev, tropP, tropT, tropZ)
1237 :
1238 : case(TROP_ALG_HYBSTOB)
1239 0 : call tropopause_hybridstobie(pstate, tropLev, tropP, tropT, tropZ)
1240 :
1241 : case(TROP_ALG_TWMO)
1242 715176 : call tropopause_twmo(pstate, tropLev, tropP, tropT, tropZ)
1243 :
1244 : case(TROP_ALG_WMO)
1245 0 : call tropopause_wmo(pstate, tropLev, tropP, tropT, tropZ)
1246 :
1247 : case(TROP_ALG_CPP)
1248 58824 : call tropopause_cpp(pstate, tropLev, tropP, tropT, tropZ)
1249 :
1250 : case default
1251 0 : write(iulog, *) 'tropopause: Invalid detection algorithm (', algorithm, ') specified.'
1252 945983 : call endrun
1253 : end select
1254 :
1255 945983 : return
1256 : end subroutine tropopause_findUsing
1257 :
1258 :
1259 : ! This routine interpolates the pressures in the physics state to
1260 : ! find the pressure at the specified tropopause altitude.
1261 915634 : function tropopause_interpolateP(pstate, icol, tropLev, tropZ)
1262 :
1263 : implicit none
1264 :
1265 : type(physics_state), intent(in) :: pstate
1266 : integer, intent(in) :: icol ! column being processed
1267 : integer, intent(in) :: tropLev ! tropopause level index
1268 : real(r8), optional, intent(in) :: tropZ ! tropopause pressure (m)
1269 : real(r8) :: tropopause_interpolateP
1270 :
1271 : ! Local Variables
1272 : real(r8) :: tropP ! tropopause pressure (Pa)
1273 : real(r8) :: dlogPdZ ! dlog(p)/dZ
1274 :
1275 : ! Interpolate the temperature linearly against log(P)
1276 :
1277 : ! Is the tropopause at the midpoint?
1278 915634 : if (tropZ == pstate%zm(icol, tropLev)) then
1279 0 : tropP = pstate%pmid(icol, tropLev)
1280 :
1281 915634 : else if (tropZ > pstate%zm(icol, tropLev)) then
1282 :
1283 : ! It is above the midpoint? Make sure we aren't at the top.
1284 446662 : if (tropLev > 1) then
1285 446662 : dlogPdZ = (log(pstate%pmid(icol, tropLev)) - log(pstate%pmid(icol, tropLev - 1))) / &
1286 893324 : (pstate%zm(icol, tropLev) - pstate%zm(icol, tropLev - 1))
1287 446662 : tropP = pstate%pmid(icol, tropLev) + exp((tropZ - pstate%zm(icol, tropLev)) * dlogPdZ)
1288 : end if
1289 : else
1290 :
1291 : ! It is below the midpoint. Make sure we aren't at the bottom.
1292 468972 : if (tropLev < pver) then
1293 468972 : dlogPdZ = (log(pstate%pmid(icol, tropLev + 1)) - log(pstate%pmid(icol, tropLev))) / &
1294 937944 : (pstate%zm(icol, tropLev + 1) - pstate%zm(icol, tropLev))
1295 468972 : tropP = pstate%pmid(icol, tropLev) + exp((tropZ - pstate%zm(icol, tropLev)) * dlogPdZ)
1296 : end if
1297 : end if
1298 :
1299 915634 : tropopause_interpolateP = tropP
1300 915634 : end function tropopause_interpolateP
1301 :
1302 :
1303 : ! This routine interpolates the temperatures in the physics state to
1304 : ! find the temperature at the specified tropopause pressure.
1305 1841751 : function tropopause_interpolateT(pstate, icol, tropLev, tropP)
1306 :
1307 : implicit none
1308 :
1309 : type(physics_state), intent(in) :: pstate
1310 : integer, intent(in) :: icol ! column being processed
1311 : integer, intent(in) :: tropLev ! tropopause level index
1312 : real(r8), optional, intent(in) :: tropP ! tropopause pressure (Pa)
1313 : real(r8) :: tropopause_interpolateT
1314 :
1315 : ! Local Variables
1316 : real(r8) :: tropT ! tropopause temperature (K)
1317 : real(r8) :: dTdlogP ! dT/dlog(P)
1318 :
1319 : ! Interpolate the temperature linearly against log(P)
1320 :
1321 : ! Is the tropopause at the midpoint?
1322 1841751 : if (tropP == pstate%pmid(icol, tropLev)) then
1323 0 : tropT = pstate%t(icol, tropLev)
1324 :
1325 1841751 : else if (tropP < pstate%pmid(icol, tropLev)) then
1326 :
1327 : ! It is above the midpoint? Make sure we aren't at the top.
1328 887910 : if (tropLev > 1) then
1329 887910 : dTdlogP = (pstate%t(icol, tropLev) - pstate%t(icol, tropLev - 1)) / &
1330 1775820 : (log(pstate%pmid(icol, tropLev)) - log(pstate%pmid(icol, tropLev - 1)))
1331 887910 : tropT = pstate%t(icol, tropLev) + (log(tropP) - log(pstate%pmid(icol, tropLev))) * dTdlogP
1332 : end if
1333 : else
1334 :
1335 : ! It is below the midpoint. Make sure we aren't at the bottom.
1336 953841 : if (tropLev < pver) then
1337 953841 : dTdlogP = (pstate%t(icol, tropLev + 1) - pstate%t(icol, tropLev)) / &
1338 1907682 : (log(pstate%pmid(icol, tropLev + 1)) - log(pstate%pmid(icol, tropLev)))
1339 953841 : tropT = pstate%t(icol, tropLev) + (log(tropP) - log(pstate%pmid(icol, tropLev))) * dTdlogP
1340 : end if
1341 : end if
1342 :
1343 1841751 : tropopause_interpolateT = tropT
1344 1841751 : end function tropopause_interpolateT
1345 :
1346 :
1347 : ! This routine interpolates the geopotential height in the physics state to
1348 : ! find the geopotential height at the specified tropopause pressure.
1349 1841751 : function tropopause_interpolateZ(pstate, icol, tropLev, tropP)
1350 : use physconst, only: rga
1351 :
1352 : implicit none
1353 :
1354 : type(physics_state), intent(in) :: pstate
1355 : integer, intent(in) :: icol ! column being processed
1356 : integer, intent(in) :: tropLev ! tropopause level index
1357 : real(r8), optional, intent(in) :: tropP ! tropopause pressure (Pa)
1358 : real(r8) :: tropopause_interpolateZ
1359 :
1360 : ! Local Variables
1361 : real(r8) :: tropZ ! tropopause geopotential height (m)
1362 : real(r8) :: dZdlogP ! dZ/dlog(P)
1363 :
1364 : ! Interpolate the geopotential height linearly against log(P)
1365 :
1366 : ! Is the tropoause at the midpoint?
1367 1841751 : if (tropP == pstate%pmid(icol, tropLev)) then
1368 0 : tropZ = pstate%zm(icol, tropLev)
1369 :
1370 1841751 : else if (tropP < pstate%pmid(icol, tropLev)) then
1371 :
1372 : ! It is above the midpoint? Make sure we aren't at the top.
1373 0 : dZdlogP = (pstate%zm(icol, tropLev) - pstate%zi(icol, tropLev)) / &
1374 887910 : (log(pstate%pmid(icol, tropLev)) - log(pstate%pint(icol, tropLev)))
1375 887910 : tropZ = pstate%zm(icol, tropLev) + (log(tropP) - log(pstate%pmid(icol, tropLev))) * dZdlogP
1376 : else
1377 :
1378 : ! It is below the midpoint. Make sure we aren't at the bottom.
1379 0 : dZdlogP = (pstate%zm(icol, tropLev) - pstate%zi(icol, tropLev+1)) / &
1380 953841 : (log(pstate%pmid(icol, tropLev)) - log(pstate%pint(icol, tropLev+1)))
1381 953841 : tropZ = pstate%zm(icol, tropLev) + (log(tropP) - log(pstate%pmid(icol, tropLev))) * dZdlogP
1382 : end if
1383 :
1384 1841751 : tropopause_interpolateZ = tropZ + pstate%phis(icol)*rga
1385 1841751 : end function tropopause_interpolateZ
1386 :
1387 :
1388 : ! Output the tropopause pressure and temperature to the history files. Two sets
1389 : ! of output will be generated, one for the default algorithm and another one
1390 : ! using the default routine, but backed by a climatology when the default
1391 : ! algorithm fails.
1392 58824 : subroutine tropopause_output(pstate)
1393 : use cam_history, only : outfld
1394 :
1395 : implicit none
1396 :
1397 : type(physics_state), intent(in) :: pstate
1398 :
1399 : ! Local Variables
1400 : integer :: i
1401 : integer :: alg
1402 : integer :: ncol ! number of cloumns in the chunk
1403 : integer :: lchnk ! chunk identifier
1404 : integer :: tropLev(pcols) ! tropopause level index
1405 : real(r8) :: tropP(pcols) ! tropopause pressure (Pa)
1406 : real(r8) :: tropT(pcols) ! tropopause temperature (K)
1407 : real(r8) :: tropZ(pcols) ! tropopause height (m)
1408 : real(r8) :: tropFound(pcols) ! tropopause found
1409 : real(r8) :: tropDZ(pcols, pver) ! relative tropopause height (m)
1410 : real(r8) :: tropPdf(pcols, pver) ! tropopause probability distribution
1411 :
1412 : ! Information about the chunk.
1413 58824 : lchnk = pstate%lchnk
1414 58824 : ncol = pstate%ncol
1415 :
1416 : ! Find the tropopause using the default algorithm backed by the climatology.
1417 58824 : call tropopause_find(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ)
1418 :
1419 58824 : tropPdf(:,:) = 0._r8
1420 58824 : tropFound(:) = 0._r8
1421 93059568 : tropDZ(:,:) = fillvalue
1422 982224 : do i = 1, ncol
1423 982224 : if (tropLev(i) /= NOTFOUND) then
1424 923400 : tropPdf(i, tropLev(i)) = 1._r8
1425 923400 : tropFound(i) = 1._r8
1426 86799600 : tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
1427 : end if
1428 : end do
1429 :
1430 58824 : call outfld('TROP_P', tropP(:ncol), ncol, lchnk)
1431 58824 : call outfld('TROP_T', tropT(:ncol), ncol, lchnk)
1432 58824 : call outfld('TROP_Z', tropZ(:ncol), ncol, lchnk)
1433 21474864 : call outfld('TROP_DZ', tropDZ(:ncol, :), ncol, lchnk)
1434 21474864 : call outfld('TROP_PD', tropPdf(:ncol, :), ncol, lchnk)
1435 58824 : call outfld('TROP_FD', tropFound(:ncol), ncol, lchnk)
1436 :
1437 :
1438 : ! Find the tropopause using just the primary algorithm.
1439 58824 : call tropopause_find(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, backup=TROP_ALG_NONE)
1440 :
1441 58824 : tropPdf(:,:) = 0._r8
1442 58824 : tropFound(:) = 0._r8
1443 93059568 : tropDZ(:,:) = fillvalue
1444 :
1445 982224 : do i = 1, ncol
1446 982224 : if (tropLev(i) /= NOTFOUND) then
1447 918351 : tropPdf(i, tropLev(i)) = 1._r8
1448 918351 : tropFound(i) = 1._r8
1449 86324994 : tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
1450 : end if
1451 : end do
1452 :
1453 58824 : call outfld('TROPP_P', tropP(:ncol), ncol, lchnk)
1454 58824 : call outfld('TROPP_T', tropT(:ncol), ncol, lchnk)
1455 58824 : call outfld('TROPP_Z', tropZ(:ncol), ncol, lchnk)
1456 21474864 : call outfld('TROPP_DZ', tropDZ(:ncol, :), ncol, lchnk)
1457 21474864 : call outfld('TROPP_PD', tropPdf(:ncol, :), ncol, lchnk)
1458 58824 : call outfld('TROPP_FD', tropFound(:ncol), ncol, lchnk)
1459 :
1460 :
1461 : ! Find the tropopause using just the cold point algorithm.
1462 58824 : call tropopause_find(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, primary=TROP_ALG_CPP, backup=TROP_ALG_NONE)
1463 :
1464 58824 : tropPdf(:,:) = 0._r8
1465 58824 : tropFound(:) = 0._r8
1466 93059568 : tropDZ(:,:) = fillvalue
1467 :
1468 982224 : do i = 1, ncol
1469 982224 : if (tropLev(i) /= NOTFOUND) then
1470 915634 : tropPdf(i, tropLev(i)) = 1._r8
1471 915634 : tropFound(i) = 1._r8
1472 86069596 : tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
1473 : end if
1474 : end do
1475 :
1476 58824 : call outfld('TROPF_P', tropP(:ncol), ncol, lchnk)
1477 58824 : call outfld('TROPF_T', tropT(:ncol), ncol, lchnk)
1478 58824 : call outfld('TROPF_Z', tropZ(:ncol), ncol, lchnk)
1479 21474864 : call outfld('TROPF_DZ', tropDZ(:ncol, :), ncol, lchnk)
1480 21474864 : call outfld('TROPF_PD', tropPdf(:ncol, :), ncol, lchnk)
1481 58824 : call outfld('TROPF_FD', tropFound(:ncol), ncol, lchnk)
1482 :
1483 :
1484 : ! If requested, do all of the algorithms.
1485 : if (output_all) then
1486 :
1487 : do alg = 2, TROP_NALG
1488 :
1489 : ! Find the tropopause using just the analytic algorithm.
1490 : call tropopause_find(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, primary=alg, backup=TROP_ALG_NONE)
1491 :
1492 : tropPdf(:,:) = 0._r8
1493 : tropFound(:) = 0._r8
1494 :
1495 : do i = 1, ncol
1496 : if (tropLev(i) /= NOTFOUND) then
1497 : tropPdf(i, tropLev(i)) = 1._r8
1498 : tropFound(i) = 1._r8
1499 : end if
1500 : end do
1501 :
1502 : call outfld('TROP' // TROP_LETTER(alg) // '_P', tropP(:ncol), ncol, lchnk)
1503 : call outfld('TROP' // TROP_LETTER(alg) // '_T', tropT(:ncol), ncol, lchnk)
1504 : call outfld('TROP' // TROP_LETTER(alg) // '_Z', tropZ(:ncol), ncol, lchnk)
1505 : call outfld('TROP' // TROP_LETTER(alg) // '_PD', tropPdf(:ncol, :), ncol, lchnk)
1506 : call outfld('TROP' // TROP_LETTER(alg) // '_FD', tropFound(:ncol), ncol, lchnk)
1507 : end do
1508 : end if
1509 :
1510 58824 : return
1511 58824 : end subroutine tropopause_output
1512 : end module tropopause
|
