Line data Source code
1 : ! This is the CAM interface to the CCPP-ized tropopause_find scheme.
2 : ! Full compatibility, bit-for-bit, to old CAM approach is achieved through
3 : ! this module, however this module will not be necessary in CAM-SIMA.
4 : !
5 : ! For science description of the underlying algorithms, refer to
6 : ! atmospheric_physics/tropopause_find/tropopause_find.F90.
7 : ! (hplin, 8/20/24)
8 :
9 : module tropopause
10 : !---------------------------------------------------------------
11 : ! ... variables for the tropopause module
12 : !---------------------------------------------------------------
13 :
14 : use shr_kind_mod, only : r8 => shr_kind_r8
15 : use shr_const_mod, only : pi => shr_const_pi
16 : use ppgrid, only : pcols, pver, pverp, begchunk, endchunk
17 : use cam_abortutils, only : endrun
18 : use cam_logfile, only : iulog
19 : use cam_history_support, only : fillvalue
20 : use physics_types, only : physics_state
21 : use spmd_utils, only : masterproc
22 :
23 : implicit none
24 :
25 : private
26 :
27 : public :: tropopause_readnl, tropopause_init, tropopause_find_cam, tropopause_output
28 : public :: tropopause_findChemTrop
29 : public :: TROP_ALG_NONE, TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
30 : public :: TROP_ALG_STOBIE, TROP_ALG_HYBSTOB, TROP_ALG_TWMO, TROP_ALG_WMO
31 : public :: TROP_ALG_CPP
32 : public :: NOTFOUND
33 :
34 : save
35 :
36 : ! These parameters define and enumeration to be used to define the primary
37 : ! and backup algorithms to be used with the tropopause_find() method. The
38 : ! backup algorithm is meant to provide a solution when the primary algorithm
39 : ! fail. The algorithms that can't fail are: TROP_ALG_ANALYTIC, TROP_ALG_CLIMATE
40 : ! and TROP_ALG_STOBIE.
41 : integer, parameter :: TROP_ALG_NONE = 1 ! Don't evaluate
42 : integer, parameter :: TROP_ALG_ANALYTIC = 2 ! Analytic Expression
43 : integer, parameter :: TROP_ALG_CLIMATE = 3 ! Climatology
44 : integer, parameter :: TROP_ALG_STOBIE = 4 ! Stobie Algorithm
45 : integer, parameter :: TROP_ALG_TWMO = 5 ! WMO Definition, Reichler et al. [2003]
46 : integer, parameter :: TROP_ALG_WMO = 6 ! WMO Definition
47 : integer, parameter :: TROP_ALG_HYBSTOB = 7 ! Hybrid Stobie Algorithm
48 : integer, parameter :: TROP_ALG_CPP = 8 ! Cold Point Parabolic
49 : integer, parameter :: TROP_ALG_CHEMTROP = 9 ! Chemical tropopause
50 :
51 : ! Note: exclude CHEMTROP here as it is a new flag added in CCPP-ized routines to unify the chemTrop routine. (hplin, 8/20/24)
52 : integer, parameter :: TROP_NALG = 8 ! Number of Algorithms
53 : character,parameter :: TROP_LETTER(TROP_NALG) = (/ ' ', 'A', 'C', 'S', 'T', 'W', 'H', 'F' /)
54 : ! unique identifier for output, don't use P
55 :
56 : ! These variables should probably be controlled by namelist entries.
57 : logical ,parameter :: output_all = .False. ! output tropopause info from all algorithms
58 : integer ,parameter :: default_primary = TROP_ALG_TWMO ! default primary algorithm
59 : integer ,parameter :: default_backup = TROP_ALG_CLIMATE ! default backup algorithm
60 :
61 : ! Namelist variables
62 : character(len=256) :: tropopause_climo_file = 'trop_climo' ! absolute filepath of climatology file
63 :
64 : ! These variables are used to store the climatology data.
65 : real(r8) :: days(12) ! days in the climatology
66 : real(r8), pointer :: tropp_p_loc(:,:,:) ! climatological tropopause pressures
67 :
68 : integer, parameter :: NOTFOUND = -1
69 :
70 : ! physical constants
71 : ! These constants are set in module variables rather than as parameters
72 : ! to support the aquaplanet mode in which the constants have values determined
73 : ! by the experiment protocol
74 : real(r8) :: cnst_kap ! = cappa
75 : real(r8) :: cnst_faktor ! = -gravit/rair
76 : real(r8) :: cnst_ka1 ! = cnst_kap - 1._r8
77 :
78 : !================================================================================================
79 : contains
80 : !================================================================================================
81 :
82 : ! Read namelist variables.
83 1490712 : subroutine tropopause_readnl(nlfile)
84 :
85 : use namelist_utils, only: find_group_name
86 : use units, only: getunit, freeunit
87 : use mpishorthand
88 :
89 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
90 :
91 : ! Local variables
92 : integer :: unitn, ierr
93 : character(len=*), parameter :: subname = 'tropopause_readnl'
94 :
95 : namelist /tropopause_nl/ tropopause_climo_file
96 : !-----------------------------------------------------------------------------
97 :
98 1536 : if (masterproc) then
99 2 : unitn = getunit()
100 2 : open( unitn, file=trim(nlfile), status='old' )
101 2 : call find_group_name(unitn, 'tropopause_nl', status=ierr)
102 2 : if (ierr == 0) then
103 2 : read(unitn, tropopause_nl, iostat=ierr)
104 2 : if (ierr /= 0) then
105 0 : call endrun(subname // ':: ERROR reading namelist')
106 : end if
107 : end if
108 2 : close(unitn)
109 2 : call freeunit(unitn)
110 : end if
111 :
112 : #ifdef SPMD
113 : ! Broadcast namelist variables
114 1536 : call mpibcast(tropopause_climo_file, len(tropopause_climo_file), mpichar, 0, mpicom)
115 : #endif
116 :
117 1536 : end subroutine tropopause_readnl
118 :
119 :
120 : ! This routine is called during intialization and must be called before the
121 : ! other methods in this module can be used. Its main tasks are to read in the
122 : ! climatology from a file and to define the output fields. Much of this code
123 : ! is taken from mo_tropopause.
124 1536 : subroutine tropopause_init()
125 :
126 : use cam_history, only: addfld, horiz_only
127 : use tropopause_find, only: tropopause_find_init
128 : use physconst, only: cappa, rair, gravit, pi
129 :
130 : character(len=512) :: errmsg
131 : integer :: errflg
132 :
133 : ! Call underlying CCPP-initialization routine.
134 1536 : call tropopause_find_init(cappa, rair, gravit, pi, errmsg, errflg)
135 :
136 : ! Define the output fields.
137 1536 : call addfld('TROP_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure', flag_xyfill=.True.)
138 1536 : call addfld('TROP_T', horiz_only, 'A', 'K', 'Tropopause Temperature', flag_xyfill=.True.)
139 1536 : call addfld('TROP_Z', horiz_only, 'A', 'm', 'Tropopause Height', flag_xyfill=.True.)
140 3072 : call addfld('TROP_DZ', (/ 'lev' /), 'A', 'm', 'Relative Tropopause Height')
141 3072 : call addfld('TROP_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Probabilty')
142 1536 : call addfld('TROP_FD', horiz_only, 'A', 'probability', 'Tropopause Found')
143 :
144 1536 : call addfld('TROPP_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (primary)', flag_xyfill=.True.)
145 1536 : call addfld('TROPP_T', horiz_only, 'A', 'K', 'Tropopause Temperature (primary)', flag_xyfill=.True.)
146 1536 : call addfld('TROPP_Z', horiz_only, 'A', 'm', 'Tropopause Height (primary)', flag_xyfill=.True.)
147 3072 : call addfld('TROPP_DZ', (/ 'lev' /), 'A', 'm', 'Relative Tropopause Height (primary)')
148 3072 : call addfld('TROPP_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (primary)')
149 1536 : call addfld('TROPP_FD', horiz_only, 'A', 'probability', 'Tropopause Found (primary)')
150 :
151 1536 : call addfld('TROPF_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (cold point)', flag_xyfill=.True.)
152 1536 : call addfld('TROPF_T', horiz_only, 'A', 'K', 'Tropopause Temperature (cold point)', flag_xyfill=.True.)
153 1536 : call addfld('TROPF_Z', horiz_only, 'A', 'm', 'Tropopause Height (cold point)', flag_xyfill=.True.)
154 3072 : call addfld('TROPF_DZ', (/ 'lev' /), 'A', 'm', 'Relative Tropopause Height (cold point)', flag_xyfill=.True.)
155 3072 : call addfld('TROPF_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (cold point)')
156 1536 : call addfld('TROPF_FD', horiz_only, 'A', 'probability', 'Tropopause Found (cold point)')
157 :
158 3072 : call addfld( 'hstobie_trop', (/ 'lev' /), 'I', 'fraction of model time', 'Lowest level with stratospheric chemsitry')
159 3072 : call addfld( 'hstobie_linoz', (/ 'lev' /), 'I', 'fraction of model time', 'Lowest possible Linoz level')
160 : call addfld( 'hstobie_tropop', (/ 'lev' /), 'I', 'fraction of model time', &
161 3072 : 'Troposphere boundary calculated in chemistry' )
162 :
163 : ! If requested, be prepared to output results from all of the methods.
164 : if (output_all) then
165 : call addfld('TROPA_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (analytic)', flag_xyfill=.True.)
166 : call addfld('TROPA_T', horiz_only, 'A', 'K', 'Tropopause Temperature (analytic)', flag_xyfill=.True.)
167 : call addfld('TROPA_Z', horiz_only, 'A', 'm', 'Tropopause Height (analytic)', flag_xyfill=.True.)
168 : call addfld('TROPA_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (analytic)')
169 : call addfld('TROPA_FD', horiz_only, 'A', 'probability', 'Tropopause Found (analytic)')
170 :
171 : call addfld('TROPC_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (climatology)', flag_xyfill=.True.)
172 : call addfld('TROPC_T', horiz_only, 'A', 'K', 'Tropopause Temperature (climatology)', flag_xyfill=.True.)
173 : call addfld('TROPC_Z', horiz_only, 'A', 'm', 'Tropopause Height (climatology)', flag_xyfill=.True.)
174 : call addfld('TROPC_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (climatology)')
175 : call addfld('TROPC_FD', horiz_only, 'A', 'probability', 'Tropopause Found (climatology)')
176 :
177 : call addfld('TROPS_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (stobie)', flag_xyfill=.True.)
178 : call addfld('TROPS_T', horiz_only, 'A', 'K', 'Tropopause Temperature (stobie)', flag_xyfill=.True.)
179 : call addfld('TROPS_Z', horiz_only, 'A', 'm', 'Tropopause Height (stobie)', flag_xyfill=.True.)
180 : call addfld('TROPS_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (stobie)')
181 : call addfld('TROPS_FD', horiz_only, 'A', 'probability', 'Tropopause Found (stobie)')
182 :
183 : call addfld('TROPT_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (twmo)', flag_xyfill=.True.)
184 : call addfld('TROPT_T', horiz_only, 'A', 'K', 'Tropopause Temperature (twmo)', flag_xyfill=.True.)
185 : call addfld('TROPT_Z', horiz_only, 'A', 'm', 'Tropopause Height (twmo)', flag_xyfill=.True.)
186 : call addfld('TROPT_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (twmo)')
187 : call addfld('TROPT_FD', horiz_only, 'A', 'probability', 'Tropopause Found (twmo)')
188 :
189 : call addfld('TROPW_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (WMO)', flag_xyfill=.True.)
190 : call addfld('TROPW_T', horiz_only, 'A', 'K', 'Tropopause Temperature (WMO)', flag_xyfill=.True.)
191 : call addfld('TROPW_Z', horiz_only, 'A', 'm', 'Tropopause Height (WMO)', flag_xyfill=.True.)
192 : call addfld('TROPW_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (WMO)')
193 : call addfld('TROPW_FD', horiz_only, 'A', 'probability', 'Tropopause Found (WMO)')
194 :
195 : call addfld('TROPH_P', horiz_only, 'A', 'Pa', 'Tropopause Pressure (Hybrid Stobie)', flag_xyfill=.True.)
196 : call addfld('TROPH_T', horiz_only, 'A', 'K', 'Tropopause Temperature (Hybrid Stobie)', flag_xyfill=.True.)
197 : call addfld('TROPH_Z', horiz_only, 'A', 'm', 'Tropopause Height (Hybrid Stobie)', flag_xyfill=.True.)
198 : call addfld('TROPH_PD', (/ 'lev' /), 'A', 'probability', 'Tropopause Distribution (Hybrid Stobie)')
199 : call addfld('TROPH_FD', horiz_only, 'A', 'probability', 'Tropopause Found (Hybrid Stobie)')
200 : end if
201 :
202 :
203 1536 : call tropopause_read_file()
204 :
205 :
206 1536 : end subroutine tropopause_init
207 :
208 :
209 1536 : subroutine tropopause_read_file
210 : !------------------------------------------------------------------
211 : ! ... initialize upper boundary values
212 : !------------------------------------------------------------------
213 1536 : use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
214 : use dyn_grid, only : get_dyn_grid_parm
215 : use phys_grid, only : get_ncols_p, get_rlat_all_p, get_rlon_all_p
216 : use ioFileMod, only : getfil
217 : use time_manager, only : get_calday
218 : use physconst, only : pi
219 : use cam_pio_utils, only: cam_pio_openfile
220 : use pio, only : file_desc_t, var_desc_t, pio_inq_dimid, pio_inq_dimlen, &
221 : pio_inq_varid, pio_get_var, pio_closefile, pio_nowrite
222 :
223 : !------------------------------------------------------------------
224 : ! ... local variables
225 : !------------------------------------------------------------------
226 : integer :: i, j, n
227 : integer :: ierr
228 : type(file_desc_t) :: pio_id
229 : integer :: dimid
230 : type(var_desc_t) :: vid
231 : integer :: nlon, nlat, ntimes
232 : integer :: start(3)
233 : integer :: count(3)
234 : integer, parameter :: dates(12) = (/ 116, 214, 316, 415, 516, 615, &
235 : 716, 816, 915, 1016, 1115, 1216 /)
236 : integer :: plon, plat
237 : type(interp_type) :: lon_wgts, lat_wgts
238 1536 : real(r8), allocatable :: tropp_p_in(:,:,:)
239 1536 : real(r8), allocatable :: lat(:)
240 1536 : real(r8), allocatable :: lon(:)
241 : real(r8) :: to_lats(pcols), to_lons(pcols)
242 : real(r8), parameter :: d2r=pi/180._r8, zero=0._r8, twopi=pi*2._r8
243 : character(len=256) :: locfn
244 : integer :: c, ncols
245 :
246 :
247 : plon = get_dyn_grid_parm('plon')
248 : plat = get_dyn_grid_parm('plat')
249 :
250 :
251 : !-----------------------------------------------------------------------
252 : ! ... open netcdf file
253 : !-----------------------------------------------------------------------
254 1536 : call getfil (tropopause_climo_file, locfn, 0)
255 1536 : call cam_pio_openfile(pio_id, trim(locfn), PIO_NOWRITE)
256 :
257 : !-----------------------------------------------------------------------
258 : ! ... get time dimension
259 : !-----------------------------------------------------------------------
260 1536 : ierr = pio_inq_dimid( pio_id, 'time', dimid )
261 1536 : ierr = pio_inq_dimlen( pio_id, dimid, ntimes )
262 1536 : if( ntimes /= 12 )then
263 0 : write(iulog,*) 'tropopause_init: number of months = ',ntimes,'; expecting 12'
264 0 : call endrun
265 : end if
266 : !-----------------------------------------------------------------------
267 : ! ... get latitudes
268 : !-----------------------------------------------------------------------
269 1536 : ierr = pio_inq_dimid( pio_id, 'lat', dimid )
270 1536 : ierr = pio_inq_dimlen( pio_id, dimid, nlat )
271 4608 : allocate( lat(nlat), stat=ierr )
272 1536 : if( ierr /= 0 ) then
273 0 : write(iulog,*) 'tropopause_init: lat allocation error = ',ierr
274 0 : call endrun
275 : end if
276 1536 : ierr = pio_inq_varid( pio_id, 'lat', vid )
277 1536 : ierr = pio_get_var( pio_id, vid, lat )
278 115200 : lat(:nlat) = lat(:nlat) * d2r
279 : !-----------------------------------------------------------------------
280 : ! ... get longitudes
281 : !-----------------------------------------------------------------------
282 1536 : ierr = pio_inq_dimid( pio_id, 'lon', dimid )
283 1536 : ierr = pio_inq_dimlen( pio_id, dimid, nlon )
284 4608 : allocate( lon(nlon), stat=ierr )
285 1536 : if( ierr /= 0 ) then
286 0 : write(iulog,*) 'tropopause_init: lon allocation error = ',ierr
287 0 : call endrun
288 : end if
289 1536 : ierr = pio_inq_varid( pio_id, 'lon', vid )
290 1536 : ierr = pio_get_var( pio_id, vid, lon )
291 224256 : lon(:nlon) = lon(:nlon) * d2r
292 :
293 : !------------------------------------------------------------------
294 : ! ... allocate arrays
295 : !------------------------------------------------------------------
296 7680 : allocate( tropp_p_in(nlon,nlat,ntimes), stat=ierr )
297 1536 : if( ierr /= 0 ) then
298 0 : write(iulog,*) 'tropopause_init: tropp_p_in allocation error = ',ierr
299 0 : call endrun
300 : end if
301 : !------------------------------------------------------------------
302 : ! ... read in the tropopause pressure
303 : !------------------------------------------------------------------
304 1536 : ierr = pio_inq_varid( pio_id, 'trop_p', vid )
305 1536 : start = (/ 1, 1, 1 /)
306 6144 : count = (/ nlon, nlat, ntimes /)
307 1536 : ierr = pio_get_var( pio_id, vid, start, count, tropp_p_in )
308 :
309 : !------------------------------------------------------------------
310 : ! ... close the netcdf file
311 : !------------------------------------------------------------------
312 1536 : call pio_closefile( pio_id )
313 :
314 : !--------------------------------------------------------------------
315 : ! ... regrid
316 : !--------------------------------------------------------------------
317 :
318 6144 : allocate( tropp_p_loc(pcols,begchunk:endchunk,ntimes), stat=ierr )
319 :
320 1536 : if( ierr /= 0 ) then
321 0 : write(iulog,*) 'tropopause_init: tropp_p_loc allocation error = ',ierr
322 0 : call endrun
323 : end if
324 :
325 7728 : do c=begchunk,endchunk
326 6192 : ncols = get_ncols_p(c)
327 6192 : call get_rlat_all_p(c, pcols, to_lats)
328 6192 : call get_rlon_all_p(c, pcols, to_lons)
329 6192 : call lininterp_init(lon, nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
330 6192 : call lininterp_init(lat, nlat, to_lats, ncols, 1, lat_wgts)
331 80496 : do n=1,ntimes
332 80496 : call lininterp(tropp_p_in(:,:,n), nlon, nlat, tropp_p_loc(1:ncols,c,n), ncols, lon_wgts, lat_wgts)
333 : end do
334 6192 : call lininterp_finish(lon_wgts)
335 7728 : call lininterp_finish(lat_wgts)
336 : end do
337 1536 : deallocate(lon)
338 1536 : deallocate(lat)
339 1536 : deallocate(tropp_p_in)
340 :
341 : !--------------------------------------------------------
342 : ! ... initialize the monthly day of year times
343 : !--------------------------------------------------------
344 :
345 19968 : do n = 1,12
346 19968 : days(n) = get_calday( dates(n), 0 )
347 : end do
348 1536 : if (masterproc) then
349 2 : write(iulog,*) 'tropopause_init : days'
350 2 : write(iulog,'(1p,5g15.8)') days(:)
351 : endif
352 :
353 3072 : end subroutine tropopause_read_file
354 :
355 : ! Searches all the columns in the chunk and attempts to identify the tropopause.
356 : ! Two routines can be specifed, a primary routine which is tried first and a
357 : ! backup routine which will be tried only if the first routine fails. If the
358 : ! tropopause can not be identified by either routine, then a NOTFOUND is returned
359 : ! for the tropopause level, temperature and pressure.
360 6702840 : subroutine tropopause_find_cam(pstate, tropLev, tropP, tropT, tropZ, primary, backup)
361 :
362 1536 : use tropopause_find, only: tropopause_findWithBackup
363 :
364 : use cam_history, only: outfld
365 : use time_manager, only: get_curr_calday
366 :
367 : implicit none
368 :
369 : type(physics_state), intent(in) :: pstate
370 : integer, optional, intent(in) :: primary ! primary detection algorithm
371 : integer, optional, intent(in) :: backup ! backup detection algorithm
372 : integer, intent(out) :: tropLev(:) ! tropopause level index
373 : real(r8), optional, intent(out) :: tropP(:) ! tropopause pressure (Pa)
374 : real(r8), optional, intent(out) :: tropT(:) ! tropopause temperature (K)
375 : real(r8), optional, intent(out) :: tropZ(:) ! tropopause height (m)
376 :
377 : ! Local Variable
378 : integer :: primAlg ! Primary algorithm
379 : integer :: backAlg ! Backup algorithm
380 :
381 : real(r8) :: calday
382 : integer :: ncol
383 :
384 : real(r8) :: hstobie_trop (pcols, pver)
385 : real(r8) :: hstobie_linoz (pcols, pver)
386 : real(r8) :: hstobie_tropop(pcols, pver)
387 :
388 : character(len=512) :: errmsg
389 : integer :: errflg
390 :
391 : ! Get compatibility variables for CCPP-ized routine
392 6702840 : ncol = pstate%ncol
393 6702840 : calday = get_curr_calday()
394 :
395 : ! Initialize the results to a missing value, so that the algorithms will
396 : ! attempt to find the tropopause for all of them. Only do this for the active columns.
397 118624680 : tropLev(:ncol) = NOTFOUND
398 81300168 : if (present(tropP)) tropP(:ncol) = fillvalue
399 81300168 : if (present(tropT)) tropT(:ncol) = fillvalue
400 81300168 : if (present(tropZ)) tropZ(:ncol) = fillvalue
401 :
402 : ! Set the algorithms to be used, either the ones provided or the defaults.
403 6702840 : if (present(primary)) then
404 1489176 : primAlg = primary
405 : else
406 5213664 : primAlg = default_primary
407 : end if
408 :
409 6702840 : if (present(backup)) then
410 2978352 : backAlg = backup
411 : else
412 3724488 : backAlg = default_backup
413 : end if
414 :
415 : ! This does not call the tropopause_find_run routine directly, because it
416 : ! computes multiple needed tropopauses simultaneously. Instead, here we
417 : ! specify the algorithm needed directly to the algorithm driver routine.
418 : call tropopause_findWithBackup( &
419 : ncol = ncol, &
420 : pver = pver, &
421 : fillvalue = fillvalue, &
422 0 : lat = pstate%lat(:ncol), &
423 0 : pint = pstate%pint(:ncol, :pverp), &
424 0 : pmid = pstate%pmid(:ncol, :pver), &
425 0 : t = pstate%t(:ncol, :pver), &
426 0 : zi = pstate%zi(:ncol, :pverp), &
427 0 : zm = pstate%zm(:ncol, :pver), &
428 0 : phis = pstate%phis(:ncol), &
429 : calday = calday, &
430 0 : tropp_p_loc = tropp_p_loc(:ncol,pstate%lchnk,:), & ! Subset into chunk as the underlying routines are no longer chunkized.
431 : tropp_days = days, &
432 6702840 : tropLev = tropLev(:ncol), &
433 : tropP = tropP, &
434 : tropT = tropT, &
435 : tropZ = tropZ, &
436 : primary = primAlg, &
437 : backup = backAlg, &
438 6702840 : hstobie_trop = hstobie_trop(:ncol, :pver), & ! Only used if TROP_ALG_HYBSTOB
439 : hstobie_linoz = hstobie_linoz(:ncol, :pver), & ! Only used if TROP_ALG_HYBSTOB
440 : hstobie_tropop = hstobie_tropop(:ncol, :pver), & ! Only used if TROP_ALG_HYBSTOB
441 : errmsg = errmsg, &
442 : errflg = errflg &
443 24579144 : )
444 :
445 : ! Output hybridstobie specific fields
446 6702840 : if(primAlg == TROP_ALG_HYBSTOB) then
447 0 : call outfld('hstobie_trop', hstobie_trop(:ncol,:), ncol, pstate%lchnk )
448 0 : call outfld('hstobie_linoz', hstobie_linoz(:ncol,:), ncol, pstate%lchnk )
449 0 : call outfld('hstobie_tropop', hstobie_tropop(:ncol,:), ncol, pstate%lchnk )
450 : endif
451 6702840 : end subroutine tropopause_find_cam
452 :
453 : ! Searches all the columns in the chunk and attempts to identify the "chemical"
454 : ! tropopause. This is the lapse rate tropopause, backed up by the climatology
455 : ! if the lapse rate fails to find the tropopause at pressures higher than a certain
456 : ! threshold. This pressure threshold depends on latitude. Between 50S and 50N,
457 : ! the climatology is used if the lapse rate tropopause is not found at P > 75 hPa.
458 : ! At high latitude (poleward of 50), the threshold is increased to 125 hPa to
459 : ! eliminate false events that are sometimes detected in the cold polar stratosphere.
460 : !
461 : ! NOTE: This routine was adapted from code in chemistry.F90 and mo_gasphase_chemdr.F90.
462 12665736 : subroutine tropopause_findChemTrop(pstate, tropLev)
463 :
464 6702840 : use tropopause_find, only: tropopause_findWithBackup
465 :
466 : use time_manager, only: get_curr_calday
467 :
468 : implicit none
469 :
470 : type(physics_state), intent(in) :: pstate
471 : integer, intent(out) :: tropLev(:) ! tropopause level index
472 :
473 : ! Local Variable
474 : real(r8) :: calday
475 : integer :: i
476 : integer :: ncol
477 :
478 : character(len=512) :: errmsg
479 : integer :: errflg
480 :
481 : ! Get compatibility variables for CCPP-ized routine
482 12665736 : ncol = pstate%ncol
483 12665736 : calday = get_curr_calday()
484 :
485 : ! Now call the unified routine with the CHEMTROP option, which has automatic
486 : ! backup fall to climatology.
487 : call tropopause_findWithBackup( &
488 : ncol = ncol, &
489 : pver = pver, &
490 : fillvalue = fillvalue, &
491 0 : lat = pstate%lat(:ncol), &
492 0 : pint = pstate%pint(:ncol, :pverp), &
493 0 : pmid = pstate%pmid(:ncol, :pver), &
494 0 : t = pstate%t(:ncol, :pver), &
495 0 : zi = pstate%zi(:ncol, :pverp), &
496 0 : zm = pstate%zm(:ncol, :pver), &
497 0 : phis = pstate%phis(:ncol), &
498 : calday = calday, &
499 0 : tropp_p_loc = tropp_p_loc(:ncol,pstate%lchnk,:), & ! Subset into chunk as the underlying routines are no longer chunkized.
500 : tropp_days = days, &
501 12665736 : tropLev = tropLev(1:ncol), &
502 : primary = TROP_ALG_CHEMTROP, &
503 : backup = TROP_ALG_CLIMATE, &
504 : errmsg = errmsg, &
505 : errflg = errflg &
506 25331472 : )
507 12665736 : end subroutine tropopause_findChemTrop
508 :
509 : ! Output the tropopause pressure and temperature to the history files. Two sets
510 : ! of output will be generated, one for the default algorithm and another one
511 : ! using the default routine, but backed by a climatology when the default
512 : ! algorithm fails.
513 1489176 : subroutine tropopause_output(pstate)
514 12665736 : use cam_history, only : outfld
515 :
516 : implicit none
517 :
518 : type(physics_state), intent(in) :: pstate
519 :
520 : ! Local Variables
521 : integer :: i
522 : integer :: alg
523 : integer :: ncol ! number of cloumns in the chunk
524 : integer :: lchnk ! chunk identifier
525 : integer :: tropLev(pcols) ! tropopause level index
526 : real(r8) :: tropP(pcols) ! tropopause pressure (Pa)
527 : real(r8) :: tropT(pcols) ! tropopause temperature (K)
528 : real(r8) :: tropZ(pcols) ! tropopause height (m)
529 : real(r8) :: tropFound(pcols) ! tropopause found
530 : real(r8) :: tropDZ(pcols, pver) ! relative tropopause height (m)
531 : real(r8) :: tropPdf(pcols, pver) ! tropopause probability distribution
532 :
533 : ! Information about the chunk.
534 1489176 : lchnk = pstate%lchnk
535 1489176 : ncol = pstate%ncol
536 :
537 : ! Find the tropopause using the default algorithm backed by the climatology.
538 1489176 : call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ)
539 :
540 1489176 : tropPdf(:,:) = 0._r8
541 1489176 : tropFound(:) = 0._r8
542 2355876432 : tropDZ(:,:) = fillvalue
543 24865776 : do i = 1, ncol
544 24865776 : if (tropLev(i) /= NOTFOUND) then
545 23376600 : tropPdf(i, tropLev(i)) = 1._r8
546 23376600 : tropFound(i) = 1._r8
547 2197400400 : tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
548 : end if
549 : end do
550 :
551 1489176 : call outfld('TROP_P', tropP(:ncol), ncol, lchnk)
552 1489176 : call outfld('TROP_T', tropT(:ncol), ncol, lchnk)
553 1489176 : call outfld('TROP_Z', tropZ(:ncol), ncol, lchnk)
554 543653136 : call outfld('TROP_DZ', tropDZ(:ncol, :), ncol, lchnk)
555 543653136 : call outfld('TROP_PD', tropPdf(:ncol, :), ncol, lchnk)
556 1489176 : call outfld('TROP_FD', tropFound(:ncol), ncol, lchnk)
557 :
558 :
559 : ! Find the tropopause using just the primary algorithm.
560 1489176 : call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, backup=TROP_ALG_NONE)
561 :
562 1489176 : tropPdf(:,:) = 0._r8
563 1489176 : tropFound(:) = 0._r8
564 2355876432 : tropDZ(:,:) = fillvalue
565 :
566 24865776 : do i = 1, ncol
567 24865776 : if (tropLev(i) /= NOTFOUND) then
568 23205946 : tropPdf(i, tropLev(i)) = 1._r8
569 23205946 : tropFound(i) = 1._r8
570 2181358924 : tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
571 : end if
572 : end do
573 :
574 1489176 : call outfld('TROPP_P', tropP(:ncol), ncol, lchnk)
575 1489176 : call outfld('TROPP_T', tropT(:ncol), ncol, lchnk)
576 1489176 : call outfld('TROPP_Z', tropZ(:ncol), ncol, lchnk)
577 543653136 : call outfld('TROPP_DZ', tropDZ(:ncol, :), ncol, lchnk)
578 543653136 : call outfld('TROPP_PD', tropPdf(:ncol, :), ncol, lchnk)
579 1489176 : call outfld('TROPP_FD', tropFound(:ncol), ncol, lchnk)
580 :
581 :
582 : ! Find the tropopause using just the cold point algorithm.
583 1489176 : call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, primary=TROP_ALG_CPP, backup=TROP_ALG_NONE)
584 :
585 1489176 : tropPdf(:,:) = 0._r8
586 1489176 : tropFound(:) = 0._r8
587 2355876432 : tropDZ(:,:) = fillvalue
588 :
589 24865776 : do i = 1, ncol
590 24865776 : if (tropLev(i) /= NOTFOUND) then
591 23101516 : tropPdf(i, tropLev(i)) = 1._r8
592 23101516 : tropFound(i) = 1._r8
593 2171542504 : tropDZ(i,:) = pstate%zm(i,:) - tropZ(i)
594 : end if
595 : end do
596 :
597 1489176 : call outfld('TROPF_P', tropP(:ncol), ncol, lchnk)
598 1489176 : call outfld('TROPF_T', tropT(:ncol), ncol, lchnk)
599 1489176 : call outfld('TROPF_Z', tropZ(:ncol), ncol, lchnk)
600 543653136 : call outfld('TROPF_DZ', tropDZ(:ncol, :), ncol, lchnk)
601 543653136 : call outfld('TROPF_PD', tropPdf(:ncol, :), ncol, lchnk)
602 1489176 : call outfld('TROPF_FD', tropFound(:ncol), ncol, lchnk)
603 :
604 :
605 : ! If requested, do all of the algorithms.
606 : if (output_all) then
607 :
608 : do alg = 2, TROP_NALG
609 :
610 : ! Find the tropopause using just the analytic algorithm.
611 : call tropopause_find_cam(pstate, tropLev, tropP=tropP, tropT=tropT, tropZ=tropZ, primary=alg, backup=TROP_ALG_NONE)
612 :
613 : tropPdf(:,:) = 0._r8
614 : tropFound(:) = 0._r8
615 :
616 : do i = 1, ncol
617 : if (tropLev(i) /= NOTFOUND) then
618 : tropPdf(i, tropLev(i)) = 1._r8
619 : tropFound(i) = 1._r8
620 : end if
621 : end do
622 :
623 : call outfld('TROP' // TROP_LETTER(alg) // '_P', tropP(:ncol), ncol, lchnk)
624 : call outfld('TROP' // TROP_LETTER(alg) // '_T', tropT(:ncol), ncol, lchnk)
625 : call outfld('TROP' // TROP_LETTER(alg) // '_Z', tropZ(:ncol), ncol, lchnk)
626 : call outfld('TROP' // TROP_LETTER(alg) // '_PD', tropPdf(:ncol, :), ncol, lchnk)
627 : call outfld('TROP' // TROP_LETTER(alg) // '_FD', tropFound(:ncol), ncol, lchnk)
628 : end do
629 : end if
630 1489176 : end subroutine tropopause_output
631 : end module tropopause
|
