Line data Source code
1 : module gw_drag
2 :
3 : !--------------------------------------------------------------------------
4 : ! CAM and WACCM gravity wave parameterizations were merged by Sean Patrick
5 : ! Santos in Summer 2013, and at the same time, gw_drag was split into
6 : ! various modules. This is the CAM interface and driver module. The below
7 : ! notes are for the old CAM and WACCM versions of gw_drag.
8 : !--------------------------------------------------------------------------
9 : ! This file came from wa17 and was modified by Fabrizio: 07-02-2004
10 : ! Standard gw_drag with modification (6) of latitude profile of gw spectrum
11 : !--------------------------------------------------------------------------
12 : ! Purpose:
13 : !
14 : ! Module to compute the forcing due to parameterized gravity waves. Both an
15 : ! orographic and an internal source spectrum are considered.
16 : !
17 : ! Author: Byron Boville
18 : !
19 : !--------------------------------------------------------------------------
20 : use shr_kind_mod, only: r8=>shr_kind_r8, cl=>shr_kind_cl
21 : use shr_log_mod, only: errMsg => shr_log_errMsg
22 : use shr_assert_mod, only: shr_assert
23 :
24 : use ppgrid, only: pcols, pver, begchunk, endchunk
25 : use constituents, only: pcnst
26 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
27 : use spmd_utils, only: masterproc
28 : use cam_history, only: outfld
29 : use cam_logfile, only: iulog
30 : use cam_abortutils, only: endrun
31 : use error_messages, only: alloc_err
32 :
33 :
34 : use ref_pres, only: do_molec_diff, nbot_molec, press_lim_idx
35 : use physconst, only: cpair
36 :
37 : ! These are the actual switches for different gravity wave sources.
38 : use phys_control, only: use_gw_oro, use_gw_front, use_gw_front_igw, &
39 : use_gw_convect_dp, use_gw_convect_sh, &
40 : use_simple_phys, use_gw_movmtn_pbl
41 :
42 : use gw_common, only: GWBand
43 : use gw_convect, only: BeresSourceDesc
44 : use gw_movmtn, only: MovMtnSourceDesc
45 : use gw_front, only: CMSourceDesc
46 :
47 : ! Typical module header
48 : implicit none
49 : private
50 : save
51 :
52 : !
53 : ! PUBLIC: interfaces
54 : !
55 : public :: gw_drag_readnl ! Read namelist
56 : public :: gw_init ! Initialization
57 : public :: gw_tend ! interface to actual parameterization
58 :
59 : !
60 : ! PRIVATE: Rest of the data and interfaces are private to this module
61 : !
62 : real(r8), parameter :: unset_r8 = huge(1._r8)
63 :
64 : ! A mid-scale "band" with only stationary waves (l = 0).
65 : type(GWBand) :: band_oro
66 : ! Medium scale waves.
67 : type(GWBand) :: band_mid
68 : ! Long scale waves for IGWs.
69 : type(GWBand) :: band_long
70 : ! Medium scale waves for moving mountain
71 : type(GWBand) :: band_movmtn
72 :
73 : ! Top level for gravity waves.
74 : integer, parameter :: ktop = 1
75 : ! Bottom level for frontal waves.
76 : integer :: kbot_front
77 :
78 : ! Factor for SH orographic waves.
79 : real(r8) :: gw_oro_south_fac = 1._r8
80 :
81 : ! Frontogenesis function critical threshold.
82 : real(r8) :: frontgfc = unset_r8
83 :
84 : ! Tendency efficiencies.
85 :
86 : ! Ridge scheme.
87 : logical :: use_gw_rdg_beta = .false.
88 : integer :: n_rdg_beta = 1
89 : real(r8) :: effgw_rdg_beta = unset_r8
90 : real(r8) :: effgw_rdg_beta_max = unset_r8
91 : real(r8) :: rdg_beta_cd_llb = unset_r8 ! Low-level obstacle drag coefficient Ridge scheme.
92 : logical :: trpd_leewv_rdg_beta = .false.
93 :
94 : logical :: use_gw_rdg_gamma = .false.
95 : integer :: n_rdg_gamma = -1
96 : real(r8) :: effgw_rdg_gamma = unset_r8
97 : real(r8) :: effgw_rdg_gamma_max = unset_r8
98 : real(r8) :: rdg_gamma_cd_llb = unset_r8
99 : logical :: trpd_leewv_rdg_gamma = .false.
100 : character(len=cl) :: bnd_rdggm = 'bnd_rdggm' ! full pathname for meso-Gamma ridge dataset
101 :
102 : ! Orography.
103 : real(r8) :: effgw_oro = unset_r8
104 : ! C&M scheme.
105 : real(r8) :: effgw_cm = unset_r8
106 : ! C&M scheme (inertial waves).
107 : real(r8) :: effgw_cm_igw = unset_r8
108 : ! Beres (deep convection).
109 : real(r8) :: effgw_beres_dp = unset_r8
110 : ! Beres (shallow convection).
111 : real(r8) :: effgw_beres_sh = unset_r8
112 :
113 : ! Horzontal wavelengths [m].
114 : real(r8), parameter :: wavelength_mid = 1.e5_r8
115 : real(r8), parameter :: wavelength_long = 3.e5_r8
116 :
117 : ! Background stress source strengths.
118 : real(r8) :: taubgnd = unset_r8
119 : real(r8) :: taubgnd_igw = unset_r8
120 :
121 : ! Whether or not to use a polar taper for frontally generated waves.
122 : logical :: gw_polar_taper = .false.
123 :
124 : ! Whether or not to enforce an upper boundary condition of tau = 0.
125 : ! (Like many variables, this is only here to hold the value between
126 : ! the readnl phase and the init phase of the CAM physics; only gw_common
127 : ! should actually use it.)
128 : logical :: tau_0_ubc = .false.
129 :
130 : ! Whether or not to limit tau *before* applying any efficiency factors.
131 : logical :: gw_limit_tau_without_eff = .false.
132 :
133 : ! Whether or not to apply tendency max
134 : logical :: gw_apply_tndmax = .true.
135 :
136 : ! Files to read Beres source spectra from.
137 : character(len=cl) :: gw_drag_file = ""
138 : character(len=cl) :: gw_drag_file_sh = ""
139 : character(len=cl) :: gw_drag_file_mm = ""
140 :
141 : ! Beres settings and table.
142 : type(BeresSourceDesc) :: beres_dp_desc
143 : type(BeresSourceDesc) :: beres_sh_desc
144 :
145 : ! Moving mountain settings and table.
146 : type(MovMtnSourceDesc) :: movmtn_desc
147 :
148 : ! Frontogenesis wave settings.
149 : type(CMSourceDesc) :: cm_desc
150 : type(CMSourceDesc) :: cm_igw_desc
151 :
152 : ! Indices into pbuf
153 : integer :: kvt_idx = -1
154 : integer :: ttend_dp_idx = -1
155 : integer :: ttend_sh_idx = -1
156 : integer :: frontgf_idx = -1
157 : integer :: frontga_idx = -1
158 : integer :: sgh_idx = -1
159 :
160 : ! From CLUBB
161 : integer :: ttend_clubb_idx = -1
162 : integer :: upwp_clubb_gw_idx = -1
163 : integer :: vpwp_clubb_gw_idx = -1
164 : integer :: thlp2_clubb_gw_idx = -1
165 : integer :: wpthlp_clubb_gw_idx = -1
166 :
167 : ! anisotropic ridge fields
168 : integer, parameter :: prdg = 16
169 :
170 : real(r8), allocatable, dimension(:,:), target :: &
171 : rdg_gbxar
172 :
173 : ! Meso Beta
174 : real(r8), allocatable, dimension(:,:,:), target :: &
175 : rdg_hwdth, &
176 : rdg_clngt, &
177 : rdg_mxdis, &
178 : rdg_anixy, &
179 : rdg_angll
180 :
181 : ! Meso Gamma
182 : real(r8), allocatable, dimension(:,:,:), target :: &
183 : rdg_hwdthg, &
184 : rdg_clngtg, &
185 : rdg_mxdisg, &
186 : rdg_anixyg, &
187 : rdg_angllg
188 :
189 : ! Water constituent indices for budget
190 : integer :: ixcldliq = -1
191 : integer :: ixcldice = -1
192 :
193 : ! Prefixes for history field names
194 : character(len=1), parameter :: cm_pf = " "
195 : character(len=1), parameter :: cm_igw_pf = "I"
196 : character(len=1), parameter :: beres_dp_pf = "B"
197 : character(len=1), parameter :: beres_sh_pf = "S"
198 :
199 : ! namelist
200 : logical :: history_amwg ! output the variables used by the AMWG diag package
201 : logical :: gw_lndscl_sgh = .true. ! scale SGH by land frac
202 : real(r8) :: gw_prndl = 0.25_r8
203 : real(r8) :: gw_qbo_hdepth_scaling = 1._r8 ! heating depth scaling factor
204 :
205 : ! Width of gaussian used to create frontogenesis tau profile [m s-1].
206 : real(r8) :: front_gaussian_width = -huge(1._r8)
207 :
208 : real(r8) :: alpha_gw_movmtn
209 :
210 : logical :: gw_top_taper=.false.
211 : real(r8), pointer :: vramp(:)=>null()
212 :
213 : !==========================================================================
214 : contains
215 : !==========================================================================
216 :
217 1536 : subroutine gw_drag_readnl(nlfile)
218 :
219 : use namelist_utils, only: find_group_name
220 : use units, only: getunit, freeunit
221 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8, &
222 : mpi_character, mpi_logical, mpi_integer
223 : use gw_rdg, only: gw_rdg_readnl
224 :
225 : ! File containing namelist input.
226 : character(len=*), intent(in) :: nlfile
227 :
228 : ! Local variables
229 : integer :: unitn, ierr
230 : character(len=*), parameter :: sub = 'gw_drag_readnl'
231 :
232 : ! Maximum wave number and width of spectrum bins.
233 : integer :: pgwv = -1
234 : real(r8) :: gw_dc = unset_r8
235 : integer :: pgwv_long = -1
236 : real(r8) :: gw_dc_long = unset_r8
237 :
238 : namelist /gw_drag_nl/ pgwv, gw_dc, pgwv_long, gw_dc_long, tau_0_ubc, &
239 : effgw_beres_dp, effgw_beres_sh, effgw_cm, effgw_cm_igw, effgw_oro, &
240 : frontgfc, gw_drag_file, gw_drag_file_sh, gw_drag_file_mm, taubgnd, &
241 : taubgnd_igw, gw_polar_taper, &
242 : use_gw_rdg_beta, n_rdg_beta, effgw_rdg_beta, effgw_rdg_beta_max, &
243 : rdg_beta_cd_llb, trpd_leewv_rdg_beta, &
244 : use_gw_rdg_gamma, n_rdg_gamma, effgw_rdg_gamma, effgw_rdg_gamma_max, &
245 : rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, bnd_rdggm, &
246 : gw_oro_south_fac, gw_limit_tau_without_eff, &
247 : gw_lndscl_sgh, gw_prndl, gw_apply_tndmax, gw_qbo_hdepth_scaling, &
248 : gw_top_taper, front_gaussian_width, alpha_gw_movmtn
249 : !----------------------------------------------------------------------
250 :
251 1536 : if (use_simple_phys) return
252 :
253 1536 : if (masterproc) then
254 2 : unitn = getunit()
255 2 : open( unitn, file=trim(nlfile), status='old' )
256 2 : call find_group_name(unitn, 'gw_drag_nl', status=ierr)
257 2 : if (ierr == 0) then
258 2 : read(unitn, gw_drag_nl, iostat=ierr)
259 2 : if (ierr /= 0) then
260 0 : call endrun(sub // ':: ERROR reading namelist')
261 : end if
262 : end if
263 2 : close(unitn)
264 2 : call freeunit(unitn)
265 : end if
266 :
267 1536 : call mpi_bcast(pgwv, 1, mpi_integer, mstrid, mpicom, ierr)
268 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: pgwv")
269 1536 : call mpi_bcast(gw_dc, 1, mpi_real8, mstrid, mpicom, ierr)
270 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_dc")
271 1536 : call mpi_bcast(pgwv_long, 1, mpi_integer, mstrid, mpicom, ierr)
272 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: pgwv_long")
273 1536 : call mpi_bcast(gw_dc_long, 1, mpi_real8, mstrid, mpicom, ierr)
274 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_dc_long")
275 1536 : call mpi_bcast(tau_0_ubc, 1, mpi_logical, mstrid, mpicom, ierr)
276 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: tau_0_ubc")
277 1536 : call mpi_bcast(effgw_beres_dp, 1, mpi_real8, mstrid, mpicom, ierr)
278 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_beres_dp")
279 1536 : call mpi_bcast(effgw_beres_sh, 1, mpi_real8, mstrid, mpicom, ierr)
280 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_beres_sh")
281 1536 : call mpi_bcast(effgw_cm, 1, mpi_real8, mstrid, mpicom, ierr)
282 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_cm")
283 1536 : call mpi_bcast(effgw_cm_igw, 1, mpi_real8, mstrid, mpicom, ierr)
284 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_cm_igw")
285 1536 : call mpi_bcast(effgw_oro, 1, mpi_real8, mstrid, mpicom, ierr)
286 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_oro")
287 :
288 1536 : call mpi_bcast(use_gw_rdg_beta, 1, mpi_logical, mstrid, mpicom, ierr)
289 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_beta")
290 1536 : call mpi_bcast(n_rdg_beta, 1, mpi_integer, mstrid, mpicom, ierr)
291 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: n_rdg_beta")
292 1536 : call mpi_bcast(effgw_rdg_beta, 1, mpi_real8, mstrid, mpicom, ierr)
293 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_beta")
294 1536 : call mpi_bcast(effgw_rdg_beta_max, 1, mpi_real8, mstrid, mpicom, ierr)
295 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_beta_max")
296 1536 : call mpi_bcast(rdg_beta_cd_llb, 1, mpi_real8, mstrid, mpicom, ierr)
297 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rdg_beta_cd_llb")
298 1536 : call mpi_bcast(trpd_leewv_rdg_beta, 1, mpi_logical, mstrid, mpicom, ierr)
299 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: trpd_leewv_rdg_beta")
300 :
301 1536 : call mpi_bcast(use_gw_rdg_gamma, 1, mpi_logical, mstrid, mpicom, ierr)
302 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_gw_rdg_gamma")
303 1536 : call mpi_bcast(n_rdg_gamma, 1, mpi_integer, mstrid, mpicom, ierr)
304 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: n_rdg_gamma")
305 1536 : call mpi_bcast(effgw_rdg_gamma, 1, mpi_real8, mstrid, mpicom, ierr)
306 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_gamma")
307 1536 : call mpi_bcast(effgw_rdg_gamma_max, 1, mpi_real8, mstrid, mpicom, ierr)
308 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: effgw_rdg_gamma_max")
309 1536 : call mpi_bcast(rdg_gamma_cd_llb, 1, mpi_real8, mstrid, mpicom, ierr)
310 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rdg_gamma_cd_llb")
311 1536 : call mpi_bcast(trpd_leewv_rdg_gamma, 1, mpi_logical, mstrid, mpicom, ierr)
312 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: trpd_leewv_rdg_gamma")
313 1536 : call mpi_bcast(bnd_rdggm, len(bnd_rdggm), mpi_character, mstrid, mpicom, ierr)
314 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: bnd_rdggm")
315 :
316 1536 : call mpi_bcast(gw_oro_south_fac, 1, mpi_real8, mstrid, mpicom, ierr)
317 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_oro_south_fac")
318 1536 : call mpi_bcast(frontgfc, 1, mpi_real8, mstrid, mpicom, ierr)
319 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frontgfc")
320 1536 : call mpi_bcast(taubgnd, 1, mpi_real8, mstrid, mpicom, ierr)
321 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: taubgnd")
322 1536 : call mpi_bcast(taubgnd_igw, 1, mpi_real8, mstrid, mpicom, ierr)
323 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: taubgnd_igw")
324 :
325 1536 : call mpi_bcast(gw_polar_taper, 1, mpi_logical, mstrid, mpicom, ierr)
326 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_polar_taper")
327 1536 : call mpi_bcast(gw_limit_tau_without_eff, 1, mpi_logical, mstrid, mpicom, ierr)
328 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_limit_tau_without_eff")
329 1536 : call mpi_bcast(gw_apply_tndmax, 1, mpi_logical, mstrid, mpicom, ierr)
330 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_apply_tndmax")
331 :
332 1536 : call mpi_bcast(gw_top_taper, 1, mpi_logical, mstrid, mpicom, ierr)
333 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_top_taper")
334 :
335 1536 : call mpi_bcast(gw_lndscl_sgh, 1, mpi_logical, mstrid, mpicom, ierr)
336 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_lndscl_sgh")
337 1536 : call mpi_bcast(gw_prndl, 1, mpi_real8, mstrid, mpicom, ierr)
338 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_prndl")
339 1536 : call mpi_bcast(gw_qbo_hdepth_scaling, 1, mpi_real8, mstrid, mpicom, ierr)
340 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_qbo_hdepth_scaling")
341 :
342 1536 : call mpi_bcast(gw_drag_file, len(gw_drag_file), mpi_character, mstrid, mpicom, ierr)
343 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file")
344 1536 : call mpi_bcast(gw_drag_file_sh, len(gw_drag_file_sh), mpi_character, mstrid, mpicom, ierr)
345 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file_sh")
346 1536 : call mpi_bcast(gw_drag_file_mm, len(gw_drag_file_mm), mpi_character, mstrid, mpicom, ierr)
347 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: gw_drag_file_mm")
348 :
349 1536 : call mpi_bcast(front_gaussian_width, 1, mpi_real8, mstrid, mpicom, ierr)
350 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: front_gaussian_width")
351 :
352 1536 : call mpi_bcast(alpha_gw_movmtn, 1, mpi_real8, mstrid, mpicom, ierr)
353 1536 : if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: alpha_gw_movmtn")
354 :
355 : ! Check if pgwv was set.
356 : call shr_assert(pgwv >= 0, &
357 : "gw_drag_readnl: pgwv must be set via the namelist and &
358 : &non-negative."// &
359 1536 : errMsg(__FILE__, __LINE__))
360 :
361 : ! Check if gw_dc was set.
362 : call shr_assert(gw_dc /= unset_r8, &
363 : "gw_drag_readnl: gw_dc must be set via the namelist."// &
364 1536 : errMsg(__FILE__, __LINE__))
365 :
366 1536 : band_oro = GWBand(0, gw_dc, 1.0_r8, wavelength_mid)
367 1536 : band_mid = GWBand(pgwv, gw_dc, 1.0_r8, wavelength_mid)
368 1536 : band_long = GWBand(pgwv_long, gw_dc_long, 1.0_r8, wavelength_long)
369 1536 : band_movmtn = GWBand(0, gw_dc, 1.0_r8, wavelength_mid)
370 :
371 1536 : if (use_gw_rdg_gamma .or. use_gw_rdg_beta) then
372 0 : call gw_rdg_readnl(nlfile)
373 : end if
374 :
375 : end subroutine gw_drag_readnl
376 :
377 : !==========================================================================
378 :
379 1536 : subroutine gw_init()
380 : !-----------------------------------------------------------------------
381 : ! Time independent initialization for multiple gravity wave
382 : ! parameterization.
383 : !-----------------------------------------------------------------------
384 :
385 : use cam_history, only: addfld, add_default, horiz_only
386 : use cam_history, only: register_vector_field
387 : use interpolate_data, only: lininterp
388 : use phys_control, only: phys_getopts
389 : use physics_buffer, only: pbuf_get_index
390 : use constituents, only: cnst_get_ind
391 :
392 : use cam_initfiles, only: topo_file_get_id
393 :
394 : ! temporary for restart with ridge scheme
395 : use cam_initfiles, only: bnd_topo
396 :
397 : use cam_pio_utils, only: cam_pio_openfile
398 : use cam_grid_support, only: cam_grid_check, cam_grid_id
399 : use cam_grid_support, only: cam_grid_get_dim_names
400 : use pio, only: file_desc_t, pio_nowrite, pio_closefile
401 : use ncdio_atm, only: infld
402 : use ioFileMod, only: getfil
403 :
404 : use ref_pres, only: pref_edge, pref_mid
405 : use physconst, only: gravit, rair, rearth
406 :
407 : use gw_common, only: gw_common_init
408 : use gw_front, only: gaussian_cm_desc
409 :
410 : !---------------------------Local storage-------------------------------
411 :
412 : integer :: i, l, k
413 : character(len=1) :: cn
414 :
415 : ! Index for levels at specific pressures.
416 : integer :: kfront
417 :
418 : ! output tendencies and state variables for CAM4 temperature,
419 : ! water vapor, cloud ice and cloud liquid budgets.
420 : logical :: history_budget
421 : ! output history file number for budget fields
422 : integer :: history_budget_histfile_num
423 : ! output variables of interest in WACCM runs
424 : logical :: history_waccm
425 :
426 : ! Interpolated Newtonian cooling coefficients.
427 : real(r8) :: alpha(pver+1)
428 :
429 : ! Levels of pre-calculated Newtonian cooling (1/day).
430 : ! The following profile is digitized from:
431 : ! Wehrbein and Leovy (JAS, 39, 1532-1544, 1982) figure 5
432 :
433 : integer, parameter :: nalph = 71
434 : real(r8) :: alpha0(nalph) = [ &
435 : 0.1_r8, 0.1_r8, 0.1_r8, 0.1_r8, &
436 : 0.1_r8, 0.1_r8, 0.1_r8, 0.1_r8, &
437 : 0.1_r8, 0.1_r8, 0.10133333_r8, 0.104_r8, &
438 : 0.108_r8, 0.112_r8, 0.116_r8, 0.12066667_r8, &
439 : 0.126_r8, 0.132_r8, 0.138_r8, 0.144_r8, &
440 : 0.15133333_r8, 0.16_r8, 0.17_r8, 0.18_r8, &
441 : 0.19_r8, 0.19933333_r8, 0.208_r8, 0.216_r8, &
442 : 0.224_r8, 0.232_r8, 0.23466667_r8, 0.232_r8, &
443 : 0.224_r8, 0.216_r8, 0.208_r8, 0.20133333_r8, &
444 : 0.196_r8, 0.192_r8, 0.188_r8, 0.184_r8, &
445 : 0.18266667_r8, 0.184_r8, 0.188_r8, 0.192_r8, &
446 : 0.196_r8, 0.19333333_r8, 0.184_r8, 0.168_r8, &
447 : 0.152_r8, 0.136_r8, 0.12133333_r8, 0.108_r8, &
448 : 0.096_r8, 0.084_r8, 0.072_r8, 0.061_r8, &
449 : 0.051_r8, 0.042_r8, 0.033_r8, 0.024_r8, &
450 : 0.017666667_r8, 0.014_r8, 0.013_r8, 0.012_r8, &
451 : 0.011_r8, 0.010333333_r8, 0.01_r8, 0.01_r8, &
452 : 0.01_r8, 0.01_r8, 0.01_r8 &
453 : ]
454 :
455 : ! Pressure levels that were used to calculate alpha0 (hPa).
456 : real(r8) :: palph(nalph) = [ &
457 : 2.06115E-06_r8, 2.74280E-06_r8, 3.64988E-06_r8, 4.85694E-06_r8, &
458 : 6.46319E-06_r8, 8.60065E-06_r8, 1.14450E-05_r8, 1.52300E-05_r8, &
459 : 2.02667E-05_r8, 2.69692E-05_r8, 3.58882E-05_r8, 4.77568E-05_r8, &
460 : 6.35507E-05_r8, 8.45676E-05_r8, 0.000112535_r8, 0.000149752_r8, &
461 : 0.000199277_r8, 0.000265180_r8, 0.000352878_r8, 0.000469579_r8, &
462 : 0.000624875_r8, 0.000831529_r8, 0.00110653_r8, 0.00147247_r8, &
463 : 0.00195943_r8, 0.00260744_r8, 0.00346975_r8, 0.00461724_r8, &
464 : 0.00614421_r8, 0.00817618_r8, 0.0108801_r8, 0.0144783_r8, &
465 : 0.0192665_r8, 0.0256382_r8, 0.0341170_r8, 0.0453999_r8, &
466 : 0.0604142_r8, 0.0803939_r8, 0.106981_r8, 0.142361_r8, &
467 : 0.189442_r8, 0.252093_r8, 0.335463_r8, 0.446404_r8, &
468 : 0.594036_r8, 0.790490_r8, 1.05192_r8, 1.39980_r8, &
469 : 1.86273_r8, 2.47875_r8, 3.29851_r8, 4.38936_r8, &
470 : 5.84098_r8, 7.77266_r8, 10.3432_r8, 13.7638_r8, &
471 : 18.3156_r8, 24.3728_r8, 32.4332_r8, 43.1593_r8, &
472 : 57.4326_r8, 76.4263_r8, 101.701_r8, 135.335_r8, &
473 : 180.092_r8, 239.651_r8, 318.907_r8, 424.373_r8, &
474 : 564.718_r8, 751.477_r8, 1000._r8 &
475 : ]
476 :
477 : ! Read data from file
478 : type(file_desc_t), pointer :: fh_topo
479 : type(file_desc_t) :: fh_rdggm
480 : integer :: grid_id
481 : character(len=8) :: dim1name, dim2name
482 : logical :: found
483 : character(len=cl) :: bnd_rdggm_loc ! filepath of topo file on local disk
484 :
485 : ! Allow reporting of error messages.
486 : character(len=128) :: errstring
487 : character(len=*), parameter :: sub = 'gw_init'
488 :
489 : ! temporary workaround for restart w/ ridge scheme
490 : character(len=cl) :: bnd_topo_loc ! filepath of topo file on local disk
491 :
492 : integer :: botndx,topndx
493 :
494 : !-----------------------------------------------------------------------
495 :
496 1536 : if (do_molec_diff) then
497 0 : kvt_idx = pbuf_get_index('kvt')
498 : end if
499 :
500 1536 : if (masterproc) then
501 2 : write(iulog,*) ' '
502 2 : write(iulog,*) "GW_DRAG: band_mid%ngwv = ", band_mid%ngwv
503 4 : do l = -band_mid%ngwv, band_mid%ngwv
504 : write (iulog,'(A,I0,A,F7.2)') &
505 4 : "GW_DRAG: band_mid%cref(",l,") = ",band_mid%cref(l)
506 : enddo
507 2 : write(iulog,*) 'GW_DRAG: band_mid%kwv = ', band_mid%kwv
508 2 : write(iulog,*) 'GW_DRAG: band_mid%fcrit2 = ', band_mid%fcrit2
509 2 : write(iulog,*) ' '
510 2 : write(iulog,*) "GW_DRAG: band_long%ngwv = ", band_long%ngwv
511 4 : do l = -band_long%ngwv, band_long%ngwv
512 : write (iulog,'(A,I2,A,F7.2)') &
513 4 : "GW_DRAG: band_long%cref(",l,") = ",band_long%cref(l)
514 : enddo
515 2 : write(iulog,*) 'GW_DRAG: band_long%kwv = ', band_long%kwv
516 2 : write(iulog,*) 'GW_DRAG: band_long%fcrit2 = ', band_long%fcrit2
517 2 : write(iulog,*) ' '
518 : end if
519 :
520 : ! pre-calculated newtonian damping:
521 : ! * convert to s-1
522 : ! * ensure it is not smaller than 1e-6
523 : ! * convert palph from hpa to pa
524 :
525 110592 : do k=1,nalph
526 109056 : alpha0(k) = alpha0(k) / 86400._r8
527 109056 : alpha0(k) = max(alpha0(k), 1.e-6_r8)
528 110592 : palph(k) = palph(k)*1.e2_r8
529 : end do
530 :
531 : ! interpolate to current vertical grid and obtain alpha
532 :
533 1536 : call lininterp (alpha0 ,palph, nalph , alpha , pref_edge , pver+1)
534 1536 : if (masterproc) then
535 2 : write (iulog,*) 'gw_init: newtonian damping (1/day):'
536 2 : write (iulog,fmt='(a4,a12,a10)') ' k ',' pref_edge ', &
537 4 : ' alpha '
538 56 : do k = 1, pver+1
539 54 : write (iulog,fmt='(i4,1e12.5,1f10.2)') k,pref_edge(k), &
540 110 : alpha(k)*86400._r8
541 : end do
542 : end if
543 :
544 1536 : if (masterproc) then
545 2 : write(iulog,*) 'KTOP =',ktop
546 : end if
547 :
548 : ! Used to decide whether temperature tendencies should be output.
549 : call phys_getopts( history_budget_out = history_budget, &
550 : history_budget_histfile_num_out = history_budget_histfile_num, &
551 : history_waccm_out = history_waccm, &
552 1536 : history_amwg_out = history_amwg )
553 :
554 : ! Initialize subordinate modules.
555 : call gw_common_init(pver,&
556 : tau_0_ubc, ktop, gravit, rair, alpha, gw_prndl, &
557 1536 : gw_qbo_hdepth_scaling, errstring )
558 : call shr_assert(trim(errstring) == "", "gw_common_init: "//errstring// &
559 1536 : errMsg(__FILE__, __LINE__))
560 :
561 1536 : if ( use_gw_oro ) then
562 :
563 1536 : if (effgw_oro == unset_r8) then
564 : call endrun("gw_init: Orographic gravity waves enabled, &
565 0 : &but effgw_oro was not set.")
566 : end if
567 : end if
568 :
569 1536 : if (use_gw_oro .or. use_gw_rdg_beta .or. use_gw_rdg_gamma) then
570 :
571 1536 : sgh_idx = pbuf_get_index('SGH')
572 :
573 : ! Declare history variables for orographic term
574 : call addfld ('TAUAORO', (/ 'ilev' /), 'I','N m-2', &
575 3072 : 'Total stress from original OGW scheme')
576 : call addfld ('TTGWORO', (/ 'lev' /), 'A','K s-1', &
577 3072 : 'T tendency - orographic gravity wave drag')
578 : call addfld ('TTGWSDFORO', (/ 'lev' /), 'A','K s-1', &
579 3072 : 'T tendency - orographic gravity wave, diffusion.')
580 : call addfld ('TTGWSKEORO', (/ 'lev' /), 'A','K s-1', &
581 3072 : 'T tendency - orographic gravity wave, breaking KE.')
582 : call addfld ('UTGWORO', (/ 'lev' /), 'A','m s-2', &
583 3072 : 'U tendency - orographic gravity wave drag')
584 : call addfld ('VTGWORO', (/ 'lev' /), 'A','m s-2', &
585 3072 : 'V tendency - orographic gravity wave drag')
586 1536 : call register_vector_field('UTGWORO', 'VTGWORO')
587 : call addfld ('TAUGWX', horiz_only, 'A','N m-2', &
588 1536 : 'Zonal gravity wave surface stress')
589 : call addfld ('TAUGWY', horiz_only, 'A','N m-2', &
590 1536 : 'Meridional gravity wave surface stress')
591 1536 : call register_vector_field('TAUGWX', 'TAUGWY')
592 :
593 1536 : if (history_amwg) then
594 1536 : call add_default('TAUGWX ', 1, ' ')
595 1536 : call add_default('TAUGWY ', 1, ' ')
596 : end if
597 :
598 1536 : if (history_waccm) then
599 0 : call add_default('UTGWORO ', 1, ' ')
600 0 : call add_default('VTGWORO ', 1, ' ')
601 0 : call add_default('TAUGWX ', 1, ' ')
602 0 : call add_default('TAUGWY ', 1, ' ')
603 : end if
604 :
605 : end if
606 :
607 1536 : if (use_gw_rdg_beta .or. use_gw_rdg_gamma) then
608 0 : grid_id = cam_grid_id('physgrid')
609 0 : if (.not. cam_grid_check(grid_id)) then
610 0 : call endrun(sub//': ERROR: no "physgrid" grid')
611 : end if
612 0 : call cam_grid_get_dim_names(grid_id, dim1name, dim2name)
613 : end if
614 :
615 1536 : if (use_gw_rdg_beta) then
616 :
617 0 : if (effgw_rdg_beta == unset_r8) then
618 : call endrun(sub//": ERROR: Anisotropic OGW enabled, &
619 0 : &but effgw_rdg_beta was not set.")
620 : end if
621 :
622 0 : fh_topo => topo_file_get_id()
623 0 : bnd_topo_loc = ' '
624 0 : if (.not. associated(fh_topo)) then
625 :
626 : ! Try to open topo file here. This workaround will not be needed
627 : ! once the refactored initialization sequence is on trunk.
628 :
629 0 : allocate(fh_topo)
630 : ! Error exit is from getfil if file not found.
631 0 : call getfil(bnd_topo, bnd_topo_loc)
632 0 : call cam_pio_openfile(fh_topo, bnd_topo_loc, PIO_NOWRITE)
633 :
634 : end if
635 :
636 : ! Get beta ridge data
637 : allocate( &
638 0 : rdg_gbxar(pcols,begchunk:endchunk), &
639 0 : rdg_hwdth(pcols,prdg,begchunk:endchunk), &
640 0 : rdg_clngt(pcols,prdg,begchunk:endchunk), &
641 0 : rdg_mxdis(pcols,prdg,begchunk:endchunk), &
642 0 : rdg_anixy(pcols,prdg,begchunk:endchunk), &
643 0 : rdg_angll(pcols,prdg,begchunk:endchunk) )
644 :
645 : call infld('GBXAR', fh_topo, dim1name, dim2name, 1, pcols, &
646 0 : begchunk, endchunk, rdg_gbxar, found, gridname='physgrid')
647 0 : if (.not. found) call endrun(sub//': ERROR: GBXAR not found on topo file')
648 0 : rdg_gbxar = rdg_gbxar * (rearth/1000._r8)*(rearth/1000._r8) ! transform to km^2
649 :
650 : call infld('HWDTH', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
651 0 : 1, prdg, begchunk, endchunk, rdg_hwdth, found, gridname='physgrid')
652 0 : if (.not. found) call endrun(sub//': ERROR: HWDTH not found on topo file')
653 :
654 : call infld('CLNGT', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
655 0 : 1, prdg, begchunk, endchunk, rdg_clngt, found, gridname='physgrid')
656 0 : if (.not. found) call endrun(sub//': ERROR: CLNGT not found on topo file')
657 :
658 : call infld('MXDIS', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
659 0 : 1, prdg, begchunk, endchunk, rdg_mxdis, found, gridname='physgrid')
660 0 : if (.not. found) call endrun(sub//': ERROR: MXDIS not found on topo file')
661 :
662 : call infld('ANIXY', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
663 0 : 1, prdg, begchunk, endchunk, rdg_anixy, found, gridname='physgrid')
664 0 : if (.not. found) call endrun(sub//': ERROR: ANIXY not found on topo file')
665 :
666 : call infld('ANGLL', fh_topo, dim1name, 'nrdg', dim2name, 1, pcols, &
667 0 : 1, prdg, begchunk, endchunk, rdg_angll, found, gridname='physgrid')
668 0 : if (.not. found) call endrun(sub//': ERROR: ANGLL not found on topo file')
669 :
670 : ! close topo file only if it was opened here
671 0 : if (len_trim(bnd_topo_loc) > 0) then
672 0 : call pio_closefile(fh_topo)
673 : end if
674 :
675 : call addfld ('WBR_HT1', horiz_only, 'I','m', &
676 0 : 'Wave breaking height for DSW')
677 : call addfld ('TLB_HT1', horiz_only, 'I','m', &
678 0 : 'Form drag layer height')
679 : call addfld ('BWV_HT1', horiz_only, 'I','m', &
680 0 : 'Bottom of freely-propagating OGW regime')
681 : call addfld ('TAUDSW1', horiz_only, 'I','Nm-2', &
682 0 : 'DSW enhanced drag')
683 : call addfld ('TAUORO1', horiz_only, 'I','Nm-2', &
684 0 : 'lower BC on propagating wave stress')
685 : call addfld ('UBMSRC1', horiz_only, 'I','ms-1', &
686 0 : 'below-peak-level on-ridge wind')
687 : call addfld ('USRC1', horiz_only, 'I','ms-1', &
688 0 : 'below-peak-level Zonal wind')
689 : call addfld ('VSRC1', horiz_only, 'I','ms-1', &
690 0 : 'below-peak-level Meridional wind')
691 : call addfld ('NSRC1', horiz_only, 'I','s-1', &
692 0 : 'below-peak-level stratification')
693 : call addfld ('MXDIS1', horiz_only, 'I','m', &
694 0 : 'Ridge/obstacle height')
695 : call addfld ('ANGLL1', horiz_only, 'I','degrees', &
696 0 : 'orientation clockwise w/resp north-south')
697 : call addfld ('ANIXY1', horiz_only, 'I','1', &
698 0 : 'Ridge quality')
699 : call addfld ('HWDTH1', horiz_only, 'I','km', &
700 0 : 'Ridge width')
701 : call addfld ('CLNGT1', horiz_only, 'I','km', &
702 0 : 'Ridge length')
703 : call addfld ('GBXAR1', horiz_only, 'I','km+2', &
704 0 : 'grid box area')
705 :
706 : call addfld ('Fr1_DIAG', horiz_only, 'I','1', &
707 0 : 'Critical Froude number for linear waves')
708 : call addfld ('Fr2_DIAG', horiz_only, 'I','1', &
709 0 : 'Critical Froude number for blocked flow')
710 : call addfld ('Frx_DIAG', horiz_only, 'I','1', &
711 0 : 'Obstacle Froude Number')
712 :
713 : call addfld('UEGW', (/ 'lev' /) , 'A' ,'s-1' , &
714 0 : 'Zonal wind profile-entry to GW ' )
715 : call addfld('VEGW', (/ 'lev' /) , 'A' ,'s-1' , &
716 0 : 'Merdional wind profile-entry to GW ' )
717 0 : call register_vector_field('UEGW','VEGW')
718 : call addfld('TEGW', (/ 'lev' /) , 'A' ,'K' , &
719 0 : 'Temperature profile-entry to GW ' )
720 : call addfld('ZEGW', (/ 'ilev' /) , 'A' ,'m' , &
721 0 : 'interface geopotential heights in GW code ' )
722 : call addfld('ZMGW', (/ 'lev' /) , 'A' ,'m' , &
723 0 : 'midlayer geopotential heights in GW code ' )
724 :
725 : call addfld('TAUM1_DIAG' , (/ 'ilev' /) , 'I' ,'N m-2' , &
726 0 : 'Ridge based momentum flux profile')
727 : call addfld('TAU1RDGBETAM' , (/ 'ilev' /) , 'I' ,'N m-2' , &
728 0 : 'Ridge based momentum flux profile')
729 : call addfld('UBM1BETA', (/ 'lev' /) , 'A' ,'s-1' , &
730 0 : 'On-ridge wind profile ' )
731 : call addfld('UBT1RDGBETA' , (/ 'lev' /) , 'I' ,'m s-1' , &
732 0 : 'On-ridge wind tendency from ridge 1 ')
733 :
734 0 : do i = 1, 6
735 0 : write(cn, '(i1)') i
736 : call addfld('TAU'//cn//'RDGBETAY' , (/ 'ilev' /), 'I', 'N m-2', &
737 0 : 'Ridge based momentum flux profile')
738 : call addfld('TAU'//cn//'RDGBETAX' , (/ 'ilev' /), 'I', 'N m-2', &
739 0 : 'Ridge based momentum flux profile')
740 0 : call register_vector_field('TAU'//cn//'RDGBETAX','TAU'//cn//'RDGBETAY')
741 : call addfld('UT'//cn//'RDGBETA', (/ 'lev' /), 'I', 'm s-1', &
742 0 : 'U wind tendency from ridge '//cn)
743 : call addfld('VT'//cn//'RDGBETA', (/ 'lev' /), 'I', 'm s-1', &
744 0 : 'V wind tendency from ridge '//cn)
745 0 : call register_vector_field('UT'//cn//'RDGBETA','VT'//cn//'RDGBETA')
746 : end do
747 :
748 : call addfld('TAUARDGBETAY' , (/ 'ilev' /) , 'I' ,'N m-2' , &
749 0 : 'Ridge based momentum flux profile')
750 : call addfld('TAUARDGBETAX' , (/ 'ilev' /) , 'I' ,'N m-2' , &
751 0 : 'Ridge based momentum flux profile')
752 0 : call register_vector_field('TAUARDGBETAX','TAUARDGBETAY')
753 :
754 0 : if (history_waccm) then
755 0 : call add_default('TAUARDGBETAX', 1, ' ')
756 0 : call add_default('TAUARDGBETAY ', 1, ' ')
757 : end if
758 :
759 : end if
760 :
761 1536 : if (use_gw_rdg_gamma) then
762 :
763 0 : if (effgw_rdg_gamma == unset_r8) then
764 0 : call endrun(sub//": ERROR: Anisotropic OGW enabled, but effgw_rdg_gamma was not set.")
765 : end if
766 :
767 0 : call getfil(bnd_rdggm, bnd_rdggm_loc, iflag=1, lexist=found)
768 0 : if (found) then
769 0 : call cam_pio_openfile(fh_rdggm, bnd_rdggm_loc, PIO_NOWRITE)
770 : else
771 : call endrun(sub//': ERROR: file for gamma ridges not found: bnd_rdggm='// &
772 0 : trim(bnd_rdggm))
773 : end if
774 :
775 0 : if (.not. allocated(rdg_gbxar)) then
776 0 : allocate(rdg_gbxar(pcols,begchunk:endchunk))
777 : call infld('GBXAR', fh_rdggm, dim1name, dim2name, 1, pcols, &
778 0 : begchunk, endchunk, rdg_gbxar, found, gridname='physgrid')
779 0 : if (.not. found) call endrun(sub//': ERROR: GBXAR not found on bnd_rdggm')
780 0 : rdg_gbxar = rdg_gbxar * (rearth/1000._r8)*(rearth/1000._r8) ! transform to km^2
781 : end if
782 :
783 : ! Get meso-gamma ridge data
784 : allocate( &
785 0 : rdg_hwdthg(pcols,prdg,begchunk:endchunk), &
786 0 : rdg_clngtg(pcols,prdg,begchunk:endchunk), &
787 0 : rdg_mxdisg(pcols,prdg,begchunk:endchunk), &
788 0 : rdg_anixyg(pcols,prdg,begchunk:endchunk), &
789 0 : rdg_angllg(pcols,prdg,begchunk:endchunk) )
790 :
791 : call infld('HWDTH', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
792 0 : 1, prdg, begchunk, endchunk, rdg_hwdthg, found, gridname='physgrid')
793 0 : if (.not. found) call endrun(sub//': ERROR: HWDTH not found on bnd_rdggm')
794 :
795 : call infld('CLNGT', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
796 0 : 1, prdg, begchunk, endchunk, rdg_clngtg, found, gridname='physgrid')
797 0 : if (.not. found) call endrun(sub//': ERROR: CLNGT not found on bnd_rdggm')
798 :
799 : call infld('MXDIS', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
800 0 : 1, prdg, begchunk, endchunk, rdg_mxdisg, found, gridname='physgrid')
801 0 : if (.not. found) call endrun(sub//': ERROR: MXDIS not found on bnd_rdggm')
802 :
803 : call infld('ANIXY', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
804 0 : 1, prdg, begchunk, endchunk, rdg_anixyg, found, gridname='physgrid')
805 0 : if (.not. found) call endrun(sub//': ERROR: ANIXY not found on bnd_rdggm')
806 :
807 : call infld('ANGLL', fh_rdggm, dim1name, 'nrdg', dim2name, 1, pcols, &
808 0 : 1, prdg, begchunk, endchunk, rdg_angllg, found, gridname='physgrid')
809 0 : if (.not. found) call endrun(sub//': ERROR: ANGLL not found on bnd_rdggm')
810 :
811 0 : call pio_closefile(fh_rdggm)
812 :
813 : call addfld ('TAU1RDGGAMMAM' , (/ 'ilev' /) , 'I' ,'N m-2' , &
814 0 : 'Ridge based momentum flux profile')
815 : call addfld ('UBM1GAMMA', (/ 'lev' /) , 'A' ,'s-1' , &
816 0 : 'On-ridge wind profile ' )
817 : call addfld ('UBT1RDGGAMMA' , (/ 'lev' /) , 'I' ,'m s-1' , &
818 0 : 'On-ridge wind tendency from ridge 1 ')
819 :
820 0 : do i = 1, 6
821 0 : write(cn, '(i1)') i
822 : call addfld('TAU'//cn//'RDGGAMMAY', (/ 'ilev' /), 'I', 'N m-2', &
823 0 : 'Ridge based momentum flux profile')
824 : call addfld('TAU'//cn//'RDGGAMMAX', (/ 'ilev' /), 'I', 'N m-2', &
825 0 : 'Ridge based momentum flux profile')
826 : call addfld('UT'//cn//'RDGGAMMA' , (/ 'lev' /), 'I', 'm s-1', &
827 0 : 'U wind tendency from ridge '//cn)
828 : call addfld('VT'//cn//'RDGGAMMA' , (/ 'lev' /), 'I', 'm s-1', &
829 0 : 'V wind tendency from ridge '//cn)
830 0 : call register_vector_field('UT'//cn//'RDGGAMMA','VT'//cn//'RDGGAMMA')
831 : end do
832 :
833 : call addfld ('TAUARDGGAMMAY' , (/ 'ilev' /) , 'I' ,'N m-2' , &
834 0 : 'Ridge based momentum flux profile')
835 : call addfld ('TAUARDGGAMMAX' , (/ 'ilev' /) , 'I' ,'N m-2' , &
836 0 : 'Ridge based momentum flux profile')
837 0 : call register_vector_field('TAUARDGGAMMAX','TAUARDGGAMMAY')
838 : call addfld ('TAURDGGMX', horiz_only, 'A','N m-2', &
839 0 : 'Zonal gravity wave surface stress')
840 : call addfld ('TAURDGGMY', horiz_only, 'A','N m-2', &
841 0 : 'Meridional gravity wave surface stress')
842 0 : call register_vector_field('TAURDGGMX','TAURDGGMY')
843 : call addfld ('UTRDGGM' , (/ 'lev' /) , 'I' ,'m s-1' , &
844 0 : 'U wind tendency from ridge 6 ')
845 : call addfld ('VTRDGGM' , (/ 'lev' /) , 'I' ,'m s-1' , &
846 0 : 'V wind tendency from ridge 6 ')
847 0 : call register_vector_field('UTRDGGM','VTRDGGM')
848 : end if
849 :
850 1536 : if (use_gw_front .or. use_gw_front_igw) then
851 :
852 0 : frontgf_idx = pbuf_get_index('FRONTGF')
853 0 : frontga_idx = pbuf_get_index('FRONTGA')
854 :
855 : call shr_assert(unset_r8 /= frontgfc, &
856 : "gw_init: Frontogenesis enabled, but frontgfc was &
857 : & not set!"// &
858 0 : errMsg(__FILE__, __LINE__))
859 :
860 0 : do k = 0, pver
861 : ! Check frontogenesis at 600 hPa.
862 0 : if (pref_edge(k+1) < 60000._r8) kfront = k+1
863 : end do
864 :
865 : ! Source waves from 500 hPa.
866 0 : kbot_front = maxloc(pref_edge, 1, (pref_edge < 50000._r8)) - 1
867 :
868 0 : if (masterproc) then
869 0 : write (iulog,*) 'KFRONT =',kfront
870 0 : write (iulog,*) 'KBOT_FRONT =',kbot_front
871 0 : write(iulog,*) ' '
872 : end if
873 :
874 : call addfld ('FRONTGF', (/ 'lev' /), 'A', 'K^2/M^2/S', &
875 0 : 'Frontogenesis function at gws src level')
876 : call addfld ('FRONTGFA', (/ 'lev' /), 'A', 'K^2/M^2/S', &
877 0 : 'Frontogenesis function at gws src level')
878 :
879 0 : if (history_waccm) then
880 0 : call add_default('FRONTGF', 1, ' ')
881 0 : call add_default('FRONTGFA', 1, ' ')
882 : end if
883 :
884 : end if
885 :
886 1536 : if (use_gw_front) then
887 :
888 : call shr_assert(all(unset_r8 /= [ effgw_cm, taubgnd ]), &
889 : "gw_init: Frontogenesis mid-scale waves enabled, but not &
890 : &all required namelist variables were set!"// &
891 0 : errMsg(__FILE__, __LINE__))
892 :
893 0 : if (masterproc) then
894 0 : write(iulog,*) 'gw_init: gw spectrum taubgnd, ', &
895 0 : 'effgw_cm = ',taubgnd, effgw_cm
896 0 : write(iulog,*) ' '
897 : end if
898 :
899 : cm_desc = gaussian_cm_desc(band_mid, kbot_front, kfront, frontgfc, &
900 0 : taubgnd, front_gaussian_width)
901 :
902 : ! Output for gravity waves from frontogenesis.
903 : call gw_spec_addflds(prefix=cm_pf, scheme="C&M", band=band_mid, &
904 0 : history_defaults=history_waccm)
905 :
906 : end if
907 :
908 1536 : if (use_gw_front_igw) then
909 :
910 : call shr_assert(all(unset_r8 /= [ effgw_cm_igw, taubgnd_igw ]), &
911 : "gw_init: Frontogenesis inertial waves enabled, but not &
912 : &all required namelist variables were set!"// &
913 0 : errMsg(__FILE__, __LINE__))
914 :
915 0 : if (masterproc) then
916 0 : write(iulog,*) 'gw_init: gw spectrum taubgnd_igw, ', &
917 0 : 'effgw_cm_igw = ',taubgnd_igw, effgw_cm_igw
918 0 : write(iulog,*) ' '
919 : end if
920 :
921 : cm_igw_desc = gaussian_cm_desc(band_long, kbot_front, kfront, frontgfc, &
922 0 : taubgnd_igw, front_gaussian_width)
923 :
924 : ! Output for gravity waves from frontogenesis.
925 : call gw_spec_addflds(prefix=cm_igw_pf, scheme="C&M IGW", &
926 0 : band=band_long, history_defaults=history_waccm)
927 :
928 : end if
929 :
930 : ! ========= Moving Mountain initialization! ==========================
931 1536 : if (use_gw_movmtn_pbl) then
932 :
933 : ! get pbuf indices for CLUBB couplings
934 0 : ttend_clubb_idx = pbuf_get_index('TTEND_CLUBB')
935 0 : thlp2_clubb_gw_idx = pbuf_get_index('THLP2_CLUBB_GW')
936 0 : upwp_clubb_gw_idx = pbuf_get_index('UPWP_CLUBB_GW')
937 0 : vpwp_clubb_gw_idx = pbuf_get_index('VPWP_CLUBB_GW')
938 0 : wpthlp_clubb_gw_idx = pbuf_get_index('WPTHLP_CLUBB_GW')
939 :
940 0 : if (masterproc) then
941 0 : write (iulog,*) 'Moving Mountain development code call init_movmtn'
942 : end if
943 :
944 :
945 : ! Confirm moving mountain file is enabled
946 : call shr_assert(trim(gw_drag_file_mm) /= "", &
947 : "gw_init: No gw_drag_file provided for DP GW moving mountain lookup &
948 : &table. Set this via namelist."// &
949 0 : errMsg(__FILE__, __LINE__))
950 :
951 0 : call gw_init_movmtn(gw_drag_file_mm, band_movmtn, movmtn_desc)
952 :
953 0 : do k = 0, pver
954 : ! 950 hPa index
955 0 : if (pref_edge(k+1) < 95000._r8) movmtn_desc%k = k+1
956 : end do
957 :
958 : ! Don't use deep convection heating depths below this limit.
959 0 : movmtn_desc%min_hdepth = 1._r8
960 0 : if (masterproc) then
961 0 : write (iulog,*) 'Moving mountain deep level =',movmtn_desc%k
962 : end if
963 :
964 : call addfld ('GWUT_MOVMTN',(/ 'lev' /), 'I','m s-2', &
965 0 : 'Mov Mtn dragforce - ubm component')
966 : call addfld ('UTGW_MOVMTN',(/ 'lev' /), 'I','m s-2', &
967 0 : 'Mov Mtn dragforce - u component')
968 : call addfld ('VTGW_MOVMTN',(/ 'lev' /), 'I','m s-2', &
969 0 : 'Mov Mtn dragforce - v component')
970 : call addfld('TAU_MOVMTN', (/ 'ilev' /), 'I', 'N m-2', &
971 0 : 'Moving Mountain momentum flux profile')
972 : call addfld('U_MOVMTN_IN', (/ 'lev' /), 'I', 'm s-1', &
973 0 : 'Moving Mountain - midpoint zonal input wind')
974 : call addfld('V_MOVMTN_IN', (/ 'lev' /), 'I', 'm s-1', &
975 0 : 'Moving Mountain - midpoint meridional input wind')
976 : call addfld('UBI_MOVMTN', (/ 'ilev' /), 'I', 'm s-1', &
977 0 : 'Moving Mountain - interface wind in direction of wave')
978 : call addfld('UBM_MOVMTN', (/ 'lev' /), 'I', 'm s-1', &
979 0 : 'Moving Mountain - midpoint wind in direction of wave')
980 : call addfld ('HDEPTH_MOVMTN',horiz_only,'I','km', &
981 0 : 'Heating Depth')
982 : call addfld ('UCELL_MOVMTN',horiz_only,'I','m s-1', &
983 0 : 'Gravity Wave Moving Mountain - Source-level X-wind')
984 : call addfld ('VCELL_MOVMTN',horiz_only,'I','m s-1', &
985 0 : 'Gravity Wave Moving Mountain - Source-level Y-wind')
986 : call addfld ('CS_MOVMTN',horiz_only,'I','m s-1', &
987 0 : 'Gravity Wave Moving Mountain - phase speed in direction of wave')
988 : call addfld ('STEER_LEVEL_MOVMTN',horiz_only,'I','1', &
989 0 : 'Gravity Wave Moving Mountain - steering level for movmtn GW')
990 : call addfld ('SRC_LEVEL_MOVMTN',horiz_only,'I','1', &
991 0 : 'Gravity Wave Moving Mountain - launch level for movmtn GW')
992 : call addfld ('TND_LEVEL_MOVMTN',horiz_only,'I','1', &
993 0 : 'Gravity Wave Moving Mountain - tendency lowest level for movmtn GW')
994 : call addfld ('NETDT_MOVMTN',(/ 'lev' /),'I','K s-1', &
995 0 : 'Gravity Wave Moving Mountain - Net heating rate')
996 : call addfld ('TTEND_CLUBB',(/ 'lev' /),'A','K s-1', &
997 0 : 'Gravity Wave Moving Mountain - CLUBB Net heating rate')
998 : call addfld ('THLP2_CLUBB_GW',(/ 'ilev' /),'A','K+2', &
999 0 : 'Gravity Wave Moving Mountain - THLP variance from CLUBB to GW')
1000 : call addfld ('WPTHLP_CLUBB_GW',(/ 'ilev' /),'A','Km s-2', &
1001 0 : 'Gravity Wave Moving Mountain - WPTHLP from CLUBB to GW')
1002 : call addfld ('UPWP_CLUBB_GW',(/ 'ilev' /),'A','m+2 s-2', &
1003 0 : 'Gravity Wave Moving Mountain - X-momflux from CLUBB to GW')
1004 : call addfld ('VPWP_CLUBB_GW',(/ 'ilev' /),'A','m+2 s-2', &
1005 0 : 'Gravity Wave Moving Mountain - Y-momflux from CLUBB to GW')
1006 : call addfld ('XPWP_SRC_MOVMTN',horiz_only,'I','m+2 s-2', &
1007 0 : 'Gravity Wave Moving Mountain - flux source for moving mtn')
1008 :
1009 : end if
1010 :
1011 1536 : if (use_gw_convect_dp) then
1012 :
1013 0 : ttend_dp_idx = pbuf_get_index('TTEND_DP')
1014 :
1015 : ! Set the deep scheme specification components.
1016 0 : beres_dp_desc%storm_shift = .true.
1017 :
1018 0 : do k = 0, pver
1019 : ! 700 hPa index
1020 0 : if (pref_edge(k+1) < 70000._r8) beres_dp_desc%k = k+1
1021 : end do
1022 :
1023 0 : if (masterproc) then
1024 0 : write (iulog,*) 'Beres deep level =',beres_dp_desc%k
1025 : end if
1026 :
1027 : ! Don't use deep convection heating depths below this limit.
1028 : ! This is important for QBO. Bad result if left at 2.5 km.
1029 0 : beres_dp_desc%min_hdepth = 1000._r8
1030 :
1031 : ! Read Beres file.
1032 :
1033 : call shr_assert(trim(gw_drag_file) /= "", &
1034 : "gw_init: No gw_drag_file provided for Beres deep &
1035 : &scheme. Set this via namelist."// &
1036 0 : errMsg(__FILE__, __LINE__))
1037 :
1038 0 : call gw_init_beres(gw_drag_file, band_mid, beres_dp_desc)
1039 :
1040 : ! Output for gravity waves from the Beres scheme (deep).
1041 : call gw_spec_addflds(prefix=beres_dp_pf, scheme="Beres (deep)", &
1042 0 : band=band_mid, history_defaults=history_waccm)
1043 :
1044 : call addfld ('NETDT',(/ 'lev' /), 'A','K s-1', &
1045 0 : 'Net heating rate')
1046 : call addfld ('MAXQ0',horiz_only , 'A','K day-1', &
1047 0 : 'Max column heating rate')
1048 : call addfld ('HDEPTH',horiz_only, 'A','km', &
1049 0 : 'Heating Depth')
1050 :
1051 0 : if (history_waccm) then
1052 0 : call add_default('NETDT ', 1, ' ')
1053 0 : call add_default('HDEPTH ', 1, ' ')
1054 0 : call add_default('MAXQ0 ', 1, ' ')
1055 : end if
1056 :
1057 : end if
1058 :
1059 1536 : if (use_gw_convect_sh) then
1060 :
1061 0 : ttend_sh_idx = pbuf_get_index('TTEND_SH')
1062 :
1063 : ! Set the shallow scheme specification components.
1064 0 : beres_sh_desc%storm_shift = .false.
1065 :
1066 0 : do k = 0, pver
1067 : ! 900 hPa index
1068 0 : if (pref_edge(k+1) < 90000._r8) beres_sh_desc%k = k+1
1069 : end do
1070 :
1071 0 : if (masterproc) then
1072 0 : write (iulog,*) 'Beres shallow level =',beres_sh_desc%k
1073 : end if
1074 :
1075 : ! Use all heating depths for shallow convection.
1076 0 : beres_sh_desc%min_hdepth = 0._r8
1077 :
1078 : ! Read Beres file.
1079 :
1080 : call shr_assert(trim(gw_drag_file_sh) /= "", &
1081 : "gw_init: No gw_drag_file_sh provided for Beres shallow &
1082 : &scheme. Set this via namelist."// &
1083 0 : errMsg(__FILE__, __LINE__))
1084 :
1085 0 : call gw_init_beres(gw_drag_file_sh, band_mid, beres_sh_desc)
1086 :
1087 : ! Output for gravity waves from the Beres scheme (shallow).
1088 : call gw_spec_addflds(prefix=beres_sh_pf, scheme="Beres (shallow)", &
1089 0 : band=band_mid, history_defaults=history_waccm)
1090 :
1091 : call addfld ('SNETDT',(/ 'lev' /), 'A','K s-1', &
1092 0 : 'Net heating rate')
1093 : call addfld ('SMAXQ0',horiz_only , 'A','K day-1', &
1094 0 : 'Max column heating rate')
1095 : call addfld ('SHDEPTH',horiz_only, 'A','km', &
1096 0 : 'Heating Depth')
1097 :
1098 0 : if (history_waccm) then
1099 0 : call add_default('SNETDT ', 1, ' ')
1100 0 : call add_default('SHDEPTH ', 1, ' ')
1101 0 : call add_default('SMAXQ0 ', 1, ' ')
1102 : end if
1103 :
1104 : end if
1105 :
1106 : call addfld ('EKGW' ,(/ 'ilev' /), 'A','M2/S', &
1107 3072 : 'Effective Kzz due to diffusion by gravity waves')
1108 :
1109 1536 : if (history_waccm) then
1110 0 : call add_default('EKGW', 1, ' ')
1111 : end if
1112 :
1113 : call addfld ('UTGW_TOTAL', (/ 'lev' /), 'A','m s-2', &
1114 3072 : 'Total U tendency due to gravity wave drag')
1115 : call addfld ('VTGW_TOTAL', (/ 'lev' /), 'A','m s-2', &
1116 3072 : 'Total V tendency due to gravity wave drag')
1117 1536 : call register_vector_field('UTGW_TOTAL', 'VTGW_TOTAL')
1118 :
1119 : ! Total temperature tendency output.
1120 : call addfld ('TTGW', (/ 'lev' /), 'A', 'K s-1', &
1121 3072 : 'T tendency - gravity wave drag')
1122 :
1123 : ! Water budget terms.
1124 : call addfld('QTGW',(/ 'lev' /), 'A','kg/kg/s', &
1125 3072 : 'Q tendency - gravity wave drag')
1126 : call addfld('CLDLIQTGW',(/ 'lev' /), 'A','kg/kg/s', &
1127 3072 : 'CLDLIQ tendency - gravity wave drag')
1128 : call addfld('CLDICETGW',(/ 'lev' /), 'A','kg/kg/s', &
1129 3072 : 'CLDICE tendency - gravity wave drag')
1130 :
1131 1536 : if ( history_budget ) then
1132 0 : call add_default('TTGW', history_budget_histfile_num, ' ')
1133 0 : call add_default('QTGW', history_budget_histfile_num, ' ')
1134 0 : call add_default('CLDLIQTGW', history_budget_histfile_num, ' ')
1135 0 : call add_default('CLDICETGW', history_budget_histfile_num, ' ')
1136 : end if
1137 :
1138 : ! Get indices to actually output the above.
1139 1536 : call cnst_get_ind("CLDLIQ", ixcldliq)
1140 1536 : call cnst_get_ind("CLDICE", ixcldice)
1141 :
1142 1536 : if (gw_top_taper) then
1143 0 : allocate(vramp(pver))
1144 0 : vramp(:) = 1._r8
1145 0 : topndx = 1
1146 0 : botndx = press_lim_idx( 0.6E-02_r8, top=.true. )
1147 0 : if (botndx>1) then
1148 0 : do k=botndx,topndx,-1
1149 0 : vramp(k) = vramp(k+1)/(pref_edge(k+1)/pref_edge(k))
1150 : end do
1151 0 : if (masterproc) then
1152 0 : write(iulog,'(A)') 'GW taper coef (vramp):'
1153 0 : do k=1,pver
1154 0 : write(iulog,"('k: ',I4,' taper coef,press(Pa): ',F12.8,E12.4)") k, vramp(k), pref_mid(k)
1155 : enddo
1156 : endif
1157 : endif
1158 : end if
1159 :
1160 3072 : end subroutine gw_init
1161 :
1162 : !==========================================================================
1163 :
1164 0 : subroutine gw_init_beres(file_name, band, desc)
1165 :
1166 1536 : use ioFileMod, only: getfil
1167 : use pio, only: file_desc_t, pio_nowrite, pio_inq_varid, pio_get_var, &
1168 : pio_closefile
1169 : use cam_pio_utils, only: cam_pio_openfile
1170 :
1171 : character(len=*), intent(in) :: file_name
1172 : type(GWBand), intent(in) :: band
1173 :
1174 : type(BeresSourceDesc), intent(inout) :: desc
1175 :
1176 : type(file_desc_t) :: gw_file_desc
1177 :
1178 : ! PIO variable ids and error code.
1179 : integer :: mfccid, hdid, stat
1180 :
1181 : ! Number of wavenumbers in the input file.
1182 : integer :: ngwv_file
1183 :
1184 : ! Full path to gw_drag_file.
1185 : character(len=cl) :: file_path
1186 :
1187 : character(len=cl) :: msg
1188 :
1189 : !----------------------------------------------------------------------
1190 : ! read in look-up table for source spectra
1191 : !-----------------------------------------------------------------------
1192 :
1193 0 : call getfil(file_name, file_path)
1194 0 : call cam_pio_openfile(gw_file_desc, file_path, pio_nowrite)
1195 :
1196 : ! Get HD (heating depth) dimension.
1197 :
1198 0 : desc%maxh = get_pio_dimlen(gw_file_desc, "HD", file_path)
1199 :
1200 : ! Get MW (mean wind) dimension.
1201 :
1202 0 : desc%maxuh = get_pio_dimlen(gw_file_desc, "MW", file_path)
1203 :
1204 : ! Get PS (phase speed) dimension.
1205 :
1206 0 : ngwv_file = get_pio_dimlen(gw_file_desc, "PS", file_path)
1207 :
1208 : ! Number in each direction is half of total (and minus phase speed of 0).
1209 0 : desc%maxuh = (desc%maxuh-1)/2
1210 0 : ngwv_file = (ngwv_file-1)/2
1211 :
1212 : call shr_assert(ngwv_file >= band%ngwv, &
1213 : "gw_init_beres: PhaseSpeed in lookup table file does not cover the whole &
1214 0 : &spectrum implied by the model's ngwv. ")
1215 :
1216 : ! Allocate hd and get data.
1217 :
1218 0 : allocate(desc%hd(desc%maxh), stat=stat, errmsg=msg)
1219 :
1220 : call shr_assert(stat == 0, &
1221 : "gw_init_beres: Allocation error (hd): "//msg// &
1222 0 : errMsg(__FILE__, __LINE__))
1223 :
1224 0 : stat = pio_inq_varid(gw_file_desc,'HD',hdid)
1225 :
1226 : call handle_pio_error(stat, &
1227 0 : 'Error finding HD in: '//trim(file_path))
1228 :
1229 : stat = pio_get_var(gw_file_desc, hdid, start=[1], count=[desc%maxh], &
1230 0 : ival=desc%hd)
1231 :
1232 : call handle_pio_error(stat, &
1233 0 : 'Error reading HD from: '//trim(file_path))
1234 :
1235 : ! While not currently documented in the file, it uses kilometers. Convert
1236 : ! to meters.
1237 0 : desc%hd = desc%hd*1000._r8
1238 :
1239 : ! Allocate mfcc. "desc%maxh" and "desc%maxuh" are from the file, but the
1240 : ! model determines wavenumber dimension.
1241 :
1242 : allocate(desc%mfcc(desc%maxh,-desc%maxuh:desc%maxuh,&
1243 0 : -band%ngwv:band%ngwv), stat=stat, errmsg=msg)
1244 :
1245 : call shr_assert(stat == 0, &
1246 : "gw_init_beres: Allocation error (mfcc): "//msg// &
1247 0 : errMsg(__FILE__, __LINE__))
1248 :
1249 : ! Get mfcc data.
1250 :
1251 0 : stat = pio_inq_varid(gw_file_desc,'mfcc',mfccid)
1252 :
1253 : call handle_pio_error(stat, &
1254 0 : 'Error finding mfcc in: '//trim(file_path))
1255 :
1256 : stat = pio_get_var(gw_file_desc, mfccid, &
1257 0 : start=[1,1,ngwv_file-band%ngwv+1], count=shape(desc%mfcc), &
1258 0 : ival=desc%mfcc)
1259 :
1260 : call handle_pio_error(stat, &
1261 0 : 'Error reading mfcc from: '//trim(file_path))
1262 :
1263 0 : call pio_closefile(gw_file_desc)
1264 :
1265 0 : if (masterproc) then
1266 :
1267 0 : write(iulog,*) "Read in source spectra from file."
1268 0 : write(iulog,*) "mfcc max, min = ", &
1269 0 : maxval(desc%mfcc),", ",minval(desc%mfcc)
1270 :
1271 : endif
1272 :
1273 0 : end subroutine gw_init_beres
1274 :
1275 : !==============================================================
1276 0 : subroutine gw_init_movmtn(file_name, band, desc)
1277 :
1278 0 : use ioFileMod, only: getfil
1279 : use pio, only: file_desc_t, pio_nowrite, pio_inq_varid, pio_get_var, &
1280 : pio_closefile
1281 : use cam_pio_utils, only: cam_pio_openfile
1282 :
1283 : character(len=*), intent(in) :: file_name
1284 : type(GWBand), intent(in) :: band
1285 :
1286 : type(MovMtnSourceDesc), intent(inout) :: desc
1287 :
1288 : type(file_desc_t) :: gw_file_desc
1289 :
1290 : ! PIO variable ids and error code.
1291 : integer :: mfccid, uhid, hdid, stat
1292 :
1293 : ! Number of wavenumbers in the input file.
1294 : integer :: ngwv_file
1295 :
1296 : ! Full path to gw_drag_file.
1297 : character(len=cl) :: file_path
1298 :
1299 : character(len=cl) :: msg
1300 :
1301 : !----------------------------------------------------------------------
1302 : ! read in look-up table for source spectra
1303 : !-----------------------------------------------------------------------
1304 :
1305 0 : call getfil(file_name, file_path)
1306 :
1307 0 : call cam_pio_openfile(gw_file_desc, file_path, pio_nowrite)
1308 :
1309 : ! Get HD (heating depth) dimension.
1310 :
1311 0 : desc%maxh = 15 !get_pio_dimlen(gw_file_desc, "HD", file_path)
1312 :
1313 : ! Get MW (mean wind) dimension.
1314 :
1315 0 : desc%maxuh = 241 ! get_pio_dimlen(gw_file_desc, "MW", file_path)
1316 :
1317 : ! Get PS (phase speed) dimension.
1318 :
1319 0 : ngwv_file = 0 !get_pio_dimlen(gw_file_desc, "PS", file_path)
1320 :
1321 : ! Number in each direction is half of total (and minus phase speed of 0).
1322 0 : desc%maxuh = (desc%maxuh-1)/2
1323 0 : ngwv_file = (ngwv_file-1)/2
1324 :
1325 : call shr_assert(ngwv_file >= band%ngwv, &
1326 0 : "gw_movmtn_init: PhaseSpeed in lookup table inconsistent with moving mountain")
1327 :
1328 : ! Allocate hd and get data.
1329 :
1330 0 : allocate(desc%hd(desc%maxh), stat=stat, errmsg=msg)
1331 :
1332 : call shr_assert(stat == 0, &
1333 : "gw_init_movmtn: Allocation error (hd): "//msg// &
1334 0 : errMsg(__FILE__, __LINE__))
1335 :
1336 0 : stat = pio_inq_varid(gw_file_desc,'HDEPTH',hdid)
1337 :
1338 : call handle_pio_error(stat, &
1339 0 : 'Error finding HD in: '//trim(file_path))
1340 :
1341 : stat = pio_get_var(gw_file_desc, hdid, start=[1], count=[desc%maxh], &
1342 0 : ival=desc%hd)
1343 :
1344 : call handle_pio_error(stat, &
1345 0 : 'Error reading HD from: '//trim(file_path))
1346 :
1347 : ! While not currently documented in the file, it uses kilometers. Convert
1348 : ! to meters.
1349 0 : desc%hd = desc%hd*1000._r8
1350 :
1351 : ! Allocate wind and get data.
1352 :
1353 0 : allocate(desc%uh(desc%maxuh), stat=stat, errmsg=msg)
1354 :
1355 : call shr_assert(stat == 0, &
1356 : "gw_init_movmtn: Allocation error (uh): "//msg// &
1357 0 : errMsg(__FILE__, __LINE__))
1358 :
1359 0 : stat = pio_inq_varid(gw_file_desc,'UARR',uhid)
1360 :
1361 : call handle_pio_error(stat, &
1362 0 : 'Error finding UH in: '//trim(file_path))
1363 :
1364 : stat = pio_get_var(gw_file_desc, uhid, start=[1], count=[desc%maxuh], &
1365 0 : ival=desc%uh)
1366 :
1367 : call handle_pio_error(stat, &
1368 0 : 'Error reading UH from: '//trim(file_path))
1369 :
1370 : ! Allocate mfcc. "desc%maxh" and "desc%maxuh" are from the file, but the
1371 : ! model determines wavenumber dimension.
1372 :
1373 : allocate(desc%mfcc(desc%maxh,-desc%maxuh:desc%maxuh,&
1374 0 : -band%ngwv:band%ngwv), stat=stat, errmsg=msg)
1375 :
1376 : call shr_assert(stat == 0, &
1377 : "gw_init_movmtn: Allocation error (mfcc): "//msg// &
1378 0 : errMsg(__FILE__, __LINE__))
1379 :
1380 : ! Get mfcc data.
1381 :
1382 0 : stat = pio_inq_varid(gw_file_desc,'NEWMF',mfccid)
1383 :
1384 : call handle_pio_error(stat, &
1385 0 : 'Error finding mfcc in: '//trim(file_path))
1386 :
1387 : stat = pio_get_var(gw_file_desc, mfccid, &
1388 0 : start=[1,1], count=shape(desc%mfcc), &
1389 0 : ival=desc%mfcc)
1390 :
1391 : call handle_pio_error(stat, &
1392 0 : 'Error reading mfcc from: '//trim(file_path))
1393 :
1394 0 : call pio_closefile(gw_file_desc)
1395 :
1396 0 : if (masterproc) then
1397 :
1398 0 : write(iulog,*) "Read in Mov Mountain source file."
1399 :
1400 : endif
1401 :
1402 0 : end subroutine gw_init_movmtn
1403 : !==========================================================================
1404 :
1405 : ! Utility to reduce the repetitiveness of reads during initialization.
1406 0 : function get_pio_dimlen(file_desc, dim_name, file_path) result(dimlen)
1407 :
1408 0 : use pio, only: file_desc_t, pio_inq_dimid, pio_inq_dimlen
1409 :
1410 : type(file_desc_t), intent(in) :: file_desc
1411 : character(len=*), intent(in) :: dim_name
1412 :
1413 : ! File path, for use in error messages only.
1414 : character(len=*), intent(in) :: file_path
1415 :
1416 : integer :: dimlen
1417 :
1418 : integer :: dimid, stat
1419 :
1420 0 : stat = pio_inq_dimid(file_desc, dim_name, dimid)
1421 :
1422 : call handle_pio_error(stat, &
1423 0 : "Error finding dimension "//dim_name//" in: "//file_path)
1424 :
1425 0 : stat = pio_inq_dimlen(file_desc, dimid, dimlen)
1426 :
1427 : call handle_pio_error(stat, &
1428 0 : "Error reading dimension "//dim_name//" from: "//file_path)
1429 :
1430 0 : end function get_pio_dimlen
1431 :
1432 : !==========================================================================
1433 :
1434 : ! In fact, we'd usually expect PIO errors to abort the run before you can
1435 : ! even check the error code. But just in case, use this little assert.
1436 0 : subroutine handle_pio_error(stat, message)
1437 : use pio, only: pio_noerr
1438 : integer, intent(in) :: stat
1439 : character(len=*) :: message
1440 :
1441 : call shr_assert(stat == pio_noerr, &
1442 : "PIO error in gw_init_beres: "//trim(message)// &
1443 0 : errMsg(__FILE__, __LINE__))
1444 :
1445 0 : end subroutine handle_pio_error
1446 :
1447 : !==========================================================================
1448 :
1449 5956704 : subroutine gw_tend(state, pbuf, dt, ptend, cam_in, flx_heat)
1450 : !-----------------------------------------------------------------------
1451 : ! Interface for multiple gravity wave drag parameterization.
1452 : !-----------------------------------------------------------------------
1453 :
1454 : use physics_types, only: physics_state_copy, set_dry_to_wet
1455 : use constituents, only: cnst_type
1456 : use physics_buffer, only: physics_buffer_desc, pbuf_get_field
1457 : use camsrfexch, only: cam_in_t
1458 : ! Location-dependent cpair
1459 : use air_composition, only: cpairv
1460 : use physconst, only: pi
1461 : use coords_1d, only: Coords1D
1462 : use gw_common, only: gw_prof, gw_drag_prof, calc_taucd
1463 : use gw_common, only: momentum_flux, momentum_fixer, energy_change
1464 : use gw_common, only: energy_fixer, coriolis_speed, adjust_inertial
1465 : use gw_oro, only: gw_oro_src
1466 : use gw_front, only: gw_cm_src
1467 : use gw_convect, only: gw_beres_src
1468 : use gw_movmtn, only: gw_movmtn_src
1469 :
1470 : !------------------------------Arguments--------------------------------
1471 : type(physics_state), intent(in) :: state ! physics state structure
1472 : type(physics_buffer_desc), pointer :: pbuf(:) ! Physics buffer
1473 : real(r8), intent(in) :: dt ! time step
1474 : ! Parameterization net tendencies.
1475 : type(physics_ptend), intent(out):: ptend
1476 : type(cam_in_t), intent(in) :: cam_in
1477 : real(r8), intent(out) :: flx_heat(pcols)
1478 :
1479 : !---------------------------Local storage-------------------------------
1480 :
1481 1489176 : type(physics_state) :: state1 ! Local copy of state variable
1482 :
1483 : integer :: lchnk ! chunk identifier
1484 : integer :: ncol ! number of atmospheric columns
1485 : integer :: istat
1486 :
1487 : integer :: i, k ! loop indices
1488 :
1489 1489176 : type(Coords1D) :: p ! Pressure coordinates
1490 :
1491 2978352 : real(r8) :: ttgw(state%ncol,pver) ! temperature tendency
1492 2978352 : real(r8) :: utgw(state%ncol,pver) ! zonal wind tendency
1493 2978352 : real(r8) :: vtgw(state%ncol,pver) ! meridional wind tendency
1494 :
1495 2978352 : real(r8) :: ni(state%ncol,pver+1) ! interface Brunt-Vaisalla frequency
1496 2978352 : real(r8) :: nm(state%ncol,pver) ! midpoint Brunt-Vaisalla frequency
1497 2978352 : real(r8) :: rhoi(state%ncol,pver+1) ! interface density
1498 1489176 : real(r8), allocatable :: tau(:,:,:) ! wave Reynolds stress
1499 2978352 : real(r8) :: tau0x(state%ncol) ! c=0 sfc. stress (zonal)
1500 2978352 : real(r8) :: tau0y(state%ncol) ! c=0 sfc. stress (meridional)
1501 2978352 : real(r8) :: ubi(state%ncol,pver+1)! projection of wind at interfaces
1502 2978352 : real(r8) :: ubm(state%ncol,pver) ! projection of wind at midpoints
1503 2978352 : real(r8) :: xv(state%ncol) ! unit vector of source wind (x)
1504 2978352 : real(r8) :: yv(state%ncol) ! unit vector of source wind (y)
1505 :
1506 : integer :: m ! dummy integers
1507 2978352 : real(r8) :: qtgw(state%ncol,pver,pcnst) ! constituents tendencies
1508 :
1509 : ! Reynolds stress for waves propagating in each cardinal direction.
1510 2978352 : real(r8) :: taucd(state%ncol,pver+1,4)
1511 :
1512 : ! gravity wave wind tendency for each wave
1513 1489176 : real(r8), allocatable :: gwut(:,:,:)
1514 :
1515 : ! Temperature tendencies from diffusion and kinetic energy.
1516 2978352 : real(r8) :: dttdf(state%ncol,pver)
1517 2978352 : real(r8) :: dttke(state%ncol,pver)
1518 :
1519 : ! Wave phase speeds for each column
1520 1489176 : real(r8), allocatable :: phase_speeds(:,:)
1521 :
1522 : ! Efficiency for a gravity wave source.
1523 2978352 : real(r8) :: effgw(state%ncol)
1524 :
1525 : ! Coriolis characteristic speed.
1526 2978352 : real(r8) :: u_coriolis(state%ncol)
1527 :
1528 : ! Adjustment for inertial gravity waves.
1529 1489176 : real(r8), allocatable :: ro_adjust(:,:,:)
1530 :
1531 : ! pbuf fields
1532 : ! Molecular diffusivity
1533 1489176 : real(r8), pointer :: kvt_in(:,:)
1534 2978352 : real(r8) :: kvtt(state%ncol,pver+1)
1535 :
1536 : ! Frontogenesis
1537 1489176 : real(r8), pointer :: frontgf(:,:)
1538 1489176 : real(r8), pointer :: frontga(:,:)
1539 :
1540 : ! Temperature change due to deep convection.
1541 1489176 : real(r8), pointer :: ttend_dp(:,:)
1542 : ! Temperature change due to shallow convection.
1543 1489176 : real(r8), pointer :: ttend_sh(:,:)
1544 :
1545 : ! New couplings from CLUBB
1546 1489176 : real(r8), pointer :: ttend_clubb(:,:)
1547 1489176 : real(r8), pointer :: thlp2_clubb_gw(:,:)
1548 1489176 : real(r8), pointer :: wpthlp_clubb_gw(:,:)
1549 1489176 : real(r8), pointer :: upwp_clubb_gw(:,:)
1550 1489176 : real(r8), pointer :: vpwp_clubb_gw(:,:)
1551 2978352 : real(r8) :: xpwp_clubb(state%ncol,pver+1)
1552 :
1553 :
1554 : ! Standard deviation of orography.
1555 1489176 : real(r8), pointer :: sgh(:)
1556 :
1557 : ! gridbox area
1558 1489176 : real(r8), pointer :: gbxar(:)
1559 :
1560 : ! Beta ridges
1561 : ! width of ridges.
1562 1489176 : real(r8), pointer :: hwdth(:,:)
1563 : ! length of ridges.
1564 1489176 : real(r8), pointer :: clngt(:,:)
1565 : ! Maximum deviations of ridges.
1566 1489176 : real(r8), pointer :: mxdis(:,:)
1567 : ! orientation of ridges.
1568 1489176 : real(r8), pointer :: angll(:,:)
1569 : ! anisotropy of ridges.
1570 1489176 : real(r8), pointer :: anixy(:,:)
1571 :
1572 : ! Gamma ridges
1573 : ! width of ridges.
1574 1489176 : real(r8), pointer :: hwdthg(:,:)
1575 : ! length of ridges.
1576 1489176 : real(r8), pointer :: clngtg(:,:)
1577 : ! Maximum deviations of ridges.
1578 1489176 : real(r8), pointer :: mxdisg(:,:)
1579 : ! orientation of ridges.
1580 1489176 : real(r8), pointer :: angllg(:,:)
1581 : ! anisotropy of ridges.
1582 1489176 : real(r8), pointer :: anixyg(:,:)
1583 :
1584 : ! Indices of gravity wave source and lowest level where wind tendencies
1585 : ! are allowed.
1586 2978352 : integer :: src_level(state%ncol)
1587 2978352 : integer :: tend_level(state%ncol)
1588 :
1589 : ! Convective source heating depth.
1590 : ! heating depth
1591 2978352 : real(r8) :: hdepth(state%ncol)
1592 : ! maximum heating rate
1593 2978352 : real(r8) :: maxq0(state%ncol)
1594 :
1595 : ! Scale sgh to account for landfrac.
1596 2978352 : real(r8) :: sgh_scaled(state%ncol)
1597 :
1598 : ! Parameters for the IGW polar taper.
1599 : real(r8), parameter :: degree2radian = pi/180._r8
1600 : real(r8), parameter :: al0 = 82.5_r8 * degree2radian
1601 : real(r8), parameter :: dlat0 = 5.0_r8 * degree2radian
1602 :
1603 : ! effective gw diffusivity at interfaces needed for output
1604 2978352 : real(r8) :: egwdffi(state%ncol,pver+1)
1605 : ! sum from the two types of spectral GW
1606 2978352 : real(r8) :: egwdffi_tot(state%ncol,pver+1)
1607 :
1608 : ! Momentum fluxes used by fixer.
1609 2978352 : real(r8) :: um_flux(state%ncol), vm_flux(state%ncol)
1610 : ! Energy change used by fixer.
1611 2978352 : real(r8) :: de(state%ncol)
1612 :
1613 : ! Which constituents are being affected by diffusion.
1614 : logical :: lq(pcnst)
1615 :
1616 : ! Contiguous copies of state arrays.
1617 2978352 : real(r8) :: dse(state%ncol,pver)
1618 2978352 : real(r8) :: t(state%ncol,pver)
1619 2978352 : real(r8) :: u(state%ncol,pver)
1620 2978352 : real(r8) :: v(state%ncol,pver)
1621 2978352 : real(r8) :: q(state%ncol,pver,pcnst)
1622 2978352 : real(r8) :: piln(state%ncol,pver+1)
1623 2978352 : real(r8) :: zm(state%ncol,pver)
1624 2978352 : real(r8) :: zi(state%ncol,pver+1)
1625 :
1626 : !------------------------------------------------------------------------
1627 :
1628 : ! Make local copy of input state.
1629 1489176 : call physics_state_copy(state, state1)
1630 :
1631 : ! constituents are all treated as wet mmr
1632 1489176 : call set_dry_to_wet(state1, convert_cnst_type='dry')
1633 :
1634 1489176 : lchnk = state1%lchnk
1635 1489176 : ncol = state1%ncol
1636 :
1637 1489176 : p = Coords1D(state1%pint(:ncol,:))
1638 :
1639 647999352 : dse = state1%s(:ncol,:)
1640 647999352 : t = state1%t(:ncol,:)
1641 647999352 : u = state1%u(:ncol,:)
1642 647999352 : v = state1%v(:ncol,:)
1643 1945487232 : q = state1%q(:ncol,:,:)
1644 672865128 : piln = state1%lnpint(:ncol,:)
1645 647999352 : zm = state1%zm(:ncol,:)
1646 672865128 : zi = state1%zi(:ncol,:)
1647 :
1648 5956704 : lq = .true.
1649 : call physics_ptend_init(ptend, state1%psetcols, "Gravity wave drag", &
1650 1489176 : ls=.true., lu=.true., lv=.true., lq=lq)
1651 :
1652 : ! Profiles of background state variables
1653 1489176 : call gw_prof(ncol, p, cpair, t, rhoi, nm, ni)
1654 :
1655 1489176 : if (do_molec_diff) then
1656 : !--------------------------------------------------------
1657 : ! Initialize and calculate local molecular diffusivity
1658 : !--------------------------------------------------------
1659 :
1660 0 : call pbuf_get_field(pbuf, kvt_idx, kvt_in) ! kvt_in(1:pcols,1:pver+1)
1661 :
1662 : ! Set kvtt from pbuf field; kvtt still needs a factor of 1/cpairv.
1663 0 : kvtt = kvt_in(:ncol,:)
1664 :
1665 : ! Use linear extrapolation of cpairv to top interface.
1666 : kvtt(:,1) = kvtt(:,1) / &
1667 0 : (1.5_r8*cpairv(:ncol,1,lchnk) - &
1668 0 : 0.5_r8*cpairv(:ncol,2,lchnk))
1669 :
1670 : ! Interpolate cpairv to other interfaces.
1671 0 : do k = 2, nbot_molec
1672 0 : kvtt(:,k) = kvtt(:,k) / &
1673 0 : (cpairv(:ncol,k+1,lchnk)+cpairv(:ncol,k,lchnk)) * 2._r8
1674 : enddo
1675 :
1676 : else
1677 :
1678 672865128 : kvtt = 0._r8
1679 :
1680 : end if
1681 :
1682 1489176 : if (use_gw_front_igw) then
1683 0 : u_coriolis = coriolis_speed(band_long, state1%lat(:ncol))
1684 : end if
1685 :
1686 : ! Totals that accumulate over different sources.
1687 672865128 : egwdffi_tot = 0._r8
1688 1489176 : flx_heat = 0._r8
1689 :
1690 1489176 : if (use_gw_movmtn_pbl) then
1691 : !------------------------------------------------------------------
1692 : !Convective moving mountain gravity waves (Beres scheme).
1693 : !------------------------------------------------------------------
1694 :
1695 0 : call outfld('U_MOVMTN_IN', u, ncol, lchnk)
1696 0 : call outfld('V_MOVMTN_IN', v, ncol, lchnk)
1697 :
1698 : ! Allocate wavenumber fields.
1699 0 : allocate(tau(ncol,-band_movmtn%ngwv:band_movmtn%ngwv,pver+1),stat=istat)
1700 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_movmtn%ngwv**2+1)*(pver+1))
1701 0 : allocate(gwut(ncol,pver,-band_movmtn%ngwv:band_movmtn%ngwv),stat=istat)
1702 0 : call alloc_err(istat,'gw_tend','gwut',ncol*pver*band_movmtn%ngwv**2+1)
1703 0 : allocate(phase_speeds(ncol,-band_movmtn%ngwv:band_movmtn%ngwv),stat=istat)
1704 0 : call alloc_err(istat,'gw_tend','phase_speeds',ncol*band_movmtn%ngwv**2+1)
1705 :
1706 : ! Set up heating
1707 0 : if (ttend_dp_idx > 0) then
1708 0 : call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)
1709 : else
1710 0 : allocate(ttend_dp(pcols,pver), stat=istat)
1711 0 : call alloc_err(istat, 'gw_tend', 'ttend_dp', pcols*pver)
1712 0 : ttend_dp = 0.0_r8
1713 : end if
1714 :
1715 : ! New couplings from CLUBB
1716 0 : call pbuf_get_field(pbuf, ttend_clubb_idx, ttend_clubb)
1717 0 : call pbuf_get_field(pbuf, thlp2_clubb_gw_idx, thlp2_clubb_gw)
1718 0 : call pbuf_get_field(pbuf, wpthlp_clubb_gw_idx, wpthlp_clubb_gw)
1719 0 : call pbuf_get_field(pbuf, upwp_clubb_gw_idx, upwp_clubb_gw)
1720 0 : call pbuf_get_field(pbuf, vpwp_clubb_gw_idx, vpwp_clubb_gw)
1721 :
1722 0 : xpwp_clubb(:ncol,:) = sqrt( upwp_clubb_gw(:ncol,:)**2 + vpwp_clubb_gw(:ncol,:)**2 )
1723 :
1724 0 : effgw = 1._r8
1725 : call gw_movmtn_src(ncol, lchnk, band_movmtn , movmtn_desc, &
1726 0 : u, v, ttend_dp(:ncol,:), ttend_clubb(:ncol,:), xpwp_clubb(:ncol,:) , &
1727 : zm, alpha_gw_movmtn, src_level, tend_level, &
1728 : tau, ubm, ubi, xv, yv, &
1729 0 : phase_speeds, hdepth)
1730 : !-------------------------------------------------------------
1731 : ! gw_movmtn_src returns wave-relative wind profiles ubm,ubi
1732 : ! and unit vector components describing direction of wavevector
1733 : ! and application of wave-drag force. I believe correct setting
1734 : ! for c is c=0, since it is incorporated in ubm and (xv,yv)
1735 : !--------------------------------------------------------------
1736 :
1737 0 : call outfld('SRC_LEVEL_MOVMTN', real(src_level,r8), ncol, lchnk)
1738 0 : call outfld('TND_LEVEL_MOVMTN', real(tend_level,r8), ncol, lchnk)
1739 0 : call outfld('UBI_MOVMTN', ubi, ncol, lchnk)
1740 0 : call outfld('UBM_MOVMTN', ubm, ncol, lchnk)
1741 :
1742 : call gw_drag_prof(ncol, band_movmtn, p, src_level, tend_level, dt, &
1743 : t, vramp, &
1744 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
1745 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
1746 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
1747 0 : lapply_effgw_in=gw_apply_tndmax )
1748 :
1749 : ! Project stress into directional components.
1750 0 : taucd = calc_taucd(ncol, band_movmtn%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
1751 :
1752 : ! add the diffusion coefficients
1753 0 : do k = 1, pver+1
1754 0 : egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
1755 : end do
1756 :
1757 : ! Store constituents tendencies
1758 0 : do m=1, pcnst
1759 0 : do k = 1, pver
1760 0 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
1761 : end do
1762 : end do
1763 :
1764 : ! Add the momentum tendencies to the output tendency arrays.
1765 0 : do k = 1, pver
1766 0 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
1767 0 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
1768 : end do
1769 :
1770 0 : do k = 1, pver
1771 0 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
1772 : end do
1773 :
1774 0 : call outfld('TAU_MOVMTN', tau(:,0,:), ncol, lchnk)
1775 0 : call outfld('GWUT_MOVMTN', gwut(:,:,0), ncol, lchnk)
1776 0 : call outfld('VTGW_MOVMTN', vtgw, ncol, lchnk)
1777 0 : call outfld('UTGW_MOVMTN', utgw, ncol, lchnk)
1778 0 : call outfld('HDEPTH_MOVMTN', hdepth/1000._r8, ncol, lchnk)
1779 0 : call outfld('NETDT_MOVMTN', ttend_dp, pcols, lchnk)
1780 0 : call outfld('TTEND_CLUBB', ttend_clubb, pcols, lchnk)
1781 0 : call outfld('THLP2_CLUBB_GW', thlp2_clubb_gw, pcols, lchnk)
1782 0 : call outfld('WPTHLP_CLUBB_GW', wpthlp_clubb_gw, pcols, lchnk)
1783 0 : call outfld('UPWP_CLUBB_GW', upwp_clubb_gw, pcols, lchnk)
1784 0 : call outfld('VPWP_CLUBB_GW', vpwp_clubb_gw, pcols, lchnk)
1785 :
1786 : !Deallocate variables that are no longer used:
1787 0 : deallocate(tau, gwut, phase_speeds)
1788 :
1789 : !Deallocate/nullify ttend_dp if not a pbuf variable:
1790 0 : if (ttend_dp_idx <= 0) then
1791 0 : deallocate(ttend_dp)
1792 : nullify(ttend_dp)
1793 : end if
1794 :
1795 : end if
1796 :
1797 1489176 : if (use_gw_convect_dp) then
1798 : !------------------------------------------------------------------
1799 : ! Convective gravity waves (Beres scheme, deep).
1800 : !------------------------------------------------------------------
1801 :
1802 : ! Allocate wavenumber fields.
1803 0 : allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1),stat=istat)
1804 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1)*(pver+1))
1805 0 : allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv),stat=istat)
1806 0 : call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_mid%ngwv**2+1))
1807 0 : allocate(phase_speeds(ncol,-band_mid%ngwv:band_mid%ngwv),stat=istat)
1808 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1))
1809 :
1810 : ! Set up heating
1811 0 : call pbuf_get_field(pbuf, ttend_dp_idx, ttend_dp)
1812 :
1813 : ! Efficiency of gravity wave momentum transfer.
1814 : ! This is really only to remove the pole points.
1815 0 : where (pi/2._r8 - abs(state1%lat(:ncol)) >= 4*epsilon(1._r8))
1816 : effgw = effgw_beres_dp
1817 : elsewhere
1818 : effgw = 0._r8
1819 : end where
1820 :
1821 : ! Determine wave sources for Beres deep scheme
1822 : call gw_beres_src(ncol, band_mid, beres_dp_desc, &
1823 0 : u, v, ttend_dp(:ncol,:), zm, src_level, tend_level, tau, &
1824 0 : ubm, ubi, xv, yv, phase_speeds, hdepth, maxq0)
1825 :
1826 : ! Solve for the drag profile with Beres source spectrum.
1827 : call gw_drag_prof(ncol, band_mid, p, src_level, tend_level, dt, &
1828 : t, vramp, &
1829 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
1830 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
1831 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
1832 0 : lapply_effgw_in=gw_apply_tndmax)
1833 :
1834 : ! Project stress into directional components.
1835 0 : taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
1836 :
1837 : ! add the diffusion coefficients
1838 0 : do k = 1, pver+1
1839 0 : egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
1840 : end do
1841 :
1842 : ! Store constituents tendencies
1843 0 : do m=1, pcnst
1844 0 : do k = 1, pver
1845 0 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
1846 : end do
1847 : end do
1848 :
1849 : ! Find momentum flux, and use it to fix the wind tendencies below
1850 : ! the gravity wave region.
1851 0 : call momentum_flux(tend_level, taucd, um_flux, vm_flux)
1852 0 : call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
1853 :
1854 : ! Add the momentum tendencies to the output tendency arrays.
1855 0 : do k = 1, pver
1856 0 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
1857 0 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
1858 : end do
1859 :
1860 : ! Find energy change in the current state, and use fixer to apply
1861 : ! the difference in lower levels.
1862 0 : call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
1863 0 : ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
1864 0 : call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)
1865 :
1866 0 : do k = 1, pver
1867 0 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
1868 : end do
1869 :
1870 : ! Change ttgw to a temperature tendency before outputing it.
1871 0 : ttgw = ttgw / cpair
1872 : call gw_spec_outflds(beres_dp_pf, lchnk, ncol, band_mid, phase_speeds, u, v, &
1873 0 : xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
1874 0 : taucd)
1875 :
1876 : ! Diagnostic outputs (convert hdepth to km).
1877 0 : call outfld('NETDT', ttend_dp, pcols, lchnk)
1878 0 : call outfld('HDEPTH', hdepth/1000._r8, ncol, lchnk)
1879 0 : call outfld('MAXQ0', maxq0, ncol, lchnk)
1880 :
1881 0 : deallocate(tau, gwut, phase_speeds)
1882 :
1883 : end if
1884 :
1885 1489176 : if (use_gw_convect_sh) then
1886 : !------------------------------------------------------------------
1887 : ! Convective gravity waves (Beres scheme, shallow).
1888 : !------------------------------------------------------------------
1889 :
1890 : ! Allocate wavenumber fields.
1891 0 : allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1),stat=istat)
1892 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1)*(pver+1))
1893 0 : allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv),stat=istat)
1894 0 : call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_mid%ngwv**2+1))
1895 0 : allocate(phase_speeds(ncol,-band_mid%ngwv:band_mid%ngwv),stat=istat)
1896 0 : call alloc_err(istat,'gw_tend','phase_speeds',ncol*(band_mid%ngwv**2+1))
1897 :
1898 : ! Set up heating
1899 0 : call pbuf_get_field(pbuf, ttend_sh_idx, ttend_sh)
1900 :
1901 : ! Efficiency of gravity wave momentum transfer.
1902 : ! This is really only to remove the pole points.
1903 0 : where (pi/2._r8 - abs(state1%lat(:ncol)) >= 4*epsilon(1._r8))
1904 : effgw = effgw_beres_sh
1905 : elsewhere
1906 : effgw = 0._r8
1907 : end where
1908 :
1909 : ! Determine wave sources for Beres shallow scheme
1910 : call gw_beres_src(ncol, band_mid, beres_sh_desc, &
1911 0 : u, v, ttend_sh(:ncol,:), zm, src_level, tend_level, tau, &
1912 0 : ubm, ubi, xv, yv, phase_speeds, hdepth, maxq0)
1913 :
1914 : ! Solve for the drag profile with Beres source spectrum.
1915 : call gw_drag_prof(ncol, band_mid, p, src_level, tend_level, dt, &
1916 : t, vramp, &
1917 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
1918 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
1919 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
1920 0 : lapply_effgw_in=gw_apply_tndmax)
1921 :
1922 : ! Project stress into directional components.
1923 0 : taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
1924 :
1925 : ! add the diffusion coefficients
1926 0 : do k = 1, pver+1
1927 0 : egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
1928 : end do
1929 :
1930 : ! Store constituents tendencies
1931 0 : do m=1, pcnst
1932 0 : do k = 1, pver
1933 0 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
1934 : end do
1935 : end do
1936 :
1937 : ! Add the momentum tendencies to the output tendency arrays.
1938 : ! Don't calculate fixers, since we are too close to the ground to
1939 : ! spread momentum/energy differences across low layers.
1940 0 : do k = 1, pver
1941 0 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
1942 0 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
1943 0 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
1944 : end do
1945 :
1946 : ! Calculate energy change for output to CAM's energy checker.
1947 : ! This is sort of cheating; we don't have a good a priori idea of the
1948 : ! energy coming from surface stress, so we just integrate what we and
1949 : ! actually have so far and overwrite flx_heat with that.
1950 0 : call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
1951 0 : ptend%v(:ncol,:), ptend%s(:ncol,:), de)
1952 0 : flx_heat(:ncol) = de
1953 :
1954 : ! Change ttgw to a temperature tendency before outputing it.
1955 0 : ttgw = ttgw / cpair
1956 : call gw_spec_outflds(beres_sh_pf, lchnk, ncol, band_mid, phase_speeds, u, v, &
1957 0 : xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
1958 0 : taucd)
1959 :
1960 : ! Diagnostic outputs (convert SHDEPTH to km).
1961 0 : call outfld ('SNETDT', ttend_sh, pcols, lchnk)
1962 0 : call outfld ('SHDEPTH', hdepth/1000._r8, ncol, lchnk)
1963 0 : call outfld ('SMAXQ0', maxq0, ncol, lchnk)
1964 :
1965 0 : deallocate(tau, gwut, phase_speeds)
1966 :
1967 : end if
1968 :
1969 1489176 : if (use_gw_front .or. use_gw_front_igw) then
1970 : ! Get frontogenesis physics buffer fields set by dynamics.
1971 0 : call pbuf_get_field(pbuf, frontgf_idx, frontgf)
1972 0 : call pbuf_get_field(pbuf, frontga_idx, frontga)
1973 :
1974 : ! Output for diagnostics.
1975 0 : call outfld ('FRONTGF', frontgf, pcols, lchnk)
1976 0 : call outfld ('FRONTGFA', frontga, pcols, lchnk)
1977 : end if
1978 :
1979 1489176 : if (use_gw_front) then
1980 : !------------------------------------------------------------------
1981 : ! Frontally generated gravity waves
1982 : !------------------------------------------------------------------
1983 :
1984 : ! Allocate wavenumber fields.
1985 0 : allocate(tau(ncol,-band_mid%ngwv:band_mid%ngwv,pver+1),stat=istat)
1986 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1)*(pver+1))
1987 0 : allocate(gwut(ncol,pver,-band_mid%ngwv:band_mid%ngwv),stat=istat)
1988 0 : call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_mid%ngwv**2+1))
1989 0 : allocate(phase_speeds(ncol,-band_mid%ngwv:band_mid%ngwv),stat=istat)
1990 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_mid%ngwv**2+1))
1991 :
1992 : ! Efficiency of gravity wave momentum transfer.
1993 0 : effgw = effgw_cm
1994 : ! Frontogenesis is too high at the poles (at least for the FV
1995 : ! dycore), so introduce a polar taper.
1996 0 : if (gw_polar_taper) effgw = effgw * cos(state1%lat(:ncol))
1997 :
1998 : ! Determine the wave source for C&M background spectrum
1999 0 : call gw_cm_src(ncol, band_mid, cm_desc, u, v, frontgf(:ncol,:), &
2000 0 : src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
2001 :
2002 : ! Solve for the drag profile with C&M source spectrum.
2003 : call gw_drag_prof(ncol, band_mid, p, src_level, tend_level, dt, &
2004 : t, vramp, &
2005 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
2006 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
2007 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
2008 0 : lapply_effgw_in=gw_apply_tndmax)
2009 :
2010 : ! Project stress into directional components.
2011 0 : taucd = calc_taucd(ncol, band_mid%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
2012 :
2013 : ! add the diffusion coefficients
2014 0 : do k = 1, pver+1
2015 0 : egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
2016 : end do
2017 :
2018 : !Add the constituent tendencies
2019 0 : do m=1, pcnst
2020 0 : do k = 1, pver
2021 0 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
2022 : end do
2023 : end do
2024 :
2025 : ! Find momentum flux, and use it to fix the wind tendencies below
2026 : ! the gravity wave region.
2027 0 : call momentum_flux(tend_level, taucd, um_flux, vm_flux)
2028 0 : call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
2029 :
2030 : ! add the momentum tendencies to the output tendency arrays
2031 0 : do k = 1, pver
2032 0 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
2033 0 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
2034 : end do
2035 :
2036 : ! Find energy change in the current state, and use fixer to apply
2037 : ! the difference in lower levels.
2038 0 : call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
2039 0 : ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
2040 0 : call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)
2041 :
2042 0 : do k = 1, pver
2043 0 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
2044 : end do
2045 :
2046 : ! Change ttgw to a temperature tendency before outputing it.
2047 0 : ttgw = ttgw / cpair
2048 : call gw_spec_outflds(cm_pf, lchnk, ncol, band_mid, phase_speeds, u, v, &
2049 0 : xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
2050 0 : taucd)
2051 :
2052 0 : deallocate(tau, gwut, phase_speeds)
2053 :
2054 : end if
2055 :
2056 1489176 : if (use_gw_front_igw) then
2057 : !------------------------------------------------------------------
2058 : ! Frontally generated inertial gravity waves
2059 : !------------------------------------------------------------------
2060 :
2061 : ! Allocate wavenumber fields.
2062 0 : allocate(tau(ncol,-band_long%ngwv:band_long%ngwv,pver+1),stat=istat)
2063 0 : call alloc_err(istat,'gw_tend','tau',ncol*(band_long%ngwv**2+1)*(pver+1))
2064 0 : allocate(gwut(ncol,pver,-band_long%ngwv:band_long%ngwv),stat=istat)
2065 0 : call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_long%ngwv**2+1))
2066 0 : allocate(phase_speeds(ncol,-band_long%ngwv:band_long%ngwv),stat=istat)
2067 0 : call alloc_err(istat,'gw_tend','phase_speeds',ncol*(band_long%ngwv**2+1))
2068 0 : allocate(ro_adjust(ncol,-band_long%ngwv:band_long%ngwv,pver+1),stat=istat)
2069 0 : call alloc_err(istat,'gw_tend','ro_adjust',ncol*(band_long%ngwv**2+1)*(pver+1))
2070 :
2071 : ! Efficiency of gravity wave momentum transfer.
2072 0 : effgw = effgw_cm_igw
2073 :
2074 : ! Frontogenesis is too high at the poles (at least for the FV
2075 : ! dycore), so introduce a polar taper.
2076 0 : if (gw_polar_taper) then
2077 0 : where (abs(state1%lat(:ncol)) <= 89._r8*degree2radian)
2078 : effgw = effgw * 0.25_r8 * &
2079 0 : (1._r8+tanh((state1%lat(:ncol)+al0)/dlat0)) * &
2080 : (1._r8-tanh((state1%lat(:ncol)-al0)/dlat0))
2081 : elsewhere
2082 : effgw = 0._r8
2083 : end where
2084 : end if
2085 :
2086 : ! Determine the wave source for C&M background spectrum
2087 0 : call gw_cm_src(ncol, band_long, cm_igw_desc, u, v, frontgf(:ncol,:), &
2088 0 : src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
2089 :
2090 : call adjust_inertial(band_long, tend_level, u_coriolis, phase_speeds, ubi, &
2091 0 : tau, ro_adjust)
2092 :
2093 : ! Solve for the drag profile with C&M source spectrum.
2094 : call gw_drag_prof(ncol, band_long, p, src_level, tend_level, dt, &
2095 : t, vramp, &
2096 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
2097 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
2098 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, ro_adjust=ro_adjust, &
2099 0 : lapply_effgw_in=gw_apply_tndmax)
2100 :
2101 : ! Project stress into directional components.
2102 0 : taucd = calc_taucd(ncol, band_long%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi)
2103 :
2104 : ! add the diffusion coefficients
2105 0 : do k = 1, pver+1
2106 0 : egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
2107 : end do
2108 :
2109 : !Add the constituent tendencies
2110 0 : do m=1, pcnst
2111 0 : do k = 1, pver
2112 0 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
2113 : end do
2114 : end do
2115 :
2116 : ! Find momentum flux, and use it to fix the wind tendencies below
2117 : ! the gravity wave region.
2118 0 : call momentum_flux(tend_level, taucd, um_flux, vm_flux)
2119 0 : call momentum_fixer(tend_level, p, um_flux, vm_flux, utgw, vtgw)
2120 :
2121 : ! add the momentum tendencies to the output tendency arrays
2122 0 : do k = 1, pver
2123 0 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
2124 0 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
2125 : end do
2126 :
2127 : ! Find energy change in the current state, and use fixer to apply
2128 : ! the difference in lower levels.
2129 0 : call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
2130 0 : ptend%v(:ncol,:), ptend%s(:ncol,:)+ttgw, de)
2131 0 : call energy_fixer(tend_level, p, de-flx_heat(:ncol), ttgw)
2132 :
2133 0 : do k = 1, pver
2134 0 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
2135 : end do
2136 :
2137 : ! Change ttgw to a temperature tendency before outputing it.
2138 0 : ttgw = ttgw / cpair
2139 : call gw_spec_outflds(cm_igw_pf, lchnk, ncol, band_long, phase_speeds, u, v, &
2140 0 : xv, yv, gwut, dttdf, dttke, tau(:,:,2:), utgw, vtgw, ttgw, &
2141 0 : taucd)
2142 :
2143 0 : deallocate(tau, gwut, phase_speeds, ro_adjust)
2144 :
2145 : end if
2146 :
2147 1489176 : if (use_gw_oro) then
2148 : !---------------------------------------------------------------------
2149 : ! Orographic stationary gravity waves
2150 : !---------------------------------------------------------------------
2151 :
2152 : ! Allocate wavenumber fields.
2153 5956704 : allocate(tau(ncol,band_oro%ngwv:band_oro%ngwv,pver+1),stat=istat)
2154 1489176 : call alloc_err(istat,'gw_tend','tau',ncol*(band_oro%ngwv**2+1)*(pver+1))
2155 4467528 : allocate(gwut(ncol,pver,band_oro%ngwv:band_oro%ngwv),stat=istat)
2156 1489176 : call alloc_err(istat,'gw_tend','gwut',ncol*pver*(band_oro%ngwv**2+1))
2157 4467528 : allocate(phase_speeds(ncol,band_oro%ngwv:band_oro%ngwv),stat=istat)
2158 1489176 : call alloc_err(istat,'gw_tend','phase_speeds',ncol*(band_oro%ngwv**2+1))
2159 :
2160 : ! Efficiency of gravity wave momentum transfer.
2161 : ! Take into account that wave sources are only over land.
2162 1489176 : call pbuf_get_field(pbuf, sgh_idx, sgh)
2163 :
2164 1489176 : if (gw_lndscl_sgh) then
2165 121350528 : where (cam_in%landfrac(:ncol) >= epsilon(1._r8))
2166 1489176 : effgw = effgw_oro * cam_in%landfrac(:ncol)
2167 1489176 : sgh_scaled = sgh(:ncol) / sqrt(cam_in%landfrac(:ncol))
2168 : elsewhere
2169 : effgw = 0._r8
2170 : sgh_scaled = 0._r8
2171 : end where
2172 :
2173 : ! Determine the orographic wave source
2174 : call gw_oro_src(ncol, band_oro, p, &
2175 : u, v, t, sgh_scaled, zm, nm, &
2176 1489176 : src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
2177 : else
2178 0 : effgw = effgw_oro
2179 :
2180 : ! Determine the orographic wave source
2181 : call gw_oro_src(ncol, band_oro, p, &
2182 0 : u, v, t, sgh(:ncol), zm, nm, &
2183 0 : src_level, tend_level, tau, ubm, ubi, xv, yv, phase_speeds)
2184 : endif
2185 24865776 : do i = 1, ncol
2186 24865776 : if (state1%lat(i) < 0._r8) then
2187 642856500 : tau(i,:,:) = tau(i,:,:) * gw_oro_south_fac
2188 : end if
2189 : end do
2190 :
2191 : ! Solve for the drag profile with orographic sources.
2192 : call gw_drag_prof(ncol, band_oro, p, src_level, tend_level, dt, &
2193 : t, vramp, &
2194 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
2195 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
2196 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
2197 1489176 : lapply_effgw_in=gw_apply_tndmax)
2198 :
2199 : ! For orographic waves, don't bother with taucd, since there are no
2200 : ! momentum conservation routines or directional diagnostics.
2201 :
2202 : ! add the diffusion coefficients
2203 41696928 : do k = 1, pver+1
2204 672865128 : egwdffi_tot(:,k) = egwdffi_tot(:,k) + egwdffi(:,k)
2205 : end do
2206 :
2207 : ! Add the orographic tendencies to the spectrum tendencies.
2208 : ! Don't calculate fixers, since we are too close to the ground to
2209 : ! spread momentum/energy differences across low layers.
2210 40207752 : do k = 1, pver
2211 646510176 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
2212 646510176 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
2213 646510176 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
2214 : ! Convert to temperature tendency for output.
2215 647999352 : ttgw(:,k) = ttgw(:,k) / cpairv(:ncol, k, lchnk)
2216 : end do
2217 :
2218 : ! Calculate energy change for output to CAM's energy checker.
2219 : ! This is sort of cheating; we don't have a good a priori idea of the
2220 : ! energy coming from surface stress, so we just integrate what we and
2221 : ! actually have so far and overwrite flx_heat with that.
2222 0 : call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
2223 1489176 : ptend%v(:ncol,:), ptend%s(:ncol,:), de)
2224 24865776 : flx_heat(:ncol) = de
2225 :
2226 5956704 : do m = 1, pcnst
2227 122112432 : do k = 1, pver
2228 1943998056 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
2229 : end do
2230 : end do
2231 :
2232 : ! Write output fields to history file
2233 1489176 : call outfld('TAUAORO', tau(:,0,:), ncol, lchnk)
2234 1489176 : call outfld('UTGWORO', utgw, ncol, lchnk)
2235 1489176 : call outfld('VTGWORO', vtgw, ncol, lchnk)
2236 1489176 : call outfld('TTGWORO', ttgw, ncol, lchnk)
2237 647999352 : call outfld('TTGWSDFORO', dttdf / cpair, ncol, lchnk)
2238 647999352 : call outfld('TTGWSKEORO', dttke / cpair, ncol, lchnk)
2239 24865776 : tau0x = tau(:,0,pver+1) * xv
2240 24865776 : tau0y = tau(:,0,pver+1) * yv
2241 1489176 : call outfld('TAUGWX', tau0x, ncol, lchnk)
2242 1489176 : call outfld('TAUGWY', tau0y, ncol, lchnk)
2243 :
2244 1489176 : deallocate(tau, gwut, phase_speeds)
2245 :
2246 : end if
2247 :
2248 1489176 : if (use_gw_rdg_beta) then
2249 : !---------------------------------------------------------------------
2250 : ! Orographic stationary gravity waves
2251 : !---------------------------------------------------------------------
2252 :
2253 : ! Efficiency of gravity wave momentum transfer.
2254 : ! Take into account that wave sources are only over land.
2255 :
2256 0 : hwdth => rdg_hwdth(:ncol,:,lchnk)
2257 0 : clngt => rdg_clngt(:ncol,:,lchnk)
2258 0 : gbxar => rdg_gbxar(:ncol,lchnk)
2259 0 : mxdis => rdg_mxdis(:ncol,:,lchnk)
2260 0 : angll => rdg_angll(:ncol,:,lchnk)
2261 0 : anixy => rdg_anixy(:ncol,:,lchnk)
2262 :
2263 0 : where(mxdis < 0._r8)
2264 : mxdis = 0._r8
2265 : end where
2266 :
2267 : ! Save state at top of routine
2268 : ! Useful for unit testing checks
2269 0 : call outfld('UEGW', u , ncol, lchnk)
2270 0 : call outfld('VEGW', v , ncol, lchnk)
2271 0 : call outfld('TEGW', t , ncol, lchnk)
2272 0 : call outfld('ZEGW', zi , ncol, lchnk)
2273 0 : call outfld('ZMGW', zm , ncol, lchnk)
2274 :
2275 : call gw_rdg_calc(&
2276 : 'BETA ', ncol, lchnk, n_rdg_beta, dt, &
2277 : u, v, t, p, piln, zm, zi, &
2278 : nm, ni, rhoi, kvtt, q, dse, &
2279 : effgw_rdg_beta, effgw_rdg_beta_max, &
2280 : hwdth, clngt, gbxar, mxdis, angll, anixy, &
2281 : rdg_beta_cd_llb, trpd_leewv_rdg_beta, &
2282 0 : ptend, flx_heat)
2283 :
2284 : end if
2285 :
2286 1489176 : if (use_gw_rdg_gamma) then
2287 : !---------------------------------------------------------------------
2288 : ! Orographic stationary gravity waves
2289 : !---------------------------------------------------------------------
2290 :
2291 : ! Efficiency of gravity wave momentum transfer.
2292 : ! Take into account that wave sources are only over land.
2293 :
2294 0 : hwdthg => rdg_hwdthg(:ncol,:,lchnk)
2295 0 : clngtg => rdg_clngtg(:ncol,:,lchnk)
2296 0 : gbxar => rdg_gbxar(:ncol,lchnk)
2297 0 : mxdisg => rdg_mxdisg(:ncol,:,lchnk)
2298 0 : angllg => rdg_angllg(:ncol,:,lchnk)
2299 0 : anixyg => rdg_anixyg(:ncol,:,lchnk)
2300 :
2301 0 : where(mxdisg < 0._r8)
2302 : mxdisg = 0._r8
2303 : end where
2304 :
2305 : call gw_rdg_calc(&
2306 : 'GAMMA', ncol, lchnk, n_rdg_gamma, dt, &
2307 : u, v, t, p, piln, zm, zi, &
2308 : nm, ni, rhoi, kvtt, q, dse, &
2309 : effgw_rdg_gamma, effgw_rdg_gamma_max, &
2310 : hwdthg, clngtg, gbxar, mxdisg, angllg, anixyg, &
2311 : rdg_gamma_cd_llb, trpd_leewv_rdg_gamma, &
2312 0 : ptend, flx_heat)
2313 :
2314 : endif
2315 :
2316 : ! Convert the tendencies for the dry constituents to dry air basis.
2317 5956704 : do m = 1, pcnst
2318 5956704 : if (cnst_type(m).eq.'dry') then
2319 0 : do k = 1, pver
2320 0 : do i = 1, ncol
2321 0 : ptend%q(i,k,m) = ptend%q(i,k,m)*state1%pdel(i,k)/state1%pdeldry(i,k)
2322 : end do
2323 : end do
2324 : end if
2325 : end do
2326 :
2327 : ! Write totals to history file.
2328 1489176 : call outfld('EKGW', egwdffi_tot , ncol, lchnk)
2329 659704968 : call outfld('TTGW', ptend%s/cpairv(:,:,lchnk), pcols, lchnk)
2330 :
2331 1489176 : call outfld('UTGW_TOTAL', ptend%u, pcols, lchnk)
2332 1489176 : call outfld('VTGW_TOTAL', ptend%v, pcols, lchnk)
2333 :
2334 1489176 : call outfld('QTGW', ptend%q(:,:,1), pcols, lchnk)
2335 1489176 : call outfld('CLDLIQTGW', ptend%q(:,:,ixcldliq), pcols, lchnk)
2336 1489176 : call outfld('CLDICETGW', ptend%q(:,:,ixcldice), pcols, lchnk)
2337 :
2338 : ! Destroy objects.
2339 1489176 : call p%finalize()
2340 :
2341 2978352 : end subroutine gw_tend
2342 :
2343 : !==========================================================================
2344 :
2345 0 : subroutine gw_rdg_calc( &
2346 0 : type, ncol, lchnk, n_rdg, dt, &
2347 0 : u, v, t, p, piln, zm, zi, &
2348 0 : nm, ni, rhoi, kvtt, q, dse, &
2349 : effgw_rdg, effgw_rdg_max, &
2350 0 : hwdth, clngt, gbxar, &
2351 0 : mxdis, angll, anixy, &
2352 : rdg_cd_llb, trpd_leewv, &
2353 : ptend, flx_heat)
2354 :
2355 1489176 : use coords_1d, only: Coords1D
2356 : use gw_rdg, only: gw_rdg_src, gw_rdg_belowpeak, gw_rdg_break_trap, gw_rdg_do_vdiff
2357 : use gw_common, only: gw_drag_prof, energy_change
2358 :
2359 : character(len=5), intent(in) :: type ! BETA or GAMMA
2360 : integer, intent(in) :: ncol ! number of atmospheric columns
2361 : integer, intent(in) :: lchnk ! chunk identifier
2362 : integer, intent(in) :: n_rdg
2363 : real(r8), intent(in) :: dt ! Time step.
2364 :
2365 : real(r8), intent(in) :: u(ncol,pver) ! Midpoint zonal winds. ( m s-1)
2366 : real(r8), intent(in) :: v(ncol,pver) ! Midpoint meridional winds. ( m s-1)
2367 : real(r8), intent(in) :: t(ncol,pver) ! Midpoint temperatures. (K)
2368 : type(Coords1D), intent(in) :: p ! Pressure coordinates.
2369 : real(r8), intent(in) :: piln(ncol,pver+1) ! Log of interface pressures.
2370 : real(r8), intent(in) :: zm(ncol,pver) ! Midpoint altitudes above ground (m).
2371 : real(r8), intent(in) :: zi(ncol,pver+1) ! Interface altitudes above ground (m).
2372 : real(r8), intent(in) :: nm(ncol,pver) ! Midpoint Brunt-Vaisalla frequencies (s-1).
2373 : real(r8), intent(in) :: ni(ncol,pver+1) ! Interface Brunt-Vaisalla frequencies (s-1).
2374 : real(r8), intent(in) :: rhoi(ncol,pver+1) ! Interface density (kg m-3).
2375 : real(r8), intent(in) :: kvtt(ncol,pver+1) ! Molecular thermal diffusivity.
2376 : real(r8), intent(in) :: q(:,:,:) ! Constituent array.
2377 : real(r8), intent(in) :: dse(ncol,pver) ! Dry static energy.
2378 :
2379 :
2380 : real(r8), intent(in) :: effgw_rdg ! Tendency efficiency.
2381 : real(r8), intent(in) :: effgw_rdg_max
2382 : real(r8), intent(in) :: hwdth(ncol,prdg) ! width of ridges.
2383 : real(r8), intent(in) :: clngt(ncol,prdg) ! length of ridges.
2384 : real(r8), intent(in) :: gbxar(ncol) ! gridbox area
2385 :
2386 : real(r8), intent(in) :: mxdis(ncol,prdg) ! Height estimate for ridge (m).
2387 : real(r8), intent(in) :: angll(ncol,prdg) ! orientation of ridges.
2388 : real(r8), intent(in) :: anixy(ncol,prdg) ! Anisotropy parameter.
2389 :
2390 : real(r8), intent(in) :: rdg_cd_llb ! Drag coefficient for low-level flow
2391 : logical, intent(in) :: trpd_leewv
2392 :
2393 : type(physics_ptend), intent(inout):: ptend ! Parameterization net tendencies.
2394 :
2395 : real(r8), intent(out) :: flx_heat(pcols)
2396 :
2397 : !---------------------------Local storage-------------------------------
2398 :
2399 : integer :: k, m, nn, istat
2400 :
2401 0 : real(r8), allocatable :: tau(:,:,:) ! wave Reynolds stress
2402 : ! gravity wave wind tendency for each wave
2403 0 : real(r8), allocatable :: gwut(:,:,:)
2404 : ! Wave phase speeds for each column
2405 0 : real(r8), allocatable :: phase_speeds(:,:)
2406 :
2407 : ! Isotropic source flag [anisotropic orography].
2408 0 : integer :: isoflag(ncol)
2409 :
2410 : ! horiz wavenumber [anisotropic orography].
2411 0 : real(r8) :: kwvrdg(ncol)
2412 :
2413 : ! Efficiency for a gravity wave source.
2414 0 : real(r8) :: effgw(ncol)
2415 :
2416 : ! Indices of top gravity wave source level and lowest level where wind
2417 : ! tendencies are allowed.
2418 0 : integer :: src_level(ncol)
2419 0 : integer :: tend_level(ncol)
2420 0 : integer :: bwv_level(ncol)
2421 0 : integer :: tlb_level(ncol)
2422 :
2423 : ! Projection of wind at midpoints and interfaces.
2424 0 : real(r8) :: ubm(ncol,pver)
2425 0 : real(r8) :: ubi(ncol,pver+1)
2426 :
2427 : ! Unit vectors of source wind (zonal and meridional components).
2428 0 : real(r8) :: xv(ncol)
2429 0 : real(r8) :: yv(ncol)
2430 :
2431 : ! Averages over source region.
2432 0 : real(r8) :: ubmsrc(ncol) ! On-ridge wind.
2433 0 : real(r8) :: usrc(ncol) ! Zonal wind.
2434 0 : real(r8) :: vsrc(ncol) ! Meridional wind.
2435 0 : real(r8) :: nsrc(ncol) ! B-V frequency.
2436 0 : real(r8) :: rsrc(ncol) ! Density.
2437 :
2438 : ! normalized wavenumber
2439 0 : real(r8) :: m2src(ncol)
2440 :
2441 : ! Top of low-level flow layer.
2442 0 : real(r8) :: tlb(ncol)
2443 :
2444 : ! Bottom of linear wave region.
2445 0 : real(r8) :: bwv(ncol)
2446 :
2447 : ! Froude numbers for flow/drag regimes
2448 0 : real(r8) :: Fr1(ncol)
2449 0 : real(r8) :: Fr2(ncol)
2450 0 : real(r8) :: Frx(ncol)
2451 :
2452 : ! Wave Reynolds stresses at source level
2453 0 : real(r8) :: tauoro(ncol)
2454 0 : real(r8) :: taudsw(ncol)
2455 :
2456 : ! Surface streamline displacement height for linear waves.
2457 0 : real(r8) :: hdspwv(ncol)
2458 :
2459 : ! Surface streamline displacement height for downslope wind regime.
2460 0 : real(r8) :: hdspdw(ncol)
2461 :
2462 : ! Wave breaking level
2463 0 : real(r8) :: wbr(ncol)
2464 :
2465 0 : real(r8) :: utgw(ncol,pver) ! zonal wind tendency
2466 0 : real(r8) :: vtgw(ncol,pver) ! meridional wind tendency
2467 0 : real(r8) :: ttgw(ncol,pver) ! temperature tendency
2468 0 : real(r8) :: qtgw(ncol,pver,pcnst) ! constituents tendencies
2469 :
2470 : ! Effective gravity wave diffusivity at interfaces.
2471 0 : real(r8) :: egwdffi(ncol,pver+1)
2472 :
2473 : ! Temperature tendencies from diffusion and kinetic energy.
2474 0 : real(r8) :: dttdf(ncol,pver)
2475 0 : real(r8) :: dttke(ncol,pver)
2476 :
2477 : ! Wave stress in zonal/meridional direction
2478 0 : real(r8) :: taurx(ncol,pver+1)
2479 0 : real(r8) :: taurx0(ncol,pver+1)
2480 0 : real(r8) :: taury(ncol,pver+1)
2481 0 : real(r8) :: taury0(ncol,pver+1)
2482 : ! Provisional absolute wave stress from gw_drag_prof
2483 0 : real(r8) :: tau_diag(ncol,pver+1)
2484 :
2485 : ! U,V tendency accumulators
2486 0 : real(r8) :: utrdg(ncol,pver)
2487 0 : real(r8) :: vtrdg(ncol,pver)
2488 0 : real(r8) :: ttrdg(ncol,pver)
2489 :
2490 : ! Energy change used by fixer.
2491 0 : real(r8) :: de(ncol)
2492 :
2493 : character(len=1) :: cn
2494 : character(len=9) :: fname(4)
2495 : !----------------------------------------------------------------------------
2496 :
2497 : ! Allocate wavenumber fields.
2498 0 : allocate(tau(ncol,band_oro%ngwv:band_oro%ngwv,pver+1),stat=istat)
2499 0 : call alloc_err(istat,'gw_rdg_calc','tau',ncol*(band_oro%ngwv**2+1)*(pver+1))
2500 0 : allocate(gwut(ncol,pver,band_oro%ngwv:band_oro%ngwv),stat=istat)
2501 0 : call alloc_err(istat,'rdg_calc','gwut',ncol*pver*(band_oro%ngwv**2+1))
2502 0 : allocate(phase_speeds(ncol,band_oro%ngwv:band_oro%ngwv),stat=istat)
2503 0 : call alloc_err(istat,'rdg_calc','phase_speeds',ncol*(band_oro%ngwv**2+1))
2504 :
2505 : ! initialize accumulated momentum fluxes and tendencies
2506 0 : taurx = 0._r8
2507 0 : taury = 0._r8
2508 0 : ttrdg = 0._r8
2509 0 : utrdg = 0._r8
2510 0 : vtrdg = 0._r8
2511 0 : tau_diag = -9999._r8
2512 :
2513 0 : do nn = 1, n_rdg
2514 0 : kwvrdg = 0.001_r8 / ( hwdth(:,nn) + 0.001_r8 ) ! this cant be done every time step !!!
2515 0 : isoflag = 0
2516 0 : effgw = effgw_rdg * ( hwdth(1:ncol,nn)* clngt(1:ncol,nn) ) / gbxar(1:ncol)
2517 0 : effgw = min( effgw_rdg_max , effgw )
2518 :
2519 : call gw_rdg_src(ncol, band_oro, p, &
2520 : u, v, t, mxdis(:,nn), angll(:,nn), anixy(:,nn), kwvrdg, isoflag, zi, nm, &
2521 : src_level, tend_level, bwv_level, tlb_level, tau, ubm, ubi, xv, yv, &
2522 0 : ubmsrc, usrc, vsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, phase_speeds)
2523 :
2524 : call gw_rdg_belowpeak(ncol, band_oro, rdg_cd_llb, &
2525 : t, mxdis(:,nn), anixy(:,nn), kwvrdg, &
2526 : zi, nm, ni, rhoi, &
2527 : src_level, tau, &
2528 : ubmsrc, nsrc, rsrc, m2src, tlb, bwv, Fr1, Fr2, Frx, &
2529 0 : tauoro, taudsw, hdspwv, hdspdw)
2530 :
2531 :
2532 : call gw_rdg_break_trap(ncol, band_oro, &
2533 : zi, nm, ni, ubm, ubi, rhoi, kwvrdg , bwv, tlb, wbr, &
2534 : src_level, tlb_level, hdspwv, hdspdw, mxdis(:,nn), &
2535 : tauoro, taudsw, tau, &
2536 0 : ldo_trapped_waves=trpd_leewv)
2537 :
2538 :
2539 : call gw_drag_prof(ncol, band_oro, p, src_level, tend_level, dt, &
2540 : t, vramp, &
2541 : piln, rhoi, nm, ni, ubm, ubi, xv, yv, &
2542 : effgw, phase_speeds, kvtt, q, dse, tau, utgw, vtgw, &
2543 : ttgw, qtgw, egwdffi, gwut, dttdf, dttke, &
2544 : kwvrdg=kwvrdg, &
2545 0 : satfac_in = 1._r8, lapply_vdiff=gw_rdg_do_vdiff , tau_diag=tau_diag )
2546 :
2547 : ! Add the tendencies from each ridge to the totals.
2548 0 : do k = 1, pver
2549 : ! diagnostics
2550 0 : utrdg(:,k) = utrdg(:,k) + utgw(:,k)
2551 0 : vtrdg(:,k) = vtrdg(:,k) + vtgw(:,k)
2552 0 : ttrdg(:,k) = ttrdg(:,k) + ttgw(:,k)
2553 : ! physics tendencies
2554 0 : ptend%u(:ncol,k) = ptend%u(:ncol,k) + utgw(:,k)
2555 0 : ptend%v(:ncol,k) = ptend%v(:ncol,k) + vtgw(:,k)
2556 0 : ptend%s(:ncol,k) = ptend%s(:ncol,k) + ttgw(:,k)
2557 : end do
2558 :
2559 0 : do m = 1, pcnst
2560 0 : do k = 1, pver
2561 0 : ptend%q(:ncol,k,m) = ptend%q(:ncol,k,m) + qtgw(:,k,m)
2562 : end do
2563 : end do
2564 :
2565 0 : do k = 1, pver+1
2566 0 : taurx0(:,k) = tau(:,0,k)*xv
2567 0 : taury0(:,k) = tau(:,0,k)*yv
2568 0 : taurx(:,k) = taurx(:,k) + taurx0(:,k)
2569 0 : taury(:,k) = taury(:,k) + taury0(:,k)
2570 : end do
2571 :
2572 0 : if (nn == 1) then
2573 0 : call outfld('BWV_HT1', bwv, ncol, lchnk)
2574 0 : call outfld('TLB_HT1', tlb, ncol, lchnk)
2575 0 : call outfld('WBR_HT1', wbr, ncol, lchnk)
2576 0 : call outfld('TAUDSW1', taudsw, ncol, lchnk)
2577 0 : call outfld('TAUORO1', tauoro, ncol, lchnk)
2578 0 : call outfld('UBMSRC1', ubmsrc, ncol, lchnk)
2579 0 : call outfld('USRC1', usrc, ncol, lchnk)
2580 0 : call outfld('VSRC1', vsrc, ncol, lchnk)
2581 0 : call outfld('NSRC1' , nsrc, ncol, lchnk)
2582 : ! Froude numbers
2583 0 : call outfld('Fr1_DIAG' , Fr1, ncol, lchnk)
2584 0 : call outfld('Fr2_DIAG' , Fr2, ncol, lchnk)
2585 0 : call outfld('Frx_DIAG' , Frx, ncol, lchnk)
2586 : ! Ridge quantities - don't change. Written for convenience
2587 0 : call outfld('MXDIS1' , mxdis(:,nn) , ncol, lchnk)
2588 0 : call outfld('ANGLL1' , angll(:,nn) , ncol, lchnk)
2589 0 : call outfld('ANIXY1' , anixy(:,nn) , ncol, lchnk)
2590 0 : call outfld('HWDTH1' , hwdth(:,nn) , ncol, lchnk)
2591 0 : call outfld('CLNGT1' , clngt(:,nn) , ncol, lchnk)
2592 0 : call outfld('GBXAR1' , gbxar , ncol, lchnk)
2593 0 : call outfld('TAUM1_DIAG' , tau_diag , ncol, lchnk)
2594 0 : call outfld('TAU1RDG'//trim(type)//'M', tau(:,0,:), ncol, lchnk)
2595 0 : call outfld('UBM1'//trim(type), ubm, ncol, lchnk)
2596 0 : call outfld('UBT1RDG'//trim(type), gwut, ncol, lchnk)
2597 : end if
2598 :
2599 0 : if (nn <= 6) then
2600 0 : write(cn, '(i1)') nn
2601 0 : call outfld('TAU'//cn//'RDG'//trim(type)//'X', taurx0, ncol, lchnk)
2602 0 : call outfld('TAU'//cn//'RDG'//trim(type)//'Y', taury0, ncol, lchnk)
2603 0 : call outfld('UT'//cn//'RDG'//trim(type), utgw, ncol, lchnk)
2604 0 : call outfld('VT'//cn//'RDG'//trim(type), vtgw, ncol, lchnk)
2605 : end if
2606 :
2607 : end do ! end of loop over multiple ridges
2608 :
2609 : ! Calculate energy change for output to CAM's energy checker.
2610 0 : call energy_change(dt, p, u, v, ptend%u(:ncol,:), &
2611 0 : ptend%v(:ncol,:), ptend%s(:ncol,:), de)
2612 0 : flx_heat(:ncol) = de
2613 :
2614 0 : call outfld('TAUARDG'//trim(type)//'X', taurx, ncol, lchnk)
2615 0 : call outfld('TAUARDG'//trim(type)//'Y', taury, ncol, lchnk)
2616 :
2617 0 : if (trim(type) == 'BETA') then
2618 0 : fname(1) = 'TAUGWX'
2619 0 : fname(2) = 'TAUGWY'
2620 0 : fname(3) = 'UTGWORO'
2621 0 : fname(4) = 'VTGWORO'
2622 0 : else if (trim(type) == 'GAMMA') then
2623 0 : fname(1) = 'TAURDGGMX'
2624 0 : fname(2) = 'TAURDGGMY'
2625 0 : fname(3) = 'UTRDGGM'
2626 0 : fname(4) = 'VTRDGGM'
2627 : else
2628 : call endrun('gw_rdg_calc: FATAL: type must be either BETA or GAMMA'&
2629 0 : //' type= '//type)
2630 : end if
2631 :
2632 0 : call outfld(fname(1), taurx(:,pver+1), ncol, lchnk)
2633 0 : call outfld(fname(2), taury(:,pver+1), ncol, lchnk)
2634 0 : call outfld(fname(3), utrdg, ncol, lchnk)
2635 0 : call outfld(fname(4), vtrdg, ncol, lchnk)
2636 0 : call outfld('TTGWORO', ttrdg / cpair, ncol, lchnk)
2637 :
2638 0 : deallocate(tau, gwut, phase_speeds)
2639 :
2640 0 : end subroutine gw_rdg_calc
2641 :
2642 : !==========================================================================
2643 :
2644 : ! Add all history fields for a gravity wave spectrum source.
2645 0 : subroutine gw_spec_addflds(prefix, scheme, band, history_defaults)
2646 : use cam_history, only: addfld, add_default, register_vector_field
2647 :
2648 : !------------------------------Arguments--------------------------------
2649 :
2650 : ! One character prefix prepended to output fields.
2651 : character(len=1), intent(in) :: prefix
2652 : ! Gravity wave scheme name prepended to output field descriptions.
2653 : character(len=*), intent(in) :: scheme
2654 : ! Wave speeds.
2655 : type(GWBand), intent(in) :: band
2656 : ! Whether or not to call add_default for fields output by WACCM.
2657 : logical, intent(in) :: history_defaults
2658 :
2659 : !---------------------------Local storage-------------------------------
2660 :
2661 : integer :: l
2662 : ! 7 chars is enough for "-100.00"
2663 : character(len=7) :: fnum
2664 : ! 10 chars is enough for "BTAUXSn32"
2665 : character(len=10) :: dumc1x, dumc1y
2666 : ! Allow 80 chars for description
2667 : character(len=80) dumc2
2668 :
2669 : !-----------------------------------------------------------------------
2670 :
2671 : ! Overall wind tendencies.
2672 : call addfld (trim(prefix)//'UTGWSPEC',(/ 'lev' /), 'A','m s-2', &
2673 0 : trim(scheme)//' U tendency - gravity wave spectrum')
2674 : call addfld (trim(prefix)//'VTGWSPEC',(/ 'lev' /), 'A','m s-2', &
2675 0 : trim(scheme)//' V tendency - gravity wave spectrum')
2676 0 : call register_vector_field(trim(prefix)//'UTGWSPEC',trim(prefix)//'VTGWSPEC')
2677 :
2678 : call addfld (trim(prefix)//'TTGWSPEC',(/ 'lev' /), 'A','K s-1', &
2679 0 : trim(scheme)//' T tendency - gravity wave spectrum')
2680 :
2681 : ! Wind tendencies broken across five spectral bins.
2682 : call addfld (trim(prefix)//'UTEND1', (/ 'lev' /), 'A','m s-2', &
2683 0 : trim(scheme)//' U tendency c < -40')
2684 : call addfld (trim(prefix)//'UTEND2', (/ 'lev' /), 'A','m s-2', &
2685 0 : trim(scheme)//' U tendency -40 < c < -15')
2686 : call addfld (trim(prefix)//'UTEND3', (/ 'lev' /), 'A','m s-2', &
2687 0 : trim(scheme)//' U tendency -15 < c < 15')
2688 : call addfld (trim(prefix)//'UTEND4', (/ 'lev' /), 'A','m s-2', &
2689 0 : trim(scheme)//' U tendency 15 < c < 40')
2690 : call addfld (trim(prefix)//'UTEND5', (/ 'lev' /), 'A','m s-2', &
2691 0 : trim(scheme)//' U tendency 40 < c ')
2692 :
2693 : ! Reynold's stress toward each cardinal direction, and net zonal stress.
2694 : call addfld (trim(prefix)//'TAUE' , (/ 'ilev' /), 'A','Pa', &
2695 0 : trim(scheme)//' Eastward Reynolds stress')
2696 : call addfld (trim(prefix)//'TAUW' , (/ 'ilev' /), 'A','Pa', &
2697 0 : trim(scheme)//' Westward Reynolds stress')
2698 : call addfld (trim(prefix)//'TAUNET' , (/ 'ilev' /), 'A','Pa', &
2699 0 : trim(scheme)//' E+W Reynolds stress')
2700 : call addfld (trim(prefix)//'TAUN' , (/ 'ilev' /), 'A','Pa', &
2701 0 : trim(scheme)//' Northward Reynolds stress')
2702 : call addfld (trim(prefix)//'TAUS' , (/ 'ilev' /), 'A','Pa', &
2703 0 : trim(scheme)//' Southward Reynolds stress')
2704 :
2705 : ! Momentum flux in each direction.
2706 : call addfld (trim(prefix)//'EMF', (/ 'lev' /), 'A','Pa', &
2707 0 : trim(scheme)//' Eastward MF')
2708 : call addfld (trim(prefix)//'WMF', (/ 'lev' /), 'A','Pa', &
2709 0 : trim(scheme)//' Westward MF')
2710 : call addfld (trim(prefix)//'NMF', (/ 'lev' /), 'A','Pa', &
2711 0 : trim(scheme)//' Northward MF')
2712 : call addfld (trim(prefix)//'SMF', (/ 'lev' /), 'A','Pa', &
2713 0 : trim(scheme)//' Southward MF')
2714 :
2715 : ! Temperature tendency terms.
2716 : call addfld (trim(prefix)//'TTGWSDF' , (/ 'lev' /), 'A','K s-1', &
2717 0 : trim(scheme)//' t tendency - diffusion term')
2718 : call addfld (trim(prefix)//'TTGWSKE' , (/ 'lev' /), 'A','K s-1', &
2719 0 : trim(scheme)//' t tendency - kinetic energy conversion term')
2720 :
2721 : ! Gravity wave source spectra by wave number.
2722 0 : do l=-band%ngwv,band%ngwv
2723 : ! String containing reference speed.
2724 0 : write (fnum,fmt='(f7.2)') band%cref(l)
2725 :
2726 0 : dumc1x = tau_fld_name(l, prefix, x_not_y=.true.)
2727 0 : dumc1y = tau_fld_name(l, prefix, x_not_y=.false.)
2728 0 : dumc2 = trim(scheme)//" tau at c= "//trim(fnum)//" m s-1"
2729 0 : call addfld (trim(dumc1x),(/ 'lev' /), 'A','Pa',dumc2)
2730 0 : call addfld (trim(dumc1y),(/ 'lev' /), 'A','Pa',dumc2)
2731 :
2732 : end do
2733 :
2734 0 : if (history_defaults) then
2735 0 : call add_default(trim(prefix)//'UTGWSPEC', 1, ' ')
2736 0 : call add_default(trim(prefix)//'VTGWSPEC', 1, ' ')
2737 0 : call add_default(trim(prefix)//'TTGWSPEC', 1, ' ')
2738 0 : call add_default(trim(prefix)//'TAUE', 1, ' ')
2739 0 : call add_default(trim(prefix)//'TAUW', 1, ' ')
2740 0 : call add_default(trim(prefix)//'TAUNET', 1, ' ')
2741 0 : call add_default(trim(prefix)//'TAUN', 1, ' ')
2742 0 : call add_default(trim(prefix)//'TAUS', 1, ' ')
2743 : end if
2744 :
2745 0 : end subroutine gw_spec_addflds
2746 :
2747 : !==========================================================================
2748 :
2749 : ! Outputs for spectral waves.
2750 0 : subroutine gw_spec_outflds(prefix, lchnk, ncol, band, phase_speeds, u, v, xv, yv, &
2751 0 : gwut, dttdf, dttke, tau, utgw, vtgw, ttgw, taucd)
2752 :
2753 0 : use gw_common, only: west, east, south, north
2754 :
2755 : ! One-character prefix prepended to output fields.
2756 : character(len=1), intent(in) :: prefix
2757 : ! Chunk and number of columns in the chunk.
2758 : integer, intent(in) :: lchnk
2759 : integer, intent(in) :: ncol
2760 : ! Wave speeds.
2761 : type(GWBand), intent(in) :: band
2762 : ! Wave phase speeds for each column.
2763 : real(r8), intent(in) :: phase_speeds(ncol,-band%ngwv:band%ngwv)
2764 : ! Winds at cell midpoints.
2765 : real(r8), intent(in) :: u(ncol,pver)
2766 : real(r8), intent(in) :: v(ncol,pver)
2767 : ! Unit vector in the direction of wind at source level.
2768 : real(r8), intent(in) :: xv(ncol)
2769 : real(r8), intent(in) :: yv(ncol)
2770 : ! Wind tendency for each wave.
2771 : real(r8), intent(in) :: gwut(ncol,pver,-band%ngwv:band%ngwv)
2772 : ! Temperature tendencies from diffusion and kinetic energy.
2773 : real(r8) :: dttdf(ncol,pver)
2774 : real(r8) :: dttke(ncol,pver)
2775 : ! Wave Reynolds stress.
2776 : real(r8), intent(in) :: tau(ncol,-band%ngwv:band%ngwv,pver)
2777 : ! Zonal and meridional total wind tendencies.
2778 : real(r8), intent(in) :: utgw(ncol,pver)
2779 : real(r8), intent(in) :: vtgw(ncol,pver)
2780 : ! Temperature tendencies.
2781 : real(r8), intent(in) :: ttgw(ncol,pver)
2782 : ! Reynolds stress for waves propagating in each cardinal direction.
2783 : real(r8), intent(in) :: taucd(ncol,pver+1,4)
2784 :
2785 : ! Indices
2786 : integer :: i, k, l
2787 0 : integer :: ix(ncol, -band%ngwv:band%ngwv), iy(ncol, -band%ngwv:band%ngwv)
2788 0 : integer :: iu(ncol), iv(ncol)
2789 :
2790 : ! Zonal wind tendency, broken up into five bins.
2791 0 : real(r8) :: utb(ncol, pver, 5)
2792 : ! Definition of the bin boundaries.
2793 : real(r8), parameter :: bounds(4) = (/ -40._r8, -15._r8, &
2794 : 15._r8, 40._r8 /)
2795 :
2796 : ! Momentum flux in the four cardinal directions.
2797 0 : real(r8) :: mf(ncol, pver, 4)
2798 :
2799 : ! Wave stress in zonal/meridional direction
2800 0 : real(r8) :: taux(ncol,-band%ngwv:band%ngwv,pver)
2801 0 : real(r8) :: tauy(ncol,-band%ngwv:band%ngwv,pver)
2802 :
2803 : ! Temporaries for output
2804 0 : real(r8) :: dummyx(ncol,pver)
2805 0 : real(r8) :: dummyy(ncol,pver)
2806 : ! Variable names
2807 : character(len=10) :: dumc1x, dumc1y
2808 :
2809 :
2810 : ! Accumulate wind tendencies binned according to phase speed.
2811 :
2812 0 : utb = 0._r8
2813 :
2814 : ! Find which output bin the phase speed corresponds to.
2815 0 : ix = find_bin(phase_speeds)
2816 :
2817 : ! Put the wind tendency in that bin.
2818 0 : do l = -band%ngwv, band%ngwv
2819 0 : do k = 1, pver
2820 0 : do i = 1, ncol
2821 0 : utb(i,k,ix(i,l)) = utb(i,k,ix(i,l)) + gwut(i,k,l)
2822 : end do
2823 : end do
2824 : end do
2825 :
2826 : ! Find just the zonal part.
2827 0 : do l = 1, 5
2828 0 : do k = 1, pver
2829 0 : utb(:, k, l) = utb(:, k, l) * xv
2830 : end do
2831 : end do
2832 :
2833 0 : call outfld(trim(prefix)//'UTEND1', utb(:,:,1), ncol, lchnk)
2834 0 : call outfld(trim(prefix)//'UTEND2', utb(:,:,2), ncol, lchnk)
2835 0 : call outfld(trim(prefix)//'UTEND3', utb(:,:,3), ncol, lchnk)
2836 0 : call outfld(trim(prefix)//'UTEND4', utb(:,:,4), ncol, lchnk)
2837 0 : call outfld(trim(prefix)//'UTEND5', utb(:,:,5), ncol, lchnk)
2838 :
2839 : ! Output temperature tendencies due to diffusion and from kinetic energy.
2840 0 : call outfld(trim(prefix)//'TTGWSDF', dttdf / cpair, ncol, lchnk)
2841 0 : call outfld(trim(prefix)//'TTGWSKE', dttke / cpair, ncol, lchnk)
2842 :
2843 :
2844 : ! Output tau broken down into zonal and meridional components.
2845 :
2846 0 : taux = 0._r8
2847 0 : tauy = 0._r8
2848 :
2849 : ! Project phase_speeds, and convert each component to a wavenumber index.
2850 : ! These are mappings from the wavenumber index of tau to those of taux
2851 : ! and tauy, respectively.
2852 0 : do l=-band%ngwv,band%ngwv
2853 0 : ix(:,l) = c_to_l(phase_speeds(:,l)*xv)
2854 0 : iy(:,l) = c_to_l(phase_speeds(:,l)*yv)
2855 : end do
2856 :
2857 : ! Find projection of tau.
2858 0 : do k = 1, pver
2859 0 : do l = -band%ngwv,band%ngwv
2860 0 : do i = 1, ncol
2861 0 : taux(i,ix(i,l),k) = taux(i,ix(i,l),k) &
2862 0 : + abs(tau(i,l,k)*xv(i))
2863 0 : tauy(i,iy(i,l),k) = tauy(i,iy(i,l),k) &
2864 0 : + abs(tau(i,l,k)*yv(i))
2865 : end do
2866 : end do
2867 : end do
2868 :
2869 0 : do l=-band%ngwv,band%ngwv
2870 :
2871 0 : dummyx = taux(:,l,:)
2872 0 : dummyy = tauy(:,l,:)
2873 :
2874 0 : dumc1x = tau_fld_name(l, prefix, x_not_y=.true.)
2875 0 : dumc1y = tau_fld_name(l, prefix, x_not_y=.false.)
2876 :
2877 0 : call outfld(dumc1x,dummyx,ncol,lchnk)
2878 0 : call outfld(dumc1y,dummyy,ncol,lchnk)
2879 :
2880 : enddo
2881 :
2882 :
2883 : ! Output momentum flux in each cardinal direction.
2884 0 : mf = 0._r8
2885 :
2886 0 : do k = 1, pver
2887 :
2888 : ! Convert wind speed components to wavenumber indices.
2889 0 : iu = c_to_l(u(:,k))
2890 0 : iv = c_to_l(v(:,k))
2891 :
2892 : ! Sum tau components in each cardinal direction.
2893 : ! Split west/east and north/south based on whether wave speed exceeds
2894 : ! wind speed.
2895 0 : do l = -band%ngwv, band%ngwv
2896 :
2897 0 : where (iu > l)
2898 0 : mf(:,k,west) = mf(:,k,west) + taux(:,l,k)
2899 : elsewhere
2900 0 : mf(:,k,east) = mf(:,k,east) + taux(:,l,k)
2901 : end where
2902 :
2903 0 : where (iv > l)
2904 0 : mf(:,k,south) = mf(:,k,south) + tauy(:,l,k)
2905 : elsewhere
2906 0 : mf(:,k,north) = mf(:,k,north) + tauy(:,l,k)
2907 : end where
2908 :
2909 : end do
2910 :
2911 : end do
2912 :
2913 0 : call outfld(trim(prefix)//'WMF',mf(:,:,west),ncol,lchnk)
2914 0 : call outfld(trim(prefix)//'EMF',mf(:,:,east),ncol,lchnk)
2915 0 : call outfld(trim(prefix)//'SMF',mf(:,:,south),ncol,lchnk)
2916 0 : call outfld(trim(prefix)//'NMF',mf(:,:,north),ncol,lchnk)
2917 :
2918 : ! Simple output fields written to history file.
2919 : ! Total wind tendencies.
2920 0 : call outfld (trim(prefix)//'UTGWSPEC', utgw , ncol, lchnk)
2921 0 : call outfld (trim(prefix)//'VTGWSPEC', vtgw , ncol, lchnk)
2922 0 : call outfld (trim(prefix)//'TTGWSPEC', ttgw , ncol, lchnk)
2923 :
2924 : ! Tau in each direction.
2925 0 : call outfld (trim(prefix)//'TAUE', taucd(:,:,east), ncol, lchnk)
2926 0 : call outfld (trim(prefix)//'TAUW', taucd(:,:,west), ncol, lchnk)
2927 0 : call outfld (trim(prefix)//'TAUN', taucd(:,:,north), ncol, lchnk)
2928 0 : call outfld (trim(prefix)//'TAUS', taucd(:,:,south), ncol, lchnk)
2929 :
2930 : call outfld (trim(prefix)//'TAUNET', taucd(:,:,east)+taucd(:,:,west), &
2931 0 : ncol, lchnk)
2932 :
2933 : contains
2934 :
2935 : ! Given a value, finds which bin marked by "bounds" the value falls
2936 : ! into.
2937 0 : elemental function find_bin(val) result(idx)
2938 : real(r8), intent(in) :: val
2939 :
2940 : integer :: idx
2941 :
2942 : ! We just have to count how many bounds are exceeded.
2943 0 : if (val >= 0._r8) then
2944 0 : idx = count(val > bounds) + 1
2945 : else
2946 0 : idx = count(val >= bounds) + 1
2947 : end if
2948 :
2949 0 : end function find_bin
2950 :
2951 : ! Convert a speed to a wavenumber between -ngwv and ngwv.
2952 0 : elemental function c_to_l(c) result(l)
2953 : real(r8), intent(in) :: c
2954 :
2955 : integer :: l
2956 :
2957 0 : l = min( max(int(c/band%dc),-band%ngwv), band%ngwv )
2958 :
2959 0 : end function c_to_l
2960 :
2961 : end subroutine gw_spec_outflds
2962 :
2963 : !==========================================================================
2964 :
2965 : ! Generates names for tau output across the wave spectrum (e.g.
2966 : ! BTAUXSn01 or TAUYSp05).
2967 : ! Probably this should use a wavenumber dimension on one field rather
2968 : ! than creating a ton of numbered fields.
2969 0 : character(len=9) pure function tau_fld_name(l, prefix, x_not_y)
2970 : ! Wavenumber
2971 : integer, intent(in) :: l
2972 : ! Single-character prefix for output
2973 : character(len=1), intent(in) :: prefix
2974 : ! X or Y?
2975 : logical, intent(in) :: x_not_y
2976 :
2977 : character(len=2) :: num_str
2978 :
2979 0 : tau_fld_name = trim(prefix)
2980 :
2981 0 : tau_fld_name = trim(tau_fld_name)//"TAU"
2982 :
2983 0 : if (x_not_y) then
2984 0 : tau_fld_name = trim(tau_fld_name)//"XS"
2985 : else
2986 0 : tau_fld_name = trim(tau_fld_name)//"YS"
2987 : end if
2988 :
2989 0 : if (l < 0) then
2990 0 : tau_fld_name = trim(tau_fld_name)//"n"
2991 : else
2992 0 : tau_fld_name = trim(tau_fld_name)//"p"
2993 : end if
2994 :
2995 0 : write(num_str,'(I2.2)') abs(l)
2996 :
2997 0 : tau_fld_name = trim(tau_fld_name)//num_str
2998 :
2999 0 : end function tau_fld_name
3000 :
3001 : !==========================================================================
3002 :
3003 : end module gw_drag
|
