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