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