FV3 Bundle
gw_drag_d.F90
Go to the documentation of this file.
1 ! Generated by TAPENADE (INRIA, Tropics team)
2 ! Tapenade 3.2 (r3024) - 06/17/2009 13:03
3 !
4 ! $Id: gw_drag_d.F90,v 1.2 2013/09/10 19:30:04 jgkim Exp $
5 MODULE gw_drag_d
6  USE mapl_constantsmod, ONLY : mapl_p00, mapl_cp, mapl_grav, &
7 & mapl_rgas, mapl_vireps, mapl_kappa
8  USE gw_drag_init
9 
10 ! USE GW_DRAG_B
11 
12  IMPLICIT NONE
13 !---------------------------------------------------------------------------------
14 ! Purpose:
15 !
16 ! Module to compute the forcing due to parameterized gravity waves. Both an
17 ! orographic and an internal source spectrum are considered.
18 !
19 ! Author: Byron Boville
20 ! In-Sun Song
21 !
22 !---------------------------------------------------------------------------------
23 !save
24 ! Make default type private to the module
25  PRIVATE
26 !
27 ! PUBLIC: interfaces
28 !
29 ! interface to actual parameterization
30 ! PUBLIC gw_intr
31 !, gw_drag_prof, gw_bgnd, gw_oro, gw_prof
32  PUBLIC gw_main_d
33 !, gw_intr_d, gw_drag_prof_d, gw_oro_d
34 ! PUBLIC gw_main_d, gw_intr_d, gw_drag_prof_d, gw_bgnd_d, gw_oro_d
35 !
36 ! PRIVATE: Rest of the data and interfaces are private to this module
37 !
38 ! effective horizontal wave number for background
39  REAL, PARAMETER, PUBLIC :: kwvb=6.28e-5
40 ! effective horizontal wave number for background
41  REAL, PARAMETER, PUBLIC :: kwvbeq=6.28e-5/7.
42 ! effective horizontal wave number for orographic
43  REAL, PARAMETER, PUBLIC :: kwvo=6.28e-5
44 ! fraction of stress deposited in low level region
45  REAL, PARAMETER, PUBLIC :: fracldv=0.0
46 ! max asymmetry between tau(c) and tau(-c)
47  REAL, PARAMETER, PUBLIC :: mxasym=0.1
48 ! max range of tau for all c
49  REAL, PARAMETER, PUBLIC :: mxrange=0.001
50 ! min value of bouyancy frequency
51  REAL, PARAMETER, PUBLIC :: n2min=1.e-8
52 ! critical froude number
53  REAL, PARAMETER, PUBLIC :: fcrit2=0.5
54 ! min surface displacment height for orographic waves
55  REAL, PARAMETER, PUBLIC :: orohmin=10.
56 ! min wind speed for orographic waves
57  REAL, PARAMETER, PUBLIC :: orovmin=2.0
58 ! background source strength (/TAUSCAL)
59  REAL, PARAMETER, PUBLIC :: taubgnd=6.4
60 ! minimum (nonzero) stress
61  REAL, PARAMETER, PUBLIC :: taumin=1.e-10
62 ! scale factor for background stress source
63  REAL, PARAMETER, PUBLIC :: tauscal=0.001
64 ! maximum wind tendency
65  REAL, PARAMETER, PUBLIC :: tndmax=500./86400.
66 ! factor to limit tendency to prevent reversing u-c
67  REAL, PARAMETER, PUBLIC :: umcfac=0.5
68 ! min (u-c)**2
69  REAL, PARAMETER, PUBLIC :: ubmc2mn=0.1
70 ! constant for determining zldv from tau0
71  REAL, PARAMETER, PUBLIC :: zldvcon=10.
72  REAL, PARAMETER, PUBLIC :: rog=mapl_rgas/mapl_grav
73 ! 1/2 * horizontal wavenumber
74  REAL, PARAMETER, PUBLIC :: oroko2=0.5*kwvo
75 ! This is *not* MAPL_PI
76  REAL, PARAMETER, PUBLIC :: pi_gwd=4.0*atan(1.0)
77 
78 CONTAINS
79 SUBROUTINE gw_main_d(pcols, pver, dt, pgwv, effgworo_dev, &
80 & effgwbkg_dev, pint_dev, t_dev, u_dev, u_devd, v_dev, v_devd, &
81 & sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, &
82 & qvt_dev, rog, mapl_vireps_, rlat_dev)
83  IMPLICIT NONE
84 !-----------------------------------------------------------------------
85 ! Interface for multiple gravity wave drag parameterization.
86 !-----------------------------------------------------------------------
87 !------------------------------Arguments--------------------------------
88 ! - number of columns
89  INTEGER :: pcols
90 ! - number of vertical layers
91  INTEGER :: pver
92 ! - time step
93  REAL :: dt
94 ! - number of waves allowed (Default = 4, 0 nullifies)
95  INTEGER :: pgwv
96 ! - tendency efficiency for background gwd (Default = 0.125)
97  REAL :: effgwbkg_dev(pcols)
98 ! - tendency efficiency for orographic gwd (Default = 0.125)
99  REAL :: effgworo_dev(pcols)
100 ! - pressure at the layer edges
101  REAL :: pint_dev(pcols, pver+1)
102 ! - temperature at layers
103  REAL :: t_dev(pcols, pver)
104 ! - zonal wind at layers
105  REAL :: u_dev(pcols, pver)
106  REAL :: u_devd(pcols, pver)
107 ! - meridional wind at layers
108  REAL :: v_dev(pcols, pver)
109  REAL :: v_devd(pcols, pver)
110 ! e standard deviation of orography
111  REAL :: sgh_dev(pcols)
112 ! e reference pressure at the layeredges
113  REAL :: pref_dev(pver+1)
114 ! c pressure at the layers
115  REAL :: pmid_dev(pcols, pver)
116 ! c pressure thickness at the layers
117  REAL :: pdel_dev(pcols, pver)
118 ! c 1.0 / pdel
119  REAL :: rpdel_dev(pcols, pver)
120 ! c log(pint)
121  REAL :: lnpint_dev(pcols, pver+1)
122 ! c height above surface at layers
123  REAL :: zm_dev(pcols, pver)
124 ! c height above surface at layers
125  REAL :: zi_dev(pcols, pver+1)
126 ! c pressure at the layers
127  REAL :: qvt_dev(pcols, pver)
128 !mg latitude in radian
129  REAL :: rlat_dev(pcols)
130 ! zonal wind tendency at layer
131  REAL :: dudt_gwd_dev(pcols, pver)
132  REAL :: dudt_gwd_devd(pcols, pver)
133 ! meridional wind tendency at layer
134  REAL :: dvdt_gwd_dev(pcols, pver)
135  REAL :: dvdt_gwd_devd(pcols, pver)
136 ! temperature tendency at layer
137  REAL :: dtdt_gwd_dev(pcols, pver)
138  REAL :: dtdt_gwd_devd(pcols, pver)
139 ! zonal wind tendency at layer due to orography GWD
140  REAL :: dudt_org_dev(pcols, pver)
141 ! meridional wind tendency at layer due to orography GWD
142  REAL :: dvdt_org_dev(pcols, pver)
143 ! temperature tendency at layer due to orography GWD
144  REAL :: dtdt_org_dev(pcols, pver)
145 ! zonal gravity wave surface stress
146  REAL :: taugwdx_dev(pcols)
147 ! meridional gravity wave surface stress
148  REAL :: taugwdy_dev(pcols)
149 ! zonal orographic gravity wave stress
150  REAL :: tauox_dev(pcols, pver+1)
151 ! meridional orographic gravity wave stress
152  REAL :: tauoy_dev(pcols, pver+1)
153 ! energy flux of orographic gravity waves
154  REAL :: feo_dev(pcols, pver+1)
155 ! pseudoenergy flux of orographic gravity waves
156  REAL :: fepo_dev(pcols, pver+1)
157 ! zonal gravity wave background stress
158  REAL :: taubkgx_dev(pcols)
159 ! meridional gravity wave background stress
160  REAL :: taubkgy_dev(pcols)
161 ! zonal background gravity wave stress
162  REAL :: taubx_dev(pcols, pver+1)
163 ! meridional background gravity wave stress
164  REAL :: tauby_dev(pcols, pver+1)
165 ! energy flux of background gravity waves
166  REAL :: feb_dev(pcols, pver+1)
167 ! pseudoenergy flux of background gravity waves
168  REAL :: fepb_dev(pcols, pver+1)
169 ! dU/dt below background launch level
170  REAL :: utbsrc_dev(pcols, pver)
171 ! dV/dt below background launch level
172  REAL :: vtbsrc_dev(pcols, pver)
173 ! dT/dt below background launch level
174  REAL :: ttbsrc_dev(pcols, pver)
175 !---------------------------Local storage-------------------------------
176 ! loop indexes
177  INTEGER :: i, k, kc
178 ! launch-level index for orographic
179  INTEGER :: kbotoro
180 ! launch-level index for background
181  INTEGER :: kbotbg
182 ! top interface of gwd region
183  INTEGER :: ktopbg, ktoporo
184 ! top interface of low level stress divergence region
185  INTEGER :: kldv
186 ! min value of kldv
187  INTEGER :: kldvmn
188 ! index of top interface of source region
189  INTEGER :: ksrc
190 ! min value of ksrc
191  INTEGER :: ksrcmn
192 ! temperature tendency
193  REAL :: ttgw(pver)
194 ! zonal wind tendency
195  REAL :: utgw(pver)
196 ! meridional wind tendency
197  REAL :: vtgw(pver)
198 ! interface Brunt-Vaisalla frequency
199  REAL :: ni(0:pver)
200 ! midpoint Brunt-Vaisalla frequency
201  REAL :: nm(pver)
202 ! 1/dp across low level divergence region
203  REAL :: rdpldv
204 ! interface density
205  REAL :: rhoi(0:pver)
206 ! c=0 sfc. stress (zonal)
207  REAL :: tau0x
208 ! c=0 sfc. stress (meridional)
209  REAL :: tau0y
210 ! interface temperature
211  REAL :: ti(0:pver)
212 ! projection of wind at interfaces
213  REAL :: ubi(0:pver)
214 ! projection of wind at midpoints
215  REAL :: ubm(pver)
216 ! unit vectors of source wind (x)
217  REAL :: xv
218 ! unit vectors of source wind (y)
219  REAL :: yv
220  REAL :: utosrc(pver)
221  REAL :: vtosrc(pver)
222  REAL :: ttosrc(pver)
223 ! newtonian cooling coefficients
224  REAL :: alpha(0:pver)
225 ! newtonian cooling coefficients
226  REAL :: dback(0:pver)
227 ! wave phase speeds
228  REAL :: c(-pgwv:pgwv)
229 ! - temperature at layers
230  REAL :: att_dev(pcols, pver)
231  REAL :: att_devd(pcols, pver)
232  REAL :: pwr1
233  REAL :: pwx2
234  REAL :: pwy2
235  REAL :: pwr2
236  REAL :: hkl, hkk, tvfac, tv, rog
237  REAL :: pwr10
238  REAL :: pwx20
239  REAL :: pwy20
240  REAL :: pwr20
241  REAL :: mapl_vireps_
242  isize_s1_tau=-pgwv;isize_s2_tau=0;isize_e1_tau=pgwv;isize_e2_tau=pver
243 
244  DO i=1,pcols
245  DO k=1,pver
246  pwr10 = mapl_p00**mapl_kappa
247  pwx20 = mapl_p00/pint_dev(i, k)
248  pwy20 = -mapl_kappa
249  pwr20 = pwx20**pwy20
250  att_dev(i, k) = t_dev(i, k)*pwr10*pwr20
251  END DO
252  END DO
253  DO i=1,pcols
254  zi_dev(i, pver+1) = 0.
255  DO k=pver,1,-1
256  hkl = lnpint_dev(i, k+1) - lnpint_dev(i, k)
257  hkk = 1. - pint_dev(i, k)*hkl*rpdel_dev(i, k)
258  tvfac = 1. + mapl_vireps_*qvt_dev(i, k)
259  tv = att_dev(i, k)*tvfac
260  zm_dev(i, k) = zi_dev(i, k+1) + rog*tv*hkk
261  zi_dev(i, k) = zi_dev(i, k+1) + rog*tv*hkl
262  END DO
263  END DO
264 
265  dvdt_gwd_devd = 0.0
266  dudt_gwd_devd = 0.0
267 
268  DO i=1,pcols
269  If (.not.bypass_bgnd) then
270  do_bgnd=.true.
271  do_oro=.false.
272 
273  CALL gw_intr_d(i, pcols, pver, dt, pgwv, effgworo_dev, effgwbkg_dev&
274 & , pint_dev, att_dev, u_dev, u_devd, v_dev, v_devd&
275 & , sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, &
276 & lnpint_dev, zm_dev, rlat_dev, dudt_gwd_dev, dudt_gwd_devd, &
277 & dvdt_gwd_dev, dvdt_gwd_devd, dtdt_gwd_dev, &
278 & dudt_org_dev, dvdt_org_dev, dtdt_org_dev, taugwdx_dev, &
279 & taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, taubkgx_dev, &
280 & taubkgy_dev, taubx_dev, tauby_dev, feb_dev, fepo_dev, &
281 & fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
282 
283  Endif
284 
285  If (.not.bypass_bgnd.and. .not.bypass_oro) then
286  dvdt_gwd_devd = 0.0
287  dudt_gwd_devd = 0.0
288  Endif
289 
290  If (.not.bypass_oro) then
291  do_bgnd=.false.
292  do_oro=.true.
293 
294  CALL gw_intr_d(i, pcols, pver, dt, pgwv, effgworo_dev, effgwbkg_dev&
295 & , pint_dev, att_dev, u_dev, u_devd, v_dev, v_devd&
296 & , sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, &
297 & lnpint_dev, zm_dev, rlat_dev, dudt_gwd_dev, dudt_gwd_devd, &
298 & dvdt_gwd_dev, dvdt_gwd_devd, dtdt_gwd_dev, &
299 & dudt_org_dev, dvdt_org_dev, dtdt_org_dev, taugwdx_dev, &
300 & taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, taubkgx_dev, &
301 & taubkgy_dev, taubx_dev, tauby_dev, feb_dev, fepo_dev, &
302 & fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
303  Endif
304  END DO
305 
306  DO i=1,pcols
307  DO k=1,pver
308  pwr10 = mapl_p00**mapl_kappa
309  pwx20 = mapl_p00/pint_dev(i, k)
310  pwy20 = -mapl_kappa
311  pwr20 = pwx20**pwy20
312  t_dev(i, k) = att_dev(i, k)/pwr10/pwr20
313  END DO
314  END DO
315 
316 END SUBROUTINE gw_main_d
317  SUBROUTINE gw_intr_d(i, pcols, pver, dt, pgwv, effgworo_dev, &
318 & effgwbkg_dev, pint_dev, t_dev, u_dev, u_devd, v_dev, v_devd, sgh_dev&
319 & , pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, &
320 & rlat_dev, dudt_gwd_dev, dudt_gwd_devd, dvdt_gwd_dev, dvdt_gwd_devd, &
321 & dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev, dtdt_org_dev, taugwdx_dev&
322 & , taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, taubkgx_dev, &
323 & taubkgy_dev, taubx_dev, tauby_dev, feb_dev, fepo_dev, fepb_dev, &
324 & utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
325  IMPLICIT NONE
326 !-----------------------------------------------------------------------
327 ! Interface for multiple gravity wave drag parameterization.
328 !-----------------------------------------------------------------------
329 !------------------------------Arguments--------------------------------
330 ! - number of columns
331  INTEGER, INTENT(IN) :: pcols
332 ! - number of vertical layers
333  INTEGER, INTENT(IN) :: pver
334 ! - time step
335  REAL, INTENT(IN) :: dt
336 ! - number of waves allowed (Default = 4, 0 nullifies)
337  INTEGER, INTENT(IN) :: pgwv
338 ! - tendency efficiency for background gwd (Default = 0.125)
339  REAL :: effgwbkg_dev(pcols)
340 ! - tendency efficiency for orographic gwd (Default = 0.125)
341  REAL :: effgworo_dev(pcols)
342 ! - pressure at the layer edges
343  REAL :: pint_dev(pcols, pver+1)
344 ! - temperature at layers
345  REAL :: t_dev(pcols, pver)
346 ! - zonal wind at layers
347  REAL :: u_dev(pcols, pver)
348  REAL :: u_devd(pcols, pver)
349 ! - meridional wind at layers
350  REAL :: v_dev(pcols, pver)
351  REAL :: v_devd(pcols, pver)
352 ! e standard deviation of orography
353  REAL :: sgh_dev(pcols)
354 ! e reference pressure at the layeredges
355  REAL :: pref_dev(pver+1)
356 ! c pressure at the layers
357  REAL :: pmid_dev(pcols, pver)
358 ! c pressure thickness at the layers
359  REAL :: pdel_dev(pcols, pver)
360 ! c 1.0 / pdel
361  REAL :: rpdel_dev(pcols, pver)
362 ! c log(pint)
363  REAL :: lnpint_dev(pcols, pver+1)
364 ! c height above surface at layers
365  REAL :: zm_dev(pcols, pver)
366 !mg latitude in radian
367  REAL :: rlat_dev(pcols)
368 ! zonal wind tendency at layer
369  REAL :: dudt_gwd_dev(pcols, pver)
370  REAL :: dudt_gwd_devd(pcols, pver)
371 ! meridional wind tendency at layer
372  REAL :: dvdt_gwd_dev(pcols, pver)
373  REAL :: dvdt_gwd_devd(pcols, pver)
374 ! temperature tendency at layer
375  REAL :: dtdt_gwd_dev(pcols, pver)
376 ! zonal wind tendency at layer due to orography GWD
377  REAL :: dudt_org_dev(pcols, pver)
378 ! meridional wind tendency at layer due to orography GWD
379  REAL :: dvdt_org_dev(pcols, pver)
380 ! temperature tendency at layer due to orography GWD
381  REAL :: dtdt_org_dev(pcols, pver)
382 ! zonal gravity wave surface stress
383  REAL :: taugwdx_dev(pcols)
384 ! meridional gravity wave surface stress
385  REAL :: taugwdy_dev(pcols)
386 ! zonal orographic gravity wave stress
387  REAL :: tauox_dev(pcols, pver+1)
388 ! meridional orographic gravity wave stress
389  REAL :: tauoy_dev(pcols, pver+1)
390 ! energy flux of orographic gravity waves
391  REAL :: feo_dev(pcols, pver+1)
392 ! pseudoenergy flux of orographic gravity waves
393  REAL :: fepo_dev(pcols, pver+1)
394 ! zonal gravity wave background stress
395  REAL :: taubkgx_dev(pcols)
396 ! meridional gravity wave background stress
397  REAL :: taubkgy_dev(pcols)
398 ! zonal background gravity wave stress
399  REAL :: taubx_dev(pcols, pver+1)
400 ! meridional background gravity wave stress
401  REAL :: tauby_dev(pcols, pver+1)
402 ! energy flux of background gravity waves
403  REAL :: feb_dev(pcols, pver+1)
404 ! pseudoenergy flux of background gravity waves
405  REAL :: fepb_dev(pcols, pver+1)
406 ! dU/dt below background launch level
407  REAL :: utbsrc_dev(pcols, pver)
408 ! dV/dt below background launch level
409  REAL :: vtbsrc_dev(pcols, pver)
410 ! dT/dt below background launch level
411  REAL :: ttbsrc_dev(pcols, pver)
412 !---------------------------Local storage-------------------------------
413 ! loop indexes
414  INTEGER :: i, ii, k, kc
415 ! launch-level index for orographic
416  INTEGER :: kbotoro
417 ! launch-level index for background
418  INTEGER :: kbotbg
419 ! top interface of gwd region
420  INTEGER :: ktopbg, ktoporo
421 ! top interface of low level stress divergence region
422  INTEGER :: kldv
423 ! min value of kldv
424  INTEGER :: kldvmn
425 ! index of top interface of source region
426  INTEGER :: ksrc
427 ! min value of ksrc
428  INTEGER :: ksrcmn
429 ! temperature tendency
430  REAL :: ttgw(pver)
431 ! zonal wind tendency
432  REAL :: utgw(pver)
433  REAL :: utgwd(pver)
434 ! meridional wind tendency
435  REAL :: vtgw(pver)
436  REAL :: vtgwd(pver)
437 ! interface Brunt-Vaisalla frequency
438  REAL :: ni(0:pver)
439 ! midpoint Brunt-Vaisalla frequency
440  REAL :: nm(pver)
441 ! 1/dp across low level divergence region
442  REAL :: rdpldv
443 ! interface density
444  REAL :: rhoi(0:pver)
445 ! c=0 sfc. stress (zonal)
446  REAL :: tau0x
447 ! c=0 sfc. stress (meridional)
448  REAL :: tau0y
449 ! interface temperature
450  REAL :: ti(0:pver)
451 ! projection of wind at interfaces
452  REAL :: ubi(0:pver)
453  REAL :: ubid(0:pver)
454 ! projection of wind at midpoints
455  REAL :: ubm(pver)
456  REAL :: ubmd(pver)
457 ! unit vectors of source wind (x)
458  REAL :: xv
459  REAL :: xvd
460 ! unit vectors of source wind (y)
461  REAL :: yv
462  REAL :: yvd
463  REAL :: utosrc(pver)
464  REAL :: utosrcd(pver)
465  REAL :: vtosrc(pver)
466  REAL :: vtosrcd(pver)
467  REAL :: ttosrc(pver)
468 ! newtonian cooling coefficients
469  REAL :: alpha(0:pver)
470 ! newtonian cooling coefficients
471  REAL :: dback(0:pver)
472 ! wave phase speeds
473  REAL :: c(-pgwv:pgwv)
474  REAL :: cd(-pgwv:pgwv)
475 ! wave Reynolds stress
476  REAL :: tau(-pgwv:pgwv, 0:pver)
477  REAL :: taud(-pgwv:pgwv, 0:pver)
478 !-----------------------------------------------------------------------------
479 ! Assign newtonian cooling coefficients
480 ! -------------------------------------
481  DO k=0,pver
482  alpha(k) = 0.0
483  dback(k) = 0.0
484  END DO
485 ! Assign wave phase speeds
486 ! ------------------------
487  DO kc=-pgwv,pgwv
488  cd(kc) = 0.0_8
489  c(kc) = 10.0*kc
490  END DO
491 ! Determine the bounds of the background and orographic stress regions
492  ktoporo = 0
493  ktopbg = 0
494  kbotoro = pver
495  DO k=0,pver
496  IF (pref_dev(k+1) .LT. p_src) kbotbg = k
497 ! spectrum source at 400 mb
498  END DO
499  DO k=0,pver
500 ! Profiles of background state variables
501  CALL gw_prof(i, k, pcols, pver, u_dev, v_dev, t_dev, pmid_dev, &
502 & pint_dev, rhoi, ni, ti, nm)
503  END DO
504 !-----------------------------------------------------------------------------
505 ! Non-orographic backgound gravity wave spectrum
506 !-----------------------------------------------------------------------------
507  IF (pgwv .GT. 0 .AND. do_bgnd) THEN
508 ! Determine the wave source for a background spectrum at ~400 mb
509  CALL gw_bgnd_d(i, pcols, pver, c, u_dev, u_devd, v_dev, v_devd, &
510 & t_dev, pmid_dev, pint_dev, pdel_dev, rpdel_dev, &
511 & lnpint_dev, rlat_dev, kldv, kldvmn, ksrc, ksrcmn, rdpldv&
512 & , tau, ubi, ubid, ubm, ubmd, xv, xvd, yv, yvd, pgwv, &
513 & kbotbg)
514 ! Solve for the drag profile
515  CALL gw_drag_prof_bgnd_d(i, pcols, pver, pgwv, pgwv, kbotbg, &
516 & ktopbg, c, u_dev, v_dev, t_dev, pint_dev, &
517 & pdel_dev, rpdel_dev, lnpint_dev, rlat_dev, rhoi&
518 & , ni, ti, nm, dt, alpha, dback, kldv, kldvmn, &
519 & ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, ubm&
520 & , xv, xvd, yv, yvd, utgw, utgwd, vtgw, vtgwd, &
521 & ttgw, taubx_dev, tauby_dev, feb_dev, fepb_dev, &
522 & utosrc, utosrcd, vtosrc, vtosrcd, ttosrc, tau0x&
523 & , tau0y, effgwbkg_dev)
524 ! Add the momentum tendencies to the output tendency arrays
525  DO k=1,pver
526  dudt_gwd_devd(i, k) = utgwd(k) + utosrcd(k)
527  dudt_gwd_dev(i, k) = utgw(k) + utosrc(k)
528  dvdt_gwd_devd(i, k) = vtgwd(k) + vtosrcd(k)
529  dvdt_gwd_dev(i, k) = vtgw(k) + vtosrc(k)
530  END DO
531  ELSE
532 ! zero net tendencies if no spectrum computed
533  DO k=1,pver
534  dudt_gwd_devd(i, k) = 0.0_8
535  dudt_gwd_dev(i, k) = 0.
536  dvdt_gwd_devd(i, k) = 0.0_8
537  dvdt_gwd_dev(i, k) = 0.
538  END DO
539  taud(:, :) = 0.0_8
540  vtgwd(:) = 0.0_8
541  utgwd(:) = 0.0_8
542  ubid(:) = 0.0_8
543  ubmd(:) = 0.0_8
544  END IF
545 ! dtdt_gwd_dev(i,k) = 0.
546 !-----------------------------------------------------------------------------
547 ! Orographic stationary gravity wave
548 !-----------------------------------------------------------------------------
549 ! Determine the orographic wave source
550  IF (do_oro) THEN
551  CALL gw_oro_d(i, pcols, pver, pgwv, u_dev, u_devd, v_dev, v_devd, &
552 & t_dev, sgh_dev, pmid_dev, pint_dev, pdel_dev, zm_dev, nm, &
553 & kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, &
554 & ubm, ubmd, xv, xvd, yv, yvd, kbotoro, rlat_dev)
555 ! Solve for the drag profile
556  CALL gw_drag_prof_d(i, pcols, pver, pgwv, 0, kbotoro, ktoporo, c, &
557 & u_dev, v_dev, t_dev, pint_dev, pdel_dev, rpdel_dev, &
558 & lnpint_dev, rlat_dev, rhoi, ni, ti, nm, dt, alpha, &
559 & dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud&
560 & , ubi, ubid, ubm, xv, xvd, yv, yvd, utgw, utgwd, &
561 & vtgw, vtgwd, ttgw, tauox_dev, tauoy_dev, feo_dev, &
562 & fepo_dev, utosrc, vtosrc, ttosrc, tau0x, tau0y, &
563 & effgworo_dev)
564 ! Add the orographic tendencies to the spectrum tendencies
565 ! Compute the temperature tendency from energy conservation (includes spectrum).
566  DO k=1,pver
567  dudt_gwd_devd(i, k) = dudt_gwd_devd(i, k) + utgwd(k)
568  dudt_gwd_dev(i, k) = dudt_gwd_dev(i, k) + utgw(k)
569  dvdt_gwd_devd(i, k) = dvdt_gwd_devd(i, k) + vtgwd(k)
570  dvdt_gwd_dev(i, k) = dvdt_gwd_dev(i, k) + vtgw(k)
571  END DO
572  END IF
573 ! dtdt_gwd_dev(i,k) = dtdt_gwd_dev(i,k) + ttgw(k)
574  DO k=1,pver
575  u_devd(i, k) = u_devd(i, k) + dt*dudt_gwd_devd(i, k)
576  u_dev(i, k) = u_dev(i, k) + dudt_gwd_dev(i, k)*dt
577  v_devd(i, k) = v_devd(i, k) + dt*dvdt_gwd_devd(i, k)
578  v_dev(i, k) = v_dev(i, k) + dvdt_gwd_dev(i, k)*dt
579  END DO
580 ! t_dev(i,k) = t_dev(i,k) + dtdt_gwd_dev(i,k)*dt
581  RETURN
582  END SUBROUTINE gw_intr_d
583  SUBROUTINE gw_intr(i, pcols, pver, dt, pgwv, effgworo_dev, &
584 & effgwbkg_dev, pint_dev, t_dev, u_dev, v_dev, sgh_dev, pref_dev, &
585 & pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, rlat_dev, &
586 & dudt_gwd_dev, dvdt_gwd_dev, dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev&
587 & , dtdt_org_dev, taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, &
588 & feo_dev, taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, &
589 & fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
590  IMPLICIT NONE
591 !-----------------------------------------------------------------------
592 ! Interface for multiple gravity wave drag parameterization.
593 !-----------------------------------------------------------------------
594 !------------------------------Arguments--------------------------------
595 ! - number of columns
596  INTEGER, INTENT(IN) :: pcols
597 ! - number of vertical layers
598  INTEGER, INTENT(IN) :: pver
599 ! - time step
600  REAL, INTENT(IN) :: dt
601 ! - number of waves allowed (Default = 4, 0 nullifies)
602  INTEGER, INTENT(IN) :: pgwv
603 ! - tendency efficiency for background gwd (Default = 0.125)
604  REAL :: effgwbkg_dev(pcols)
605 ! - tendency efficiency for orographic gwd (Default = 0.125)
606  REAL :: effgworo_dev(pcols)
607 ! - pressure at the layer edges
608  REAL :: pint_dev(pcols, pver+1)
609 ! - temperature at layers
610  REAL :: t_dev(pcols, pver)
611 ! - zonal wind at layers
612  REAL :: u_dev(pcols, pver)
613 ! - meridional wind at layers
614  REAL :: v_dev(pcols, pver)
615 ! e standard deviation of orography
616  REAL :: sgh_dev(pcols)
617 ! e reference pressure at the layeredges
618  REAL :: pref_dev(pver+1)
619 ! c pressure at the layers
620  REAL :: pmid_dev(pcols, pver)
621 ! c pressure thickness at the layers
622  REAL :: pdel_dev(pcols, pver)
623 ! c 1.0 / pdel
624  REAL :: rpdel_dev(pcols, pver)
625 ! c log(pint)
626  REAL :: lnpint_dev(pcols, pver+1)
627 ! c height above surface at layers
628  REAL :: zm_dev(pcols, pver)
629 !mg latitude in radian
630  REAL :: rlat_dev(pcols)
631 ! zonal wind tendency at layer
632  REAL :: dudt_gwd_dev(pcols, pver)
633 ! meridional wind tendency at layer
634  REAL :: dvdt_gwd_dev(pcols, pver)
635 ! temperature tendency at layer
636  REAL :: dtdt_gwd_dev(pcols, pver)
637 ! zonal wind tendency at layer due to orography GWD
638  REAL :: dudt_org_dev(pcols, pver)
639 ! meridional wind tendency at layer due to orography GWD
640  REAL :: dvdt_org_dev(pcols, pver)
641 ! temperature tendency at layer due to orography GWD
642  REAL :: dtdt_org_dev(pcols, pver)
643 ! zonal gravity wave surface stress
644  REAL :: taugwdx_dev(pcols)
645 ! meridional gravity wave surface stress
646  REAL :: taugwdy_dev(pcols)
647 ! zonal orographic gravity wave stress
648  REAL :: tauox_dev(pcols, pver+1)
649 ! meridional orographic gravity wave stress
650  REAL :: tauoy_dev(pcols, pver+1)
651 ! energy flux of orographic gravity waves
652  REAL :: feo_dev(pcols, pver+1)
653 ! pseudoenergy flux of orographic gravity waves
654  REAL :: fepo_dev(pcols, pver+1)
655 ! zonal gravity wave background stress
656  REAL :: taubkgx_dev(pcols)
657 ! meridional gravity wave background stress
658  REAL :: taubkgy_dev(pcols)
659 ! zonal background gravity wave stress
660  REAL :: taubx_dev(pcols, pver+1)
661 ! meridional background gravity wave stress
662  REAL :: tauby_dev(pcols, pver+1)
663 ! energy flux of background gravity waves
664  REAL :: feb_dev(pcols, pver+1)
665 ! pseudoenergy flux of background gravity waves
666  REAL :: fepb_dev(pcols, pver+1)
667 ! dU/dt below background launch level
668  REAL :: utbsrc_dev(pcols, pver)
669 ! dV/dt below background launch level
670  REAL :: vtbsrc_dev(pcols, pver)
671 ! dT/dt below background launch level
672  REAL :: ttbsrc_dev(pcols, pver)
673 !---------------------------Local storage-------------------------------
674 ! loop indexes
675  INTEGER :: i, ii, k, kc
676 ! launch-level index for orographic
677  INTEGER :: kbotoro
678 ! launch-level index for background
679  INTEGER :: kbotbg
680 ! top interface of gwd region
681  INTEGER :: ktopbg, ktoporo
682 ! top interface of low level stress divergence region
683  INTEGER :: kldv
684 ! min value of kldv
685  INTEGER :: kldvmn
686 ! index of top interface of source region
687  INTEGER :: ksrc
688 ! min value of ksrc
689  INTEGER :: ksrcmn
690 ! temperature tendency
691  REAL :: ttgw(pver)
692 ! zonal wind tendency
693  REAL :: utgw(pver)
694 ! meridional wind tendency
695  REAL :: vtgw(pver)
696 ! interface Brunt-Vaisalla frequency
697  REAL :: ni(0:pver)
698 ! midpoint Brunt-Vaisalla frequency
699  REAL :: nm(pver)
700 ! 1/dp across low level divergence region
701  REAL :: rdpldv
702 ! interface density
703  REAL :: rhoi(0:pver)
704 ! c=0 sfc. stress (zonal)
705  REAL :: tau0x
706 ! c=0 sfc. stress (meridional)
707  REAL :: tau0y
708 ! interface temperature
709  REAL :: ti(0:pver)
710 ! projection of wind at interfaces
711  REAL :: ubi(0:pver)
712 ! projection of wind at midpoints
713  REAL :: ubm(pver)
714 ! unit vectors of source wind (x)
715  REAL :: xv
716 ! unit vectors of source wind (y)
717  REAL :: yv
718  REAL :: utosrc(pver)
719  REAL :: vtosrc(pver)
720  REAL :: ttosrc(pver)
721 ! newtonian cooling coefficients
722  REAL :: alpha(0:pver)
723 ! newtonian cooling coefficients
724  REAL :: dback(0:pver)
725 ! wave phase speeds
726  REAL :: c(-pgwv:pgwv)
727 ! wave Reynolds stress
728  REAL :: tau(-pgwv:pgwv, 0:pver)
729 !-----------------------------------------------------------------------------
730 ! Assign newtonian cooling coefficients
731 ! -------------------------------------
732  DO k=0,pver
733  alpha(k) = 0.0
734  dback(k) = 0.0
735  END DO
736 ! Assign wave phase speeds
737 ! ------------------------
738  DO kc=-pgwv,pgwv
739  c(kc) = 10.0*kc
740  END DO
741 ! Determine the bounds of the background and orographic stress regions
742  ktoporo = 0
743  ktopbg = 0
744  kbotoro = pver
745  DO k=0,pver
746  IF (pref_dev(k+1) .LT. p_src) kbotbg = k
747 ! spectrum source at 400 mb
748  END DO
749  DO k=0,pver
750 ! Profiles of background state variables
751  CALL gw_prof(i, k, pcols, pver, u_dev, v_dev, t_dev, pmid_dev, &
752 & pint_dev, rhoi, ni, ti, nm)
753  END DO
754 !-----------------------------------------------------------------------------
755 ! Non-orographic backgound gravity wave spectrum
756 !-----------------------------------------------------------------------------
757  IF (pgwv .GT. 0 .AND. do_bgnd) THEN
758 ! Determine the wave source for a background spectrum at ~400 mb
759  CALL gw_bgnd(i, pcols, pver, c, u_dev, v_dev, t_dev, pmid_dev, &
760 & pint_dev, pdel_dev, rpdel_dev, lnpint_dev, rlat_dev, kldv, &
761 & kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, pgwv, &
762 & kbotbg)
763 ! Solve for the drag profile
764  CALL gw_drag_prof_bgnd(i, pcols, pver, pgwv, pgwv, kbotbg, ktopbg&
765 & , c, u_dev, v_dev, t_dev, pint_dev, pdel_dev, &
766 & rpdel_dev, lnpint_dev, rlat_dev, rhoi, ni, ti, nm&
767 & , dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, &
768 & rdpldv, tau, ubi, ubm, xv, yv, utgw, vtgw, ttgw, &
769 & taubx_dev, tauby_dev, feb_dev, fepb_dev, utosrc, &
770 & vtosrc, ttosrc, tau0x, tau0y, effgwbkg_dev)
771 ! Add the momentum tendencies to the output tendency arrays
772  DO k=1,pver
773  dudt_gwd_dev(i, k) = utgw(k) + utosrc(k)
774  dvdt_gwd_dev(i, k) = vtgw(k) + vtosrc(k)
775  END DO
776  ELSE
777 ! dudt_gwd_dev(i,k) = utgw(k) + utosrc(k)
778 ! dvdt_gwd_dev(i,k) = vtgw(k) + vtosrc(k)
779 ! dtdt_gwd_dev(i,k) = ttgw(k) + ttosrc(k)
780 ! zero net tendencies if no spectrum computed
781  DO k=1,pver
782  dudt_gwd_dev(i, k) = 0.
783  dvdt_gwd_dev(i, k) = 0.
784  END DO
785  END IF
786 ! dtdt_gwd_dev(i,k) = 0.
787 !-----------------------------------------------------------------------------
788 ! Orographic stationary gravity wave
789 !-----------------------------------------------------------------------------
790 ! Determine the orographic wave source
791  IF (do_oro) THEN
792  CALL gw_oro(i, pcols, pver, pgwv, u_dev, v_dev, t_dev, sgh_dev, &
793 & pmid_dev, pint_dev, pdel_dev, zm_dev, nm, kldv, kldvmn, ksrc&
794 & , ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, kbotoro, rlat_dev)
795 ! Solve for the drag profile
796  CALL gw_drag_prof(i, pcols, pver, pgwv, 0, kbotoro, ktoporo, c, &
797 & u_dev, v_dev, t_dev, pint_dev, pdel_dev, rpdel_dev, &
798 & lnpint_dev, rlat_dev, rhoi, ni, ti, nm, dt, alpha, &
799 & dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, &
800 & ubm, xv, yv, utgw, vtgw, ttgw, tauox_dev, tauoy_dev, &
801 & feo_dev, fepo_dev, utosrc, vtosrc, ttosrc, tau0x, &
802 & tau0y, effgworo_dev)
803 ! Add the orographic tendencies to the spectrum tendencies
804 ! Compute the temperature tendency from energy conservation (includes spectrum).
805  DO k=1,pver
806  dudt_gwd_dev(i, k) = dudt_gwd_dev(i, k) + utgw(k)
807  dvdt_gwd_dev(i, k) = dvdt_gwd_dev(i, k) + vtgw(k)
808  END DO
809  END IF
810 ! dtdt_gwd_dev(i,k) = dtdt_gwd_dev(i,k) + ttgw(k)
811  DO k=1,pver
812  u_dev(i, k) = u_dev(i, k) + dudt_gwd_dev(i, k)*dt
813  v_dev(i, k) = v_dev(i, k) + dvdt_gwd_dev(i, k)*dt
814  END DO
815 ! t_dev(i,k) = t_dev(i,k) + dtdt_gwd_dev(i,k)*dt
816  RETURN
817  END SUBROUTINE gw_intr
818 !================================================================================
819  SUBROUTINE gw_prof(i, k, pcols, pver, u, v, t, pm, pi, rhoi, ni, ti, &
820 & nm)
821  IMPLICIT NONE
822 !-----------------------------------------------------------------------
823 ! Compute profiles of background state quantities for the multiple
824 ! gravity wave drag parameterization.
825 !
826 ! The parameterization is assumed to operate only where water vapor
827 ! concentrations are negligible in determining the density.
828 !-----------------------------------------------------------------------
829 !------------------------------Arguments--------------------------------
830 ! current atmospheric column
831  INTEGER, INTENT(IN) :: i
832 ! current atmospheric layer
833  INTEGER, INTENT(IN) :: k
834 ! number of atmospheric columns
835  INTEGER, INTENT(IN) :: pcols
836 ! number of vertical layers
837  INTEGER, INTENT(IN) :: pver
838 ! midpoint zonal wind
839  REAL :: u(pcols, pver)
840 ! midpoint meridional wind
841  REAL :: v(pcols, pver)
842 ! midpoint temperatures
843  REAL :: t(pcols, pver)
844 ! midpoint pressures
845  REAL :: pm(pcols, pver)
846 ! interface pressures
847  REAL :: pi(pcols, 0:pver)
848 ! interface density
849  REAL :: rhoi(0:pver)
850 ! interface Brunt-Vaisalla frequency
851  REAL :: ni(0:pver)
852 ! interface temperature
853  REAL :: ti(0:pver)
854 ! midpoint Brunt-Vaisalla frequency
855  REAL :: nm(pver)
856 !---------------------------Local storage-------------------------------
857  REAL :: dtdp
858 ! Brunt-Vaisalla frequency squared
859  REAL :: n2
860 ! Brunt-Vaisalla frequency squared
861  REAL :: t_new, t_new_
862  REAL :: arg1
863  INTRINSIC max
864  INTRINSIC sqrt
865  REAL :: max1
866 !-----------------------------------------------------------------------------
867 ! Determine the interface densities and Brunt-Vaisala frequencies.
868 !-----------------------------------------------------------------------------
869 ! The top interface values are calculated assuming an isothermal atmosphere
870 ! above the top level.
871  IF (k .EQ. 0) THEN
872  ti(k) = t(i, k+1)
873  CALL get_ti(t_new, ti(k))
874  rhoi(k) = pi(i, k)/(mapl_rgas*t_new)
875  arg1 = mapl_grav*mapl_grav/(mapl_cp*t_new)
876  ni(k) = sqrt(arg1)
877 ! Interior points use centered differences
878  ELSE IF (k .GT. 0 .AND. k .LT. pver) THEN
879  ti(k) = 0.5*(t(i, k)+t(i, k+1))
880  CALL get_ti(t_new, ti(k))
881  rhoi(k) = pi(i, k)/(mapl_rgas*t_new)
882  dtdp = (t(i, k+1)-t(i, k))/(pm(i, k+1)-pm(i, k))
883  CALL get_ti(t_new_, dtdp)
884  n2 = mapl_grav*mapl_grav/t_new*(1./mapl_cp-rhoi(k)*t_new_)
885  IF (n2min .LT. n2) THEN
886  max1 = n2
887  ELSE
888  max1 = n2min
889  END IF
890  ni(k) = sqrt(max1)
891 ! Bottom interface uses bottom level temperature, density; next interface
892 ! B-V frequency.
893  ELSE IF (k .EQ. pver) THEN
894  ti(k) = t(i, k)
895  CALL get_ti(t_new, ti(k))
896  rhoi(k) = pi(i, k)/(mapl_rgas*t_new)
897  ni(k) = ni(k-1)
898  END IF
899 !-----------------------------------------------------------------------------
900 ! Determine the midpoint Brunt-Vaisala frequencies.
901 !-----------------------------------------------------------------------------
902  IF (k .GT. 0) nm(k) = 0.5*(ni(k-1)+ni(k))
903  RETURN
904  END SUBROUTINE gw_prof
905 ! Differentiation of gw_oro in forward (tangent) mode:
906 ! variations of output variables: tau xv yv ubi
907 ! with respect to input variables: tau u v ubi ubm
908 !================================================================================
909  SUBROUTINE gw_oro_d(i, pcols, pver, pgwv, u, ud, v, vd, t, sgh, pm, pi&
910 & , dpm, zm, nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, &
911 & ubid, ubm, ubmd, xv, xvd, yv, yvd, kbot, rlat)
912  IMPLICIT NONE
913 !-----------------------------------------------------------------------
914 ! Orographic source for multiple gravity wave drag parameterization.
915 !
916 ! The stress is returned for a single wave with c=0, over orography.
917 ! For points where the orographic variance is small (including ocean),
918 ! the returned stress is zero.
919 !------------------------------Arguments--------------------------------
920 ! number of current column
921  INTEGER, INTENT(IN) :: i
922 ! number of atmospheric columns
923  INTEGER, INTENT(IN) :: pcols
924 ! number of atmospheric columns
925  INTEGER, INTENT(IN) :: pver
926 ! number of waves allowed
927  INTEGER, INTENT(IN) :: pgwv
928 ! midpoint zonal wind
929  REAL :: u(pcols, pver)
930  REAL :: ud(pcols, pver)
931 ! midpoint meridional wind
932  REAL :: v(pcols, pver)
933  REAL :: vd(pcols, pver)
934 ! midpoint temperatures
935  REAL :: t(pcols, pver)
936 ! standard deviation of orography
937  REAL :: sgh(pcols)
938 ! midpoint pressures
939  REAL :: pm(pcols, pver)
940 ! interface pressures
941  REAL :: pi(pcols, 0:pver)
942 ! midpoint delta p (pi(k)-pi(k-1))
943  REAL :: dpm(pcols, pver)
944 ! midpoint heights
945  REAL :: zm(pcols, pver)
946 ! midpoint Brunt-Vaisalla frequency
947  REAL :: nm(pver)
948 ! top interface of low level stress div region
949  INTEGER :: kldv
950 ! min value of kldv
951  INTEGER :: kldvmn
952 ! index of top interface of source region
953  INTEGER :: ksrc
954 ! min value of ksrc
955  INTEGER :: ksrcmn
956 ! 1/dp across low level divergence region
957  REAL :: rdpldv
958 ! wave Reynolds stress
959  REAL :: tau(-pgwv:pgwv, 0:pver)
960  REAL :: taud(-pgwv:pgwv, 0:pver)
961 ! projection of wind at interfaces
962  REAL :: ubi(0:pver)
963  REAL :: ubid(0:pver)
964 ! projection of wind at midpoints
965  REAL :: ubm(pver)
966  REAL :: ubmd(pver)
967 ! unit vectors of source wind (x)
968  REAL :: xv
969  REAL :: xvd
970 ! unit vectors of source wind (y)
971  REAL :: yv
972  REAL :: yvd
973  INTEGER :: kbot
974  REAL :: rlat(pcols)
975 !---------------------------Local storage-------------------------------
976 ! loop indexes
977  INTEGER :: k
978 ! Source-layer basic-state wind
979  REAL :: ubsrc
980  REAL :: ubsrcd
981 ! surface streamline displacment height (2*sgh)
982  REAL :: hdsp
983 ! max orographic sdv to use
984  REAL :: sghmax
985  REAL :: sghmaxd
986 ! c=0 stress from orography
987  REAL :: tauoro
988  REAL :: tauorod
989 ! real :: zldv ! top of the low level stress divergence region
990 ! b-f frequency averaged over source region
991  REAL :: nsrc
992 ! interface pressure at top of source region
993  REAL :: psrc
994 ! density averaged over source region
995  REAL :: rsrc
996 ! u wind averaged over source region
997  REAL :: usrc
998  REAL :: usrcd
999 ! v wind averaged over source region
1000  REAL :: vsrc
1001  REAL :: vsrcd
1002  REAL :: u_new, v_new, t_new
1003  REAL :: u_newd, v_newd
1004  REAL :: arg1
1005  REAL :: arg1d
1006  REAL :: result1
1007  REAL :: min1
1008  REAL :: min1d
1009  INTRINSIC min
1010  INTRINSIC sqrt
1011 ! Begins
1012 !---------------------------------------------------------------------------
1013 ! Average the basic state variables for the wave source over the depth of
1014 ! the orographic standard deviation. Here we assume that the apropiate
1015 ! values of wind, stability, etc. for determining the wave source are
1016 ! averages over the depth of the atmosphere pentrated by the typical mountain.
1017 ! Reduces to the bottom midpoint values when sgh=0, such as over ocean.
1018 !
1019 ! Also determine the depth of the low level stress divergence region, as
1020 ! the max of the boundary layer depth and the source region depth. This
1021 ! can be done here if the stress magnitude does not determine the depth,
1022 ! otherwise it must be done below.
1023 !---------------------------------------------------------------------------
1024  ksrc = pver - 1
1025  kldv = pver - 1
1026  psrc = pi(i, pver-1)
1027  CALL get_ti(t_new, t(i, pver))
1028  rsrc = pm(i, pver)/(mapl_rgas*t_new)*dpm(i, pver)
1029  CALL get_uv_d(u_new, u_newd, u(i, pver), ud(i, pver))
1030  CALL get_uv_d(v_new, v_newd, v(i, pver), vd(i, pver))
1031  usrcd = dpm(i, pver)*u_newd
1032  usrc = u_new*dpm(i, pver)
1033  vsrcd = dpm(i, pver)*v_newd
1034  vsrc = v_new*dpm(i, pver)
1035  nsrc = nm(pver)*dpm(i, pver)
1036  hdsp = 2.0*sgh(i)
1037  DO k=pver-1,pver/2,-1
1038  arg1 = zm(i, k)*zm(i, k+1)
1039  result1 = sqrt(arg1)
1040  IF (hdsp .GT. result1) THEN
1041  ksrc = k - 1
1042  kldv = k - 1
1043  psrc = pi(i, k-1)
1044  CALL get_ti(t_new, t(i, k))
1045  rsrc = rsrc + pm(i, k)/(mapl_rgas*t_new)*dpm(i, k)
1046  CALL get_uv_d(u_new, u_newd, u(i, k), ud(i, k))
1047  CALL get_uv_d(v_new, v_newd, v(i, k), vd(i, k))
1048  usrcd = usrcd + dpm(i, k)*u_newd
1049  usrc = usrc + u_new*dpm(i, k)
1050  vsrcd = vsrcd + dpm(i, k)*v_newd
1051  vsrc = vsrc + v_new*dpm(i, k)
1052  nsrc = nsrc + nm(k)*dpm(i, k)
1053  END IF
1054  END DO
1055  rsrc = rsrc/(pi(i, pver)-psrc)
1056  usrcd = usrcd/(pi(i, pver)-psrc)
1057  usrc = usrc/(pi(i, pver)-psrc)
1058  vsrcd = vsrcd/(pi(i, pver)-psrc)
1059  vsrc = vsrc/(pi(i, pver)-psrc)
1060  nsrc = nsrc/(pi(i, pver)-psrc)
1061  IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
1062  ubsrc = sqrt(ubmc2mn)
1063  xv = 1.
1064  yv = 0.
1065  xvd = 0.0_8
1066  yvd = 0.0_8
1067  ubsrcd = 0.0_8
1068  ELSE
1069  arg1d = 2*usrc*usrcd + 2*vsrc*vsrcd
1070  arg1 = usrc**2 + vsrc**2
1071  IF (arg1 .EQ. 0.0) THEN
1072  ubsrcd = 0.0_8
1073  ELSE
1074  ubsrcd = arg1d/(2.0*sqrt(arg1))
1075  END IF
1076  ubsrc = sqrt(arg1)
1077  xvd = (usrcd*ubsrc-usrc*ubsrcd)/ubsrc**2
1078  xv = usrc/ubsrc
1079  yvd = (vsrcd*ubsrc-vsrc*ubsrcd)/ubsrc**2
1080  yv = vsrc/ubsrc
1081  END IF
1082 ! Project the local wind at midpoints onto the source wind.
1083  DO k=1,pver
1084  ubmd(k) = ud(i, k)*xv + u(i, k)*xvd + vd(i, k)*yv + v(i, k)*yvd
1085  ubm(k) = u(i, k)*xv + v(i, k)*yv
1086  END DO
1087 ! Compute the interface wind projection by averaging the midpoint winds.
1088 ! Use the top level wind at the top interface.
1089  ubid(0) = ubmd(1)
1090  ubi(0) = ubm(1)
1091  DO k=1,pver
1092  ubid(k) = ubmd(k)
1093  ubi(k) = ubm(k)
1094  END DO
1095 ! Determine the orographic c=0 source term following McFarlane (1987).
1096 ! Set the source top interface index to pver, if the orographic term is zero.
1097  IF (ubsrc .GT. orovmin .AND. hdsp .GT. orohmin) THEN
1098  sghmaxd = fcrit2*2*ubsrc*ubsrcd/nsrc**2
1099  sghmax = fcrit2*(ubsrc/nsrc)**2
1100  IF (hdsp**2 .GT. sghmax) THEN
1101  min1d = sghmaxd
1102  min1 = sghmax
1103  ELSE
1104  min1 = hdsp**2
1105  min1d = 0.0_8
1106  END IF
1107  tauorod = oroko2*rsrc*nsrc*(min1d*ubsrc+min1*ubsrcd)
1108  tauoro = oroko2*min1*rsrc*nsrc*ubsrc
1109  ELSE
1110  tauoro = 0.
1111  ksrc = pver
1112  kldv = pver
1113  tauorod = 0.0_8
1114  END IF
1115 ! tauoro is nontrivial when ubsrc is positive. However, if ubi(ksrc) is negative
1116 ! [note that the sign of ubsrc is irrelevant to the sign of ubi(ksrc)], orographic
1117 ! GWs can propagative upward passing through the negative basic-state wind.
1118 ! The following is to prevent this physically unjustified simulation.
1119 ! In addition, even if ubi(ksrc) > 0, if ubm(ksrc) < 0 .and. ubi(ksrc-1) < 0,
1120 ! negative wave stress leads to the acceleration of the negative ubm(ksrc).
1121 ! This result is physically inconsistent. However, if ubm(ksrc) > 0 .and.
1122 ! ubi(ksrc-1) < 0, GWs are filtered in physically consistent way, and
1123 ! decelerate the positive ubm(ksrc). Therefore, GWs are also assumed not to be
1124 ! launched when ubm(ksrc) < 0.
1125  IF (ubi(kbot) .LT. 0. .OR. ubm(kbot) .LT. 0.) THEN
1126  tauoro = 0.
1127  ksrc = pver
1128  kldv = pver
1129  tauorod = 0.0_8
1130  END IF
1131 ! Sets kbot equal to ksrc
1132  kbot = ksrc
1133 ! Set the phase speeds and wave numbers in the direction of the source wind.
1134 ! Set the source stress magnitude (positive only, note that the sign of the
1135 ! stress is the same as (c-u).
1136  taud(0, kbot) = tauorod
1137  tau(0, kbot) = tauoro
1138 ! Determine the min value of kldv and ksrc for limiting later loops
1139 ! and the pressure at the top interface of the low level stress divergence
1140 ! region.
1141  ksrcmn = pver
1142  kldvmn = pver
1143  IF (ksrcmn .GT. ksrc) THEN
1144  ksrcmn = ksrc
1145  ELSE
1146  ksrcmn = ksrcmn
1147  END IF
1148  IF (kldvmn .GT. kldv) THEN
1149  kldvmn = kldv
1150  ELSE
1151  kldvmn = kldvmn
1152  END IF
1153  IF (kldv .NE. pver) rdpldv = 1./(pi(i, kldv)-pi(i, pver))
1154 ! kldvmn is always pver because FRACLDV == 0.
1155  IF (fracldv .LE. 0.) kldvmn = pver
1156  RETURN
1157  END SUBROUTINE gw_oro_d
1158 !================================================================================
1159  SUBROUTINE gw_oro(i, pcols, pver, pgwv, u, v, t, sgh, pm, pi, dpm, zm&
1160 & , nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
1161 & kbot, rlat)
1162  IMPLICIT NONE
1163 !-----------------------------------------------------------------------
1164 ! Orographic source for multiple gravity wave drag parameterization.
1165 !
1166 ! The stress is returned for a single wave with c=0, over orography.
1167 ! For points where the orographic variance is small (including ocean),
1168 ! the returned stress is zero.
1169 !------------------------------Arguments--------------------------------
1170 ! number of current column
1171  INTEGER, INTENT(IN) :: i
1172 ! number of atmospheric columns
1173  INTEGER, INTENT(IN) :: pcols
1174 ! number of atmospheric columns
1175  INTEGER, INTENT(IN) :: pver
1176 ! number of waves allowed
1177  INTEGER, INTENT(IN) :: pgwv
1178 ! midpoint zonal wind
1179  REAL :: u(pcols, pver)
1180 ! midpoint meridional wind
1181  REAL :: v(pcols, pver)
1182 ! midpoint temperatures
1183  REAL :: t(pcols, pver)
1184 ! standard deviation of orography
1185  REAL :: sgh(pcols)
1186 ! midpoint pressures
1187  REAL :: pm(pcols, pver)
1188 ! interface pressures
1189  REAL :: pi(pcols, 0:pver)
1190 ! midpoint delta p (pi(k)-pi(k-1))
1191  REAL :: dpm(pcols, pver)
1192 ! midpoint heights
1193  REAL :: zm(pcols, pver)
1194 ! midpoint Brunt-Vaisalla frequency
1195  REAL :: nm(pver)
1196 ! top interface of low level stress div region
1197  INTEGER :: kldv
1198 ! min value of kldv
1199  INTEGER :: kldvmn
1200 ! index of top interface of source region
1201  INTEGER :: ksrc
1202 ! min value of ksrc
1203  INTEGER :: ksrcmn
1204 ! 1/dp across low level divergence region
1205  REAL :: rdpldv
1206 ! wave Reynolds stress
1207  REAL :: tau(-pgwv:pgwv, 0:pver)
1208 ! projection of wind at interfaces
1209  REAL :: ubi(0:pver)
1210 ! projection of wind at midpoints
1211  REAL :: ubm(pver)
1212 ! unit vectors of source wind (x)
1213  REAL :: xv
1214 ! unit vectors of source wind (y)
1215  REAL :: yv
1216  INTEGER :: kbot
1217  REAL :: rlat(pcols)
1218 !---------------------------Local storage-------------------------------
1219 ! loop indexes
1220  INTEGER :: k
1221 ! Source-layer basic-state wind
1222  REAL :: ubsrc
1223 ! surface streamline displacment height (2*sgh)
1224  REAL :: hdsp
1225 ! max orographic sdv to use
1226  REAL :: sghmax
1227 ! c=0 stress from orography
1228  REAL :: tauoro
1229 ! real :: zldv ! top of the low level stress divergence region
1230 ! b-f frequency averaged over source region
1231  REAL :: nsrc
1232 ! interface pressure at top of source region
1233  REAL :: psrc
1234 ! density averaged over source region
1235  REAL :: rsrc
1236 ! u wind averaged over source region
1237  REAL :: usrc
1238 ! v wind averaged over source region
1239  REAL :: vsrc
1240  REAL :: u_new, v_new, t_new
1241  REAL :: arg1
1242  REAL :: result1
1243  REAL :: min1
1244  INTRINSIC min
1245  INTRINSIC sqrt
1246 ! Begins
1247 !---------------------------------------------------------------------------
1248 ! Average the basic state variables for the wave source over the depth of
1249 ! the orographic standard deviation. Here we assume that the apropiate
1250 ! values of wind, stability, etc. for determining the wave source are
1251 ! averages over the depth of the atmosphere pentrated by the typical mountain.
1252 ! Reduces to the bottom midpoint values when sgh=0, such as over ocean.
1253 !
1254 ! Also determine the depth of the low level stress divergence region, as
1255 ! the max of the boundary layer depth and the source region depth. This
1256 ! can be done here if the stress magnitude does not determine the depth,
1257 ! otherwise it must be done below.
1258 !---------------------------------------------------------------------------
1259  ksrc = pver - 1
1260  kldv = pver - 1
1261  psrc = pi(i, pver-1)
1262  CALL get_ti(t_new, t(i, pver))
1263  rsrc = pm(i, pver)/(mapl_rgas*t_new)*dpm(i, pver)
1264  CALL get_uv(u_new, u(i, pver))
1265  CALL get_uv(v_new, v(i, pver))
1266  usrc = u_new*dpm(i, pver)
1267  vsrc = v_new*dpm(i, pver)
1268  nsrc = nm(pver)*dpm(i, pver)
1269  hdsp = 2.0*sgh(i)
1270  DO k=pver-1,pver/2,-1
1271  arg1 = zm(i, k)*zm(i, k+1)
1272  result1 = sqrt(arg1)
1273  IF (hdsp .GT. result1) THEN
1274  ksrc = k - 1
1275  kldv = k - 1
1276  psrc = pi(i, k-1)
1277  CALL get_ti(t_new, t(i, k))
1278  rsrc = rsrc + pm(i, k)/(mapl_rgas*t_new)*dpm(i, k)
1279  CALL get_uv(u_new, u(i, k))
1280  CALL get_uv(v_new, v(i, k))
1281  usrc = usrc + u_new*dpm(i, k)
1282  vsrc = vsrc + v_new*dpm(i, k)
1283  nsrc = nsrc + nm(k)*dpm(i, k)
1284  END IF
1285  END DO
1286  rsrc = rsrc/(pi(i, pver)-psrc)
1287  usrc = usrc/(pi(i, pver)-psrc)
1288  vsrc = vsrc/(pi(i, pver)-psrc)
1289  nsrc = nsrc/(pi(i, pver)-psrc)
1290  IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
1291  ubsrc = sqrt(ubmc2mn)
1292  xv = 1.
1293  yv = 0.
1294  ELSE
1295  arg1 = usrc**2 + vsrc**2
1296  ubsrc = sqrt(arg1)
1297  xv = usrc/ubsrc
1298  yv = vsrc/ubsrc
1299  END IF
1300 ! Project the local wind at midpoints onto the source wind.
1301  DO k=1,pver
1302  ubm(k) = u(i, k)*xv + v(i, k)*yv
1303  END DO
1304 ! Compute the interface wind projection by averaging the midpoint winds.
1305 ! Use the top level wind at the top interface.
1306  ubi(0) = ubm(1)
1307  DO k=1,pver
1308  ubi(k) = ubm(k)
1309  END DO
1310 ! Determine the orographic c=0 source term following McFarlane (1987).
1311 ! Set the source top interface index to pver, if the orographic term is zero.
1312  IF (ubsrc .GT. orovmin .AND. hdsp .GT. orohmin) THEN
1313  sghmax = fcrit2*(ubsrc/nsrc)**2
1314  IF (hdsp**2 .GT. sghmax) THEN
1315  min1 = sghmax
1316  ELSE
1317  min1 = hdsp**2
1318  END IF
1319  tauoro = oroko2*min1*rsrc*nsrc*ubsrc
1320  ELSE
1321  tauoro = 0.
1322  ksrc = pver
1323  kldv = pver
1324  END IF
1325 ! tauoro is nontrivial when ubsrc is positive. However, if ubi(ksrc) is negative
1326 ! [note that the sign of ubsrc is irrelevant to the sign of ubi(ksrc)], orographic
1327 ! GWs can propagative upward passing through the negative basic-state wind.
1328 ! The following is to prevent this physically unjustified simulation.
1329 ! In addition, even if ubi(ksrc) > 0, if ubm(ksrc) < 0 .and. ubi(ksrc-1) < 0,
1330 ! negative wave stress leads to the acceleration of the negative ubm(ksrc).
1331 ! This result is physically inconsistent. However, if ubm(ksrc) > 0 .and.
1332 ! ubi(ksrc-1) < 0, GWs are filtered in physically consistent way, and
1333 ! decelerate the positive ubm(ksrc). Therefore, GWs are also assumed not to be
1334 ! launched when ubm(ksrc) < 0.
1335  IF (ubi(kbot) .LT. 0. .OR. ubm(kbot) .LT. 0.) THEN
1336  tauoro = 0.
1337  ksrc = pver
1338  kldv = pver
1339  END IF
1340 ! Sets kbot equal to ksrc
1341  kbot = ksrc
1342 ! Set the phase speeds and wave numbers in the direction of the source wind.
1343 ! Set the source stress magnitude (positive only, note that the sign of the
1344 ! stress is the same as (c-u).
1345  tau(0, kbot) = tauoro
1346 ! Determine the min value of kldv and ksrc for limiting later loops
1347 ! and the pressure at the top interface of the low level stress divergence
1348 ! region.
1349  ksrcmn = pver
1350  kldvmn = pver
1351  IF (ksrcmn .GT. ksrc) THEN
1352  ksrcmn = ksrc
1353  ELSE
1354  ksrcmn = ksrcmn
1355  END IF
1356  IF (kldvmn .GT. kldv) THEN
1357  kldvmn = kldv
1358  ELSE
1359  kldvmn = kldvmn
1360  END IF
1361  IF (kldv .NE. pver) rdpldv = 1./(pi(i, kldv)-pi(i, pver))
1362 ! kldvmn is always pver because FRACLDV == 0.
1363  IF (fracldv .LE. 0.) kldvmn = pver
1364  RETURN
1365  END SUBROUTINE gw_oro
1366 ! Differentiation of gw_bgnd in forward (tangent) mode:
1367 ! variations of output variables: xv yv ubi ubm
1368 ! with respect to input variables: u v
1369 !===============================================================================
1370  SUBROUTINE gw_bgnd_d(i, pcols, pver, c, u, ud, v, vd, t, pm, pi, dpm, &
1371 & rdpm, piln, rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubid&
1372 & , ubm, ubmd, xv, xvd, yv, yvd, ngwv, kbot)
1373  IMPLICIT NONE
1374 !-----------------------------------------------------------------------
1375 ! Driver for multiple gravity wave drag parameterization.
1376 !
1377 ! The parameterization is assumed to operate only where water vapor
1378 ! concentrations are negligible in determining the density.
1379 !-----------------------------------------------------------------------
1380 !------------------------------Arguments--------------------------------
1381 ! number of current column
1382  INTEGER :: i
1383 ! number of atmospheric columns
1384  INTEGER :: pcols
1385 ! number of atmospheric columns
1386  INTEGER :: pver
1387 ! index of bottom (source) interface
1388  INTEGER :: kbot
1389 ! number of gravity waves to use
1390  INTEGER :: ngwv
1391 ! wave phase speeds
1392  REAL :: c(-ngwv:ngwv)
1393 ! midpoint zonal wind
1394  REAL :: u(pcols, pver)
1395  REAL :: ud(pcols, pver)
1396 ! midpoint meridional wind
1397  REAL :: v(pcols, pver)
1398  REAL :: vd(pcols, pver)
1399 ! midpoint temperatures
1400  REAL :: t(pcols, pver)
1401 ! midpoint pressures
1402  REAL :: pm(pcols, pver)
1403 ! interface pressures
1404  REAL :: pi(pcols, 0:pver)
1405 ! midpoint delta p (pi(k)-pi(k-1))
1406  REAL :: dpm(pcols, pver)
1407 ! 1. / (pi(k)-pi(k-1))
1408  REAL :: rdpm(pcols, pver)
1409 ! ln(interface pressures)
1410  REAL :: piln(pcols, 0:pver)
1411 ! latitude in radians for columns
1412  REAL :: rlat(pcols)
1413 ! top interface of low level stress divergence region
1414  INTEGER :: kldv
1415 ! min value of kldv
1416  INTEGER :: kldvmn
1417 ! index of top interface of source region
1418  INTEGER :: ksrc
1419 ! min value of ksrc
1420  INTEGER :: ksrcmn
1421 ! 1/dp across low level divergence region
1422  REAL :: rdpldv
1423 ! wave Reynolds stress
1424  REAL :: tau(-ngwv:ngwv, 0:pver)
1425 ! projection of wind at interfaces
1426  REAL :: ubi(0:pver)
1427  REAL :: ubid(0:pver)
1428 ! projection of wind at midpoints
1429  REAL :: ubm(pver)
1430  REAL :: ubmd(pver)
1431 ! unit vectors of source wind (x)
1432  REAL :: xv
1433  REAL :: xvd
1434 ! unit vectors of source wind (y)
1435  REAL :: yv
1436  REAL :: yvd
1437 !---------------------------Local storage-------------------------------
1438 ! loop indexes
1439  INTEGER :: k, l
1440 ! background stress at c=0
1441  REAL :: tauback
1442 ! u wind averaged over source region
1443  REAL :: usrc
1444  REAL :: usrcd
1445 ! v wind averaged over source region
1446  REAL :: vsrc
1447  REAL :: vsrcd
1448  REAL :: ubsrc
1449  REAL :: ubsrcd
1450 ! Used in lat dependence of GW spec.
1451  REAL :: al0
1452 ! Used in lat dependence of GW spec.
1453  REAL :: dlat0
1454  REAL :: latdeg
1455 ! The actual lat dependence of GW spec.
1456  REAL :: flat_gw
1457  REAL :: u_new, v_new, tau_new
1458  REAL :: arg1
1459  REAL :: arg1d
1460  INTRINSIC dexp
1461  INTRINSIC exp
1462  INTRINSIC max
1463  INTRINSIC abs
1464  REAL :: x1
1465  REAL :: x1d
1466  INTRINSIC dble
1467  REAL :: abs6
1468  REAL :: abs5
1469  REAL :: abs4
1470  REAL :: abs3
1471  REAL :: abs2
1472  REAL :: abs1
1473  INTRINSIC sqrt
1474  REAL :: y1
1475 !---------------------------------------------------------------------------
1476 ! Determine the source layer wind and unit vectors, then project winds.
1477 !---------------------------------------------------------------------------
1478 ! Just use the source level interface values for the source
1479 ! wind speed and direction (unit vector).
1480  ubm = 0.
1481  ubi = 0.
1482  xv = 0.
1483  yv = 0.
1484  tau = 0.
1485  ksrc = kbot
1486  kldv = kbot
1487  usrcd = 0.5*(ud(i, kbot+1)+ud(i, kbot))
1488  usrc = 0.5*(u(i, kbot+1)+u(i, kbot))
1489  vsrcd = 0.5*(vd(i, kbot+1)+vd(i, kbot))
1490  vsrc = 0.5*(v(i, kbot+1)+v(i, kbot))
1491  arg1d = 2*usrc*usrcd + 2*vsrc*vsrcd
1492  arg1 = usrc**2 + vsrc**2
1493  IF (arg1 .EQ. 0.0) THEN
1494  x1d = 0.0_8
1495  ELSE
1496  x1d = arg1d/(2.0*sqrt(arg1))
1497  END IF
1498  x1 = sqrt(arg1)
1499  y1 = sqrt(ubmc2mn)
1500  IF (x1 .LT. y1) THEN
1501  ubsrc = y1
1502  ubsrcd = 0.0_8
1503  ELSE
1504  ubsrcd = x1d
1505  ubsrc = x1
1506  END IF
1507  IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
1508  xv = 1.0
1509  yv = 0.0
1510  xvd = 0.0_8
1511  yvd = 0.0_8
1512  ubmd(:) = 0.0_8
1513  ELSE
1514  xvd = (usrcd*ubsrc-usrc*ubsrcd)/ubsrc**2
1515  xv = usrc/ubsrc
1516  yvd = (vsrcd*ubsrc-vsrc*ubsrcd)/ubsrc**2
1517  yv = vsrc/ubsrc
1518  ubmd(:) = 0.0_8
1519  END IF
1520 ! Project the local wind at midpoints onto the source wind.
1521  DO k=1,pver
1522  ubmd(k) = ud(i, k)*xv + u(i, k)*xvd + vd(i, k)*yv + v(i, k)*yvd
1523  ubm(k) = u(i, k)*xv + v(i, k)*yv
1524  END DO
1525 ! Compute the bottom interface wind projection using the midpoint winds.
1526  ubid(0) = ubmd(1)
1527  ubi(0) = ubm(1)
1528  DO k=1,pver
1529  ubid(k) = ubmd(k)
1530  ubi(k) = ubm(k)
1531  END DO
1532 !-----------------------------------------------------------------------
1533 ! Gravity wave sources
1534 !-----------------------------------------------------------------------
1535 ! Determine the background stress at c=0
1536  tauback = taubgnd*tauscal
1537 ! Include dependence on latitude:
1538  latdeg = rlat(i)*180./pi_gwd
1539 !
1540  IF (-15.3 .LT. latdeg .AND. latdeg .LT. 15.3) THEN
1541  IF (latdeg .GE. 0.) THEN
1542  abs1 = latdeg
1543  ELSE
1544  abs1 = -latdeg
1545  END IF
1546  arg1 = -(dble((abs1-3.)/8.0)**2)
1547  flat_gw = 1.2*exp(arg1)
1548  IF (latdeg .GE. 0.) THEN
1549  abs2 = latdeg
1550  ELSE
1551  abs2 = -latdeg
1552  END IF
1553  IF (flat_gw .LT. 1.2 .AND. abs2 .LE. 3.) flat_gw = 1.2
1554  ELSE IF (latdeg .GT. -31. .AND. latdeg .LE. -15.3) THEN
1555  flat_gw = 0.10
1556  ELSE IF (latdeg .LT. 31. .AND. latdeg .GE. 15.3) THEN
1557  flat_gw = 0.10
1558  ELSE IF (latdeg .GT. -60. .AND. latdeg .LE. -31.) THEN
1559  IF (latdeg .GE. 0.) THEN
1560  abs3 = latdeg
1561  ELSE
1562  abs3 = -latdeg
1563  END IF
1564  arg1 = -(dble((abs3-60.)/23.)**2)
1565  flat_gw = 0.50*exp(arg1)
1566  ELSE IF (latdeg .LT. 60. .AND. latdeg .GE. 31.) THEN
1567  IF (latdeg .GE. 0.) THEN
1568  abs4 = latdeg
1569  ELSE
1570  abs4 = -latdeg
1571  END IF
1572  arg1 = -(dble((abs4-60.)/23.)**2)
1573  flat_gw = 0.50*exp(arg1)
1574  ELSE IF (latdeg .LE. -60.) THEN
1575  IF (latdeg .GE. 0.) THEN
1576  abs5 = latdeg
1577  ELSE
1578  abs5 = -latdeg
1579  END IF
1580  arg1 = -(dble((abs5-60.)/70.)**2)
1581  flat_gw = 0.50*exp(arg1)
1582  ELSE IF (latdeg .GE. 60.) THEN
1583  IF (latdeg .GE. 0.) THEN
1584  abs6 = latdeg
1585  ELSE
1586  abs6 = -latdeg
1587  END IF
1588  arg1 = -(dble((abs6-60.)/70.)**2)
1589  flat_gw = 0.50*exp(arg1)
1590  END IF
1591  tauback = tauback*flat_gw
1592 ! Set the phase speeds and wave numbers in the direction of the source wind.
1593 ! Set the source stress magnitude (positive only, note that the sign of the
1594 ! stress is the same as (c-u).
1595  DO l=1,ngwv
1596  arg1 = -((c(l)/30.)**2)
1597  tau(l, kbot) = tauback*exp(arg1)
1598  tau(-l, kbot) = tau(l, kbot)
1599  END DO
1600  tau(0, kbot) = tauback
1601 ! Determine the min value of kldv and ksrc for limiting later loops
1602 ! and the pressure at the top interface of the low level stress divergence
1603 ! region.
1604  ksrcmn = pver
1605  kldvmn = pver
1606  RETURN
1607  END SUBROUTINE gw_bgnd_d
1608 !===============================================================================
1609  SUBROUTINE gw_bgnd(i, pcols, pver, c, u, v, t, pm, pi, dpm, rdpm, piln&
1610 & , rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, &
1611 & ngwv, kbot)
1612  IMPLICIT NONE
1613 !-----------------------------------------------------------------------
1614 ! Driver for multiple gravity wave drag parameterization.
1615 !
1616 ! The parameterization is assumed to operate only where water vapor
1617 ! concentrations are negligible in determining the density.
1618 !-----------------------------------------------------------------------
1619 !------------------------------Arguments--------------------------------
1620 ! number of current column
1621  INTEGER :: i
1622 ! number of atmospheric columns
1623  INTEGER :: pcols
1624 ! number of atmospheric columns
1625  INTEGER :: pver
1626 ! index of bottom (source) interface
1627  INTEGER :: kbot
1628 ! number of gravity waves to use
1629  INTEGER :: ngwv
1630 ! wave phase speeds
1631  REAL :: c(-ngwv:ngwv)
1632 ! midpoint zonal wind
1633  REAL :: u(pcols, pver)
1634 ! midpoint meridional wind
1635  REAL :: v(pcols, pver)
1636 ! midpoint temperatures
1637  REAL :: t(pcols, pver)
1638 ! midpoint pressures
1639  REAL :: pm(pcols, pver)
1640 ! interface pressures
1641  REAL :: pi(pcols, 0:pver)
1642 ! midpoint delta p (pi(k)-pi(k-1))
1643  REAL :: dpm(pcols, pver)
1644 ! 1. / (pi(k)-pi(k-1))
1645  REAL :: rdpm(pcols, pver)
1646 ! ln(interface pressures)
1647  REAL :: piln(pcols, 0:pver)
1648 ! latitude in radians for columns
1649  REAL :: rlat(pcols)
1650 ! top interface of low level stress divergence region
1651  INTEGER :: kldv
1652 ! min value of kldv
1653  INTEGER :: kldvmn
1654 ! index of top interface of source region
1655  INTEGER :: ksrc
1656 ! min value of ksrc
1657  INTEGER :: ksrcmn
1658 ! 1/dp across low level divergence region
1659  REAL :: rdpldv
1660 ! wave Reynolds stress
1661  REAL :: tau(-ngwv:ngwv, 0:pver)
1662 ! projection of wind at interfaces
1663  REAL :: ubi(0:pver)
1664 ! projection of wind at midpoints
1665  REAL :: ubm(pver)
1666 ! unit vectors of source wind (x)
1667  REAL :: xv
1668 ! unit vectors of source wind (y)
1669  REAL :: yv
1670 !---------------------------Local storage-------------------------------
1671 ! loop indexes
1672  INTEGER :: k, l
1673 ! background stress at c=0
1674  REAL :: tauback
1675 ! u wind averaged over source region
1676  REAL :: usrc
1677 ! v wind averaged over source region
1678  REAL :: vsrc
1679  REAL :: ubsrc
1680 ! Used in lat dependence of GW spec.
1681  REAL :: al0
1682 ! Used in lat dependence of GW spec.
1683  REAL :: dlat0
1684  REAL :: latdeg
1685 ! The actual lat dependence of GW spec.
1686  REAL :: flat_gw
1687  REAL :: u_new, v_new, tau_new
1688  REAL :: arg1
1689  INTRINSIC exp
1690  INTRINSIC max
1691  INTRINSIC abs
1692  REAL :: x1
1693  INTRINSIC dble
1694  REAL :: abs6
1695  REAL :: abs5
1696  REAL :: abs4
1697  REAL :: abs3
1698  REAL :: abs2
1699  REAL :: abs1
1700  INTRINSIC sqrt
1701  REAL :: y1
1702 !---------------------------------------------------------------------------
1703 ! Determine the source layer wind and unit vectors, then project winds.
1704 !---------------------------------------------------------------------------
1705 ! Just use the source level interface values for the source
1706 ! wind speed and direction (unit vector).
1707  ubm = 0.
1708  ubi = 0.
1709  xv = 0.
1710  yv = 0.
1711  tau = 0.
1712  ksrc = kbot
1713  kldv = kbot
1714  usrc = 0.5*(u(i, kbot+1)+u(i, kbot))
1715  vsrc = 0.5*(v(i, kbot+1)+v(i, kbot))
1716  arg1 = usrc**2 + vsrc**2
1717  x1 = sqrt(arg1)
1718  y1 = sqrt(ubmc2mn)
1719  IF (x1 .LT. y1) THEN
1720  ubsrc = y1
1721  ELSE
1722  ubsrc = x1
1723  END IF
1724  IF (usrc .EQ. 0. .AND. vsrc .EQ. 0.) THEN
1725  xv = 1.0
1726  yv = 0.0
1727  ELSE
1728  xv = usrc/ubsrc
1729  yv = vsrc/ubsrc
1730  END IF
1731 ! Project the local wind at midpoints onto the source wind.
1732  DO k=1,pver
1733  ubm(k) = u(i, k)*xv + v(i, k)*yv
1734  END DO
1735 ! Compute the bottom interface wind projection using the midpoint winds.
1736  ubi(0) = ubm(1)
1737  DO k=1,pver
1738  ubi(k) = ubm(k)
1739  END DO
1740 !-----------------------------------------------------------------------
1741 ! Gravity wave sources
1742 !-----------------------------------------------------------------------
1743 ! Determine the background stress at c=0
1744  tauback = taubgnd*tauscal
1745 ! Include dependence on latitude:
1746  latdeg = rlat(i)*180./pi_gwd
1747 !
1748  IF (-15.3 .LT. latdeg .AND. latdeg .LT. 15.3) THEN
1749  IF (latdeg .GE. 0.) THEN
1750  abs1 = latdeg
1751  ELSE
1752  abs1 = -latdeg
1753  END IF
1754  arg1 = -(dble((abs1-3.)/8.0)**2)
1755  flat_gw = 1.2*exp(arg1)
1756  IF (latdeg .GE. 0.) THEN
1757  abs2 = latdeg
1758  ELSE
1759  abs2 = -latdeg
1760  END IF
1761  IF (flat_gw .LT. 1.2 .AND. abs2 .LE. 3.) flat_gw = 1.2
1762  ELSE IF (latdeg .GT. -31. .AND. latdeg .LE. -15.3) THEN
1763  flat_gw = 0.10
1764  ELSE IF (latdeg .LT. 31. .AND. latdeg .GE. 15.3) THEN
1765  flat_gw = 0.10
1766  ELSE IF (latdeg .GT. -60. .AND. latdeg .LE. -31.) THEN
1767  IF (latdeg .GE. 0.) THEN
1768  abs3 = latdeg
1769  ELSE
1770  abs3 = -latdeg
1771  END IF
1772  arg1 = -(dble((abs3-60.)/23.)**2)
1773  flat_gw = 0.50*exp(arg1)
1774  ELSE IF (latdeg .LT. 60. .AND. latdeg .GE. 31.) THEN
1775  IF (latdeg .GE. 0.) THEN
1776  abs4 = latdeg
1777  ELSE
1778  abs4 = -latdeg
1779  END IF
1780  arg1 = -(dble((abs4-60.)/23.)**2)
1781  flat_gw = 0.50*exp(arg1)
1782  ELSE IF (latdeg .LE. -60.) THEN
1783  IF (latdeg .GE. 0.) THEN
1784  abs5 = latdeg
1785  ELSE
1786  abs5 = -latdeg
1787  END IF
1788  arg1 = -(dble((abs5-60.)/70.)**2)
1789  flat_gw = 0.50*exp(arg1)
1790  ELSE IF (latdeg .GE. 60.) THEN
1791  IF (latdeg .GE. 0.) THEN
1792  abs6 = latdeg
1793  ELSE
1794  abs6 = -latdeg
1795  END IF
1796  arg1 = -(dble((abs6-60.)/70.)**2)
1797  flat_gw = 0.50*exp(arg1)
1798  END IF
1799  tauback = tauback*flat_gw
1800 ! Set the phase speeds and wave numbers in the direction of the source wind.
1801 ! Set the source stress magnitude (positive only, note that the sign of the
1802 ! stress is the same as (c-u).
1803  DO l=1,ngwv
1804  arg1 = -((c(l)/30.)**2)
1805  tau(l, kbot) = tauback*exp(arg1)
1806  tau(-l, kbot) = tau(l, kbot)
1807  END DO
1808  tau(0, kbot) = tauback
1809 ! Determine the min value of kldv and ksrc for limiting later loops
1810 ! and the pressure at the top interface of the low level stress divergence
1811 ! region.
1812  ksrcmn = pver
1813  kldvmn = pver
1814  RETURN
1815  END SUBROUTINE gw_bgnd
1816 ! Differentiation of gw_drag_prof in forward (tangent) mode:
1817 ! variations of output variables: ut vt
1818 ! with respect to input variables: tau xv ut yv vt ubi
1819 !===============================================================================
1820  SUBROUTINE gw_drag_prof_d(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u&
1821 & , v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, &
1822 & dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, ubm&
1823 & , xv, xvd, yv, yvd, ut, utd, vt, vtd, tt, taugwx, taugwy, fegw, &
1824 & fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y, effgw)
1825  IMPLICIT NONE
1826 !-----------------------------------------------------------------------
1827 ! Solve for the drag profile from the multiple gravity wave drag
1828 ! parameterization.
1829 ! 1. scan up from the wave source to determine the stress profile
1830 ! 2. scan down the stress profile to determine the tendencies
1831 ! => apply bounds to the tendency
1832 ! a. from wkb solution
1833 ! b. from computational stability constraints
1834 ! => adjust stress on interface below to reflect actual bounded tendency
1835 !-----------------------------------------------------------------------
1836 !------------------------------Arguments--------------------------------
1837 ! number of current column
1838  INTEGER, INTENT(IN) :: i
1839 ! number of atmospheric columns
1840  INTEGER, INTENT(IN) :: pcols
1841 ! number of atmospheric columns
1842  INTEGER, INTENT(IN) :: pver
1843 ! index of bottom (source) interface
1844  INTEGER, INTENT(IN) :: kbot
1845 ! index of top interface of gwd region
1846  INTEGER, INTENT(IN) :: ktop
1847 ! number of gravity waves possible
1848  INTEGER, INTENT(IN) :: pgwv
1849 ! number of gravity waves to use
1850  INTEGER, INTENT(IN) :: ngwv
1851 ! top interface of low level stress divergence region
1852  INTEGER, INTENT(IN) :: kldv
1853 ! min value of kldv
1854  INTEGER, INTENT(IN) :: kldvmn
1855 ! index of top interface of source region
1856  INTEGER, INTENT(IN) :: ksrc
1857 ! min value of ksrc
1858  INTEGER, INTENT(IN) :: ksrcmn
1859 ! wave phase speeds
1860  REAL :: c(-pgwv:pgwv)
1861 ! midpoint zonal wind
1862  REAL :: u(pcols, pver)
1863 ! midpoint meridional wind
1864  REAL :: v(pcols, pver)
1865 ! midpoint temperatures
1866  REAL :: t(pcols, pver)
1867 ! interface pressures
1868  REAL :: pi(pcols, 0:pver)
1869 ! midpoint delta p (pi(k)-pi(k-1))
1870  REAL :: dpm(pcols, pver)
1871 ! 1. / (pi(k)-pi(k-1))
1872  REAL :: rdpm(pcols, pver)
1873 ! ln(interface pressures)
1874  REAL :: piln(pcols, 0:pver)
1875  REAL :: rlat(pcols)
1876 ! interface density
1877  REAL :: rhoi(0:pver)
1878 ! interface Brunt-Vaisalla frequency
1879  REAL :: ni(0:pver)
1880 ! interface temperature
1881  REAL :: ti(0:pver)
1882 ! midpoint Brunt-Vaisalla frequency
1883  REAL :: nm(pver)
1884 ! time step
1885  REAL :: dt
1886 ! newtonian cooling coefficients
1887  REAL :: alpha(0:pver)
1888 ! newtonian cooling coefficients
1889  REAL :: dback(0:pver)
1890 ! 1/dp across low level divergence region
1891  REAL :: rdpldv
1892 ! projection of wind at interfaces
1893  REAL :: ubi(0:pver)
1894  REAL :: ubid(0:pver)
1895 ! projection of wind at midpoints
1896  REAL :: ubm(pver)
1897 ! unit vectors of source wind (x)
1898  REAL :: xv
1899  REAL :: xvd
1900 ! unit vectors of source wind (y)
1901  REAL :: yv
1902  REAL :: yvd
1903 ! tendency efficiency for gwd
1904  REAL :: effgw(pcols)
1905 ! wave Reynolds stress
1906  REAL :: tau(-pgwv:pgwv, 0:pver)
1907  REAL :: taud(-pgwv:pgwv, 0:pver)
1908 ! zonal wind tendency
1909  REAL :: ut(pver)
1910  REAL :: utd(pver)
1911 ! meridional wind tendency
1912  REAL :: vt(pver)
1913  REAL :: vtd(pver)
1914 ! temperature tendency
1915  REAL :: tt(pver)
1916 ! Total zonal GW momentum flux
1917  REAL :: taugwx(pcols, 0:pver)
1918 ! Total meridional GW momentum flux
1919  REAL :: taugwy(pcols, 0:pver)
1920 ! Total GW energy flux
1921  REAL :: fegw(pcols, 0:pver)
1922 ! Total GW pseudo energy flux
1923  REAL :: fepgw(pcols, 0:pver)
1924 ! Total U tendency below launch level
1925  REAL :: dusrc(pver)
1926 ! Total V tendency below launch level
1927  REAL :: dvsrc(pver)
1928 ! Total V tendency below launch level
1929  REAL :: dtsrc(pver)
1930 ! c=0 sfc. stress (zonal)
1931  REAL :: tau0x
1932 ! c=0 sfc. stress (meridional)
1933  REAL :: tau0y
1934 !---------------------------Local storage-------------------------------
1935 ! loop indexes
1936  INTEGER :: k, l
1937 ! fraction of dsat to use
1938  REAL :: dscal
1939 ! imaginary part of vertical wavenumber
1940  REAL :: mi
1941 ! stress after damping
1942  REAL :: taudmp
1943 !real :: taumax(pcols) ! max(tau) for any l
1944 ! saturation stress
1945  REAL :: tausat
1946  REAL :: tausatd
1947 ! (ub-c)
1948  REAL :: ubmc
1949  REAL :: ubmcd
1950 ! (ub-c)**2
1951  REAL :: ubmc2
1952 ! ubar tendency
1953  REAL :: ubt
1954  REAL :: ubtd
1955 ! tbar tendency
1956  REAL :: tbt
1957 ! ubar tendency from wave l
1958  REAL :: ubtl
1959  REAL :: ubtld
1960 ! saturation tendency
1961  REAL :: ubtlsat
1962 ! layer pressure
1963  REAL :: pm
1964 ! layer density
1965  REAL :: rhom
1966 ! launch level height
1967  REAL :: zlb
1968 ! c-u
1969  REAL :: cmu
1970  REAL :: dzm, hscal, tautmp
1971  REAL :: utl
1972  REAL :: utld
1973  REAL :: ttl
1974 ! zonal pseudomomentum flux spectrum
1975  REAL :: fpmx
1976 ! meridional pseudomomentum flux spectrum
1977  REAL :: fpmy
1978 ! energy flux (p'w') spectrum
1979  REAL :: fe
1980 ! pseudoenergy flux (p'w'+U rho u'w') spectrum
1981  REAL :: fpe
1982  REAL :: fpml
1983  REAL :: fpmt
1984  REAL :: fpel
1985  REAL :: fpet
1986  REAL :: dusrcl
1987  REAL :: dvsrcl
1988  REAL :: dtsrcl
1989  REAL :: zi
1990  REAL :: effkwvmap
1991  REAL :: zfac
1992  REAL :: uhtmax
1993  REAL :: uhtmaxd
1994  REAL :: utfac
1995  REAL :: utfacd
1996  REAL :: tau_min, cmu_, t_new
1997  REAL :: tau_mind
1998  REAL :: arg1
1999  REAL :: arg1d
2000  INTRINSIC max
2001  INTRINSIC sign
2002  INTRINSIC abs
2003  REAL :: x2
2004  REAL :: x2d
2005  REAL :: x1
2006  REAL :: x1d
2007  INTRINSIC sqrt
2008 ! Initialize gravity wave drag tendencies to zero
2009  DO l=-ngwv,ngwv
2010  DO k=pver-1,ktop,-1
2011  IF (k .LE. kbot - 1) THEN
2012  taud(l, k) = 0.0_8
2013  tau(l, k) = 0.
2014  END IF
2015  END DO
2016  END DO
2017  DO k=1,pver
2018  utd(k) = 0.0_8
2019  ut(k) = 0.
2020  vtd(k) = 0.0_8
2021  vt(k) = 0.
2022  tt(k) = 0.
2023  dusrc(k) = 0.
2024  dvsrc(k) = 0.
2025  dtsrc(k) = 0.
2026  END DO
2027 ! Initialize total momentum and energy fluxes to zero
2028  DO k=0,pver
2029  taugwx(i, k) = 0.
2030  taugwy(i, k) = 0.
2031  fegw(i, k) = 0.
2032  fepgw(i, k) = 0.
2033  END DO
2034 ! Initialize surface wave stress at c = 0 to zero
2035  tau0x = 0.
2036  tau0y = 0.
2037 !---------------------------------------------------------------------------
2038 ! Compute the stress profiles and diffusivities
2039 !---------------------------------------------------------------------------
2040 ! Determine the absolute value of the saturation stress and the diffusivity
2041 ! for each wave.
2042 ! Define critical levels where the sign of (u-c) changes between interfaces.
2043 ! Loop from bottom to top to get stress profiles
2044  DO l=-ngwv,ngwv
2045  DO k=pver-1,ktop,-1
2046  IF (k .LE. kbot - 1) THEN
2047 !!!!!!!!jk d = dback(k)
2048  ubmcd = ubid(k)
2049  ubmc = ubi(k) - c(l)
2050  IF (ngwv .GT. 0) THEN
2051  CALL get_effkwvmap_1(effkwvmap, rlat(i))
2052  IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
2053 & .LT. 15.0) effkwvmap = fcrit2*kwvbeq
2054  ELSE
2055  IF (pi(i, k) .LT. 1000.0) THEN
2056  zfac = (pi(i, k)/1000.0)**3
2057  ELSE
2058  zfac = 1.0
2059  END IF
2060  IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
2061  CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
2062  ELSE
2063  CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
2064  END IF
2065  END IF
2066  x1d = effkwvmap*rhoi(k)*3*ubmc**2*ubmcd/(2.*ni(k))
2067  x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
2068  IF (x1 .GE. 0.) THEN
2069  tausatd = x1d
2070  tausat = x1
2071  ELSE
2072  tausatd = -x1d
2073  tausat = -x1
2074  END IF
2075  IF (tausat .LE. taumin) THEN
2076  tausat = 0.0
2077  tausatd = 0.0_8
2078  END IF
2079  IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) THEN
2080  tausat = 0.0
2081  tausatd = 0.0_8
2082  END IF
2083  IF (k .EQ. ktop) THEN
2084  tausat = 0.
2085  tausatd = 0.0_8
2086  END IF
2087 !
2088  IF (k .EQ. ktop + 1) THEN
2089  tausatd = 0.02*tausatd
2090  tausat = tausat*0.02
2091  END IF
2092  IF (k .EQ. ktop + 2) THEN
2093  tausatd = 0.05*tausatd
2094  tausat = tausat*0.05
2095  END IF
2096  IF (k .EQ. ktop + 3) THEN
2097  tausatd = 0.10*tausatd
2098  tausat = tausat*0.10
2099  END IF
2100  IF (k .EQ. ktop + 4) THEN
2101  tausatd = 0.20*tausatd
2102  tausat = tausat*0.20
2103  END IF
2104  IF (k .EQ. ktop + 5) THEN
2105  tausatd = 0.50*tausatd
2106  tausat = tausat*0.50
2107  END IF
2108 !
2109  tau_mind = tausatd
2110  tau_min = tausat
2111  IF (tau_min .GT. tau(l, k+1)) THEN
2112  tau_mind = taud(l, k+1)
2113  tau_min = tau(l, k+1)
2114  END IF
2115  taud(l, k) = tau_mind
2116  tau(l, k) = tau_min
2117  END IF
2118  END DO
2119  END DO
2120 !---------------------------------------------------------------------------
2121 ! Compute the tendencies from the stress divergence.
2122 !---------------------------------------------------------------------------
2123 ! Accumulate the mean wind tendency over wavenumber.
2124 ! Loop over levels from top to bottom
2125  DO k=ktop+1,pver
2126  ubt = 0.0
2127  tbt = 0.0
2128  ubtd = 0.0_8
2129  DO l=-ngwv,ngwv
2130  IF (k .LE. kbot) THEN
2131 ! Determine the wind tendency including excess stress carried down from above.
2132  ubtld = mapl_grav*rdpm(i, k)*(taud(l, k)-taud(l, k-1))
2133  ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
2134 ! Calculate the sign of wind tendency
2135  utld = ubtld*sign(1.d0, ubtl*(c(l)-ubi(k)))
2136  utl = sign(ubtl, c(l) - ubi(k))
2137 ! Accumulate the mean wind tendency over wavenumber.
2138  ubtd = ubtd + utld
2139  ubt = ubt + utl
2140 ! Calculate irreversible temperature tendency associated with gravity wave breaking.
2141  ttl = (c(l)-ubm(k))*utl/mapl_cp
2142 ! Adding frictional heating associated with the GW momentum forcing
2143  tbt = tbt + ttl
2144  END IF
2145  END DO
2146 ! Project the mean wind tendency onto the components and scale by "efficiency".
2147  IF (k .LE. kbot) THEN
2148  utd(k) = effgw(i)*(ubtd*xv+ubt*xvd)
2149  ut(k) = ubt*xv*effgw(i)
2150  vtd(k) = effgw(i)*(ubtd*yv+ubt*yvd)
2151  vt(k) = ubt*yv*effgw(i)
2152  tt(k) = tbt*effgw(i)
2153  END IF
2154  END DO
2155 !-----------------------------------------------------------------------
2156 ! Calculates wind and temperature tendencies below launch level for
2157 ! energy and momentum conservation (does nothing for orographic GWs).
2158 !-----------------------------------------------------------------------
2159 !----kimmmmm here
2160 ! Calculate launch level height
2161  zlb = 0.
2162  DO k=ktop+1,pver
2163  IF (k .GE. kbot + 1) THEN
2164 ! Define layer pressure and density
2165  pm = (pi(i, k-1)+pi(i, k))*0.5
2166  CALL get_ti(t_new, t(i, k))
2167  rhom = pm/(mapl_rgas*t_new)
2168  zlb = zlb + dpm(i, k)/mapl_grav/rhom
2169  END IF
2170  END DO
2171 !-----------------------------------------------------------------------
2172 ! Calculates energy and momentum flux profiles
2173 !-----------------------------------------------------------------------
2174  DO l=-ngwv,ngwv
2175  DO k=ktop,pver
2176  IF (k .LE. kbot) THEN
2177  cmu = c(l) - ubi(k)
2178  CALL get_cmu(cmu_, cmu)
2179  fpmx = cmu_*tau(l, k)*xv*effgw(i)
2180  fpmy = cmu_*tau(l, k)*yv*effgw(i)
2181  fe = cmu*cmu_*tau(l, k)*effgw(i)
2182  fpe = c(l)*cmu_*tau(l, k)*effgw(i)
2183  IF (k .EQ. kbot) THEN
2184  fpml = fpmx*xv + fpmy*yv
2185  fpel = fpe
2186  END IF
2187  IF (k .EQ. ktop) THEN
2188  fpmt = fpmx*xv + fpmy*yv
2189  fpet = fpe
2190  END IF
2191 ! Record outputs for GW fluxes
2192  taugwx(i, k) = taugwx(i, k) + fpmx
2193  taugwy(i, k) = taugwy(i, k) + fpmy
2194  fegw(i, k) = fegw(i, k) + fe
2195  fepgw(i, k) = fepgw(i, k) + fpe
2196  END IF
2197  END DO
2198  DO k=ktop+1,pver
2199  IF (k .GE. kbot + 1) THEN
2200 ! Define layer pressure and density
2201  pm = (pi(i, k-1)+pi(i, k))*0.5
2202  CALL get_ti(t_new, t(i, k))
2203  rhom = pm/(mapl_rgas*t_new)
2204  dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
2205  dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
2206  dtsrcl = -((fpel-fpet-ubm(k)*(fpml-fpmt))/(rhom*zlb*mapl_cp))
2207 ! Add sub-source wind and temperature tendencies
2208  dusrc(k) = dusrc(k) + dusrcl
2209  dvsrc(k) = dvsrc(k) + dvsrcl
2210  dtsrc(k) = dtsrc(k) + dtsrcl
2211  END IF
2212  END DO
2213  END DO
2214 ! For orographic waves, sub-source tendencies are set equal to zero.
2215  DO k=1,pver
2216  IF (ngwv .EQ. 0) THEN
2217  dusrc(k) = 0.0
2218  dvsrc(k) = 0.0
2219  dtsrc(k) = 0.0
2220  END IF
2221  END DO
2222 !-----------------------------------------------------------------------
2223 ! Adjust efficiency factor to prevent unrealistically strong forcing
2224 !-----------------------------------------------------------------------
2225  uhtmax = 0.0
2226  utfac = 1.0
2227  uhtmaxd = 0.0_8
2228  DO k=1,pver
2229  arg1d = 2*ut(k)*utd(k) + 2*vt(k)*vtd(k)
2230  arg1 = ut(k)**2 + vt(k)**2
2231  IF (arg1 .EQ. 0.0) THEN
2232  x2d = 0.0_8
2233  ELSE
2234  x2d = arg1d/(2.0*sqrt(arg1))
2235  END IF
2236  x2 = sqrt(arg1)
2237  IF (x2 .LT. uhtmax) THEN
2238  uhtmax = uhtmax
2239  ELSE
2240  uhtmaxd = x2d
2241  uhtmax = x2
2242  END IF
2243  END DO
2244  IF (uhtmax .GT. tndmax) THEN
2245  utfacd = -(tndmax*uhtmaxd/uhtmax**2)
2246  utfac = tndmax/uhtmax
2247  ELSE
2248  utfacd = 0.0_8
2249  END IF
2250 !jkim effgw(i) = effgw(i)*utfac
2251  DO k=1,pver
2252  utd(k) = utd(k)*utfac + ut(k)*utfacd
2253  ut(k) = ut(k)*utfac
2254  vtd(k) = vtd(k)*utfac + vt(k)*utfacd
2255  vt(k) = vt(k)*utfac
2256  tt(k) = tt(k)*utfac
2257  dusrc(k) = dusrc(k)*utfac
2258  dvsrc(k) = dvsrc(k)*utfac
2259  dtsrc(k) = dtsrc(k)*utfac
2260  END DO
2261 !-----------------------------------------------------------------------
2262 ! Project the c=0 stress (scaled) in the direction of the source wind
2263 ! for recording on the output file.
2264 !-----------------------------------------------------------------------
2265  tau0x = tau(0, kbot)*xv*effgw(i)*utfac
2266  tau0y = tau(0, kbot)*yv*effgw(i)*utfac
2267  RETURN
2268  END SUBROUTINE gw_drag_prof_d
2269 !===============================================================================
2270  SUBROUTINE gw_drag_prof(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u, &
2271 & v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, dback&
2272 & , kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut, vt&
2273 & , tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y&
2274 & , effgw)
2275  IMPLICIT NONE
2276 !-----------------------------------------------------------------------
2277 ! Solve for the drag profile from the multiple gravity wave drag
2278 ! parameterization.
2279 ! 1. scan up from the wave source to determine the stress profile
2280 ! 2. scan down the stress profile to determine the tendencies
2281 ! => apply bounds to the tendency
2282 ! a. from wkb solution
2283 ! b. from computational stability constraints
2284 ! => adjust stress on interface below to reflect actual bounded tendency
2285 !-----------------------------------------------------------------------
2286 !------------------------------Arguments--------------------------------
2287 ! number of current column
2288  INTEGER, INTENT(IN) :: i
2289 ! number of atmospheric columns
2290  INTEGER, INTENT(IN) :: pcols
2291 ! number of atmospheric columns
2292  INTEGER, INTENT(IN) :: pver
2293 ! index of bottom (source) interface
2294  INTEGER, INTENT(IN) :: kbot
2295 ! index of top interface of gwd region
2296  INTEGER, INTENT(IN) :: ktop
2297 ! number of gravity waves possible
2298  INTEGER, INTENT(IN) :: pgwv
2299 ! number of gravity waves to use
2300  INTEGER, INTENT(IN) :: ngwv
2301 ! top interface of low level stress divergence region
2302  INTEGER, INTENT(IN) :: kldv
2303 ! min value of kldv
2304  INTEGER, INTENT(IN) :: kldvmn
2305 ! index of top interface of source region
2306  INTEGER, INTENT(IN) :: ksrc
2307 ! min value of ksrc
2308  INTEGER, INTENT(IN) :: ksrcmn
2309 ! wave phase speeds
2310  REAL :: c(-pgwv:pgwv)
2311 ! midpoint zonal wind
2312  REAL :: u(pcols, pver)
2313 ! midpoint meridional wind
2314  REAL :: v(pcols, pver)
2315 ! midpoint temperatures
2316  REAL :: t(pcols, pver)
2317 ! interface pressures
2318  REAL :: pi(pcols, 0:pver)
2319 ! midpoint delta p (pi(k)-pi(k-1))
2320  REAL :: dpm(pcols, pver)
2321 ! 1. / (pi(k)-pi(k-1))
2322  REAL :: rdpm(pcols, pver)
2323 ! ln(interface pressures)
2324  REAL :: piln(pcols, 0:pver)
2325  REAL :: rlat(pcols)
2326 ! interface density
2327  REAL :: rhoi(0:pver)
2328 ! interface Brunt-Vaisalla frequency
2329  REAL :: ni(0:pver)
2330 ! interface temperature
2331  REAL :: ti(0:pver)
2332 ! midpoint Brunt-Vaisalla frequency
2333  REAL :: nm(pver)
2334 ! time step
2335  REAL :: dt
2336 ! newtonian cooling coefficients
2337  REAL :: alpha(0:pver)
2338 ! newtonian cooling coefficients
2339  REAL :: dback(0:pver)
2340 ! 1/dp across low level divergence region
2341  REAL :: rdpldv
2342 ! projection of wind at interfaces
2343  REAL :: ubi(0:pver)
2344 ! projection of wind at midpoints
2345  REAL :: ubm(pver)
2346 ! unit vectors of source wind (x)
2347  REAL :: xv
2348 ! unit vectors of source wind (y)
2349  REAL :: yv
2350 ! tendency efficiency for gwd
2351  REAL :: effgw(pcols)
2352 ! wave Reynolds stress
2353  REAL :: tau(-pgwv:pgwv, 0:pver)
2354 ! zonal wind tendency
2355  REAL :: ut(pver)
2356 ! meridional wind tendency
2357  REAL :: vt(pver)
2358 ! temperature tendency
2359  REAL :: tt(pver)
2360 ! Total zonal GW momentum flux
2361  REAL :: taugwx(pcols, 0:pver)
2362 ! Total meridional GW momentum flux
2363  REAL :: taugwy(pcols, 0:pver)
2364 ! Total GW energy flux
2365  REAL :: fegw(pcols, 0:pver)
2366 ! Total GW pseudo energy flux
2367  REAL :: fepgw(pcols, 0:pver)
2368 ! Total U tendency below launch level
2369  REAL :: dusrc(pver)
2370 ! Total V tendency below launch level
2371  REAL :: dvsrc(pver)
2372 ! Total V tendency below launch level
2373  REAL :: dtsrc(pver)
2374 ! c=0 sfc. stress (zonal)
2375  REAL :: tau0x
2376 ! c=0 sfc. stress (meridional)
2377  REAL :: tau0y
2378 !---------------------------Local storage-------------------------------
2379 ! loop indexes
2380  INTEGER :: k, l
2381 ! fraction of dsat to use
2382  REAL :: dscal
2383 ! imaginary part of vertical wavenumber
2384  REAL :: mi
2385 ! stress after damping
2386  REAL :: taudmp
2387 !real :: taumax(pcols) ! max(tau) for any l
2388 ! saturation stress
2389  REAL :: tausat
2390 ! (ub-c)
2391  REAL :: ubmc
2392 ! (ub-c)**2
2393  REAL :: ubmc2
2394 ! ubar tendency
2395  REAL :: ubt
2396 ! tbar tendency
2397  REAL :: tbt
2398 ! ubar tendency from wave l
2399  REAL :: ubtl
2400 ! saturation tendency
2401  REAL :: ubtlsat
2402 ! layer pressure
2403  REAL :: pm
2404 ! layer density
2405  REAL :: rhom
2406 ! launch level height
2407  REAL :: zlb
2408 ! c-u
2409  REAL :: cmu
2410  REAL :: dzm, hscal, tautmp
2411  REAL :: utl
2412  REAL :: ttl
2413 ! zonal pseudomomentum flux spectrum
2414  REAL :: fpmx
2415 ! meridional pseudomomentum flux spectrum
2416  REAL :: fpmy
2417 ! energy flux (p'w') spectrum
2418  REAL :: fe
2419 ! pseudoenergy flux (p'w'+U rho u'w') spectrum
2420  REAL :: fpe
2421  REAL :: fpml
2422  REAL :: fpmt
2423  REAL :: fpel
2424  REAL :: fpet
2425  REAL :: dusrcl
2426  REAL :: dvsrcl
2427  REAL :: dtsrcl
2428  REAL :: zi
2429  REAL :: effkwvmap
2430  REAL :: zfac
2431  REAL :: uhtmax
2432  REAL :: utfac
2433  REAL :: tau_min, cmu_, t_new
2434  REAL :: arg1
2435  INTRINSIC max
2436  INTRINSIC sign
2437  INTRINSIC abs
2438  REAL :: x2
2439  REAL :: x1
2440  INTRINSIC sqrt
2441 ! Initialize gravity wave drag tendencies to zero
2442  DO l=-ngwv,ngwv
2443  DO k=pver-1,ktop,-1
2444  IF (k .LE. kbot - 1) tau(l, k) = 0.
2445  END DO
2446  END DO
2447  DO k=1,pver
2448  ut(k) = 0.
2449  vt(k) = 0.
2450  tt(k) = 0.
2451  dusrc(k) = 0.
2452  dvsrc(k) = 0.
2453  dtsrc(k) = 0.
2454  END DO
2455 ! Initialize total momentum and energy fluxes to zero
2456  DO k=0,pver
2457  taugwx(i, k) = 0.
2458  taugwy(i, k) = 0.
2459  fegw(i, k) = 0.
2460  fepgw(i, k) = 0.
2461  END DO
2462 ! Initialize surface wave stress at c = 0 to zero
2463  tau0x = 0.
2464  tau0y = 0.
2465 !---------------------------------------------------------------------------
2466 ! Compute the stress profiles and diffusivities
2467 !---------------------------------------------------------------------------
2468 ! Determine the absolute value of the saturation stress and the diffusivity
2469 ! for each wave.
2470 ! Define critical levels where the sign of (u-c) changes between interfaces.
2471 ! Loop from bottom to top to get stress profiles
2472  DO l=-ngwv,ngwv
2473  DO k=pver-1,ktop,-1
2474  IF (k .LE. kbot - 1) THEN
2475 !!!!!!!!jk d = dback(k)
2476  ubmc = ubi(k) - c(l)
2477  IF (ngwv .GT. 0) THEN
2478  CALL get_effkwvmap_1(effkwvmap, rlat(i))
2479  IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
2480 & .LT. 15.0) effkwvmap = fcrit2*kwvbeq
2481  ELSE
2482  IF (pi(i, k) .LT. 1000.0) THEN
2483  zfac = (pi(i, k)/1000.0)**3
2484  ELSE
2485  zfac = 1.0
2486  END IF
2487  IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
2488  CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
2489  ELSE
2490  CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
2491  END IF
2492  END IF
2493  x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
2494  IF (x1 .GE. 0.) THEN
2495  tausat = x1
2496  ELSE
2497  tausat = -x1
2498  END IF
2499  IF (tausat .LE. taumin) tausat = 0.0
2500  IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) tausat = 0.0
2501  IF (k .EQ. ktop) tausat = 0.
2502 !
2503  IF (k .EQ. ktop + 1) tausat = tausat*0.02
2504  IF (k .EQ. ktop + 2) tausat = tausat*0.05
2505  IF (k .EQ. ktop + 3) tausat = tausat*0.10
2506  IF (k .EQ. ktop + 4) tausat = tausat*0.20
2507  IF (k .EQ. ktop + 5) tausat = tausat*0.50
2508 !
2509  tau_min = tausat
2510  IF (tau_min .GT. tau(l, k+1)) tau_min = tau(l, k+1)
2511  tau(l, k) = tau_min
2512  END IF
2513  END DO
2514  END DO
2515 !---------------------------------------------------------------------------
2516 ! Compute the tendencies from the stress divergence.
2517 !---------------------------------------------------------------------------
2518 ! Accumulate the mean wind tendency over wavenumber.
2519 ! Loop over levels from top to bottom
2520  DO k=ktop+1,pver
2521  ubt = 0.0
2522  tbt = 0.0
2523  DO l=-ngwv,ngwv
2524  IF (k .LE. kbot) THEN
2525 ! Determine the wind tendency including excess stress carried down from above.
2526  ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
2527 ! Calculate the sign of wind tendency
2528  utl = sign(ubtl, c(l) - ubi(k))
2529 ! Accumulate the mean wind tendency over wavenumber.
2530  ubt = ubt + utl
2531 ! Calculate irreversible temperature tendency associated with gravity wave breaking.
2532  ttl = (c(l)-ubm(k))*utl/mapl_cp
2533 ! Adding frictional heating associated with the GW momentum forcing
2534  tbt = tbt + ttl
2535  END IF
2536  END DO
2537 ! Project the mean wind tendency onto the components and scale by "efficiency".
2538  IF (k .LE. kbot) THEN
2539  ut(k) = ubt*xv*effgw(i)
2540  vt(k) = ubt*yv*effgw(i)
2541  tt(k) = tbt*effgw(i)
2542  END IF
2543  END DO
2544 !-----------------------------------------------------------------------
2545 ! Calculates wind and temperature tendencies below launch level for
2546 ! energy and momentum conservation (does nothing for orographic GWs).
2547 !-----------------------------------------------------------------------
2548 !----kimmmmm here
2549 ! Calculate launch level height
2550  zlb = 0.
2551  DO k=ktop+1,pver
2552  IF (k .GE. kbot + 1) THEN
2553 ! Define layer pressure and density
2554  pm = (pi(i, k-1)+pi(i, k))*0.5
2555  CALL get_ti(t_new, t(i, k))
2556  rhom = pm/(mapl_rgas*t_new)
2557  zlb = zlb + dpm(i, k)/mapl_grav/rhom
2558  END IF
2559  END DO
2560 !-----------------------------------------------------------------------
2561 ! Calculates energy and momentum flux profiles
2562 !-----------------------------------------------------------------------
2563  DO l=-ngwv,ngwv
2564  DO k=ktop,pver
2565  IF (k .LE. kbot) THEN
2566  cmu = c(l) - ubi(k)
2567  CALL get_cmu(cmu_, cmu)
2568  fpmx = cmu_*tau(l, k)*xv*effgw(i)
2569  fpmy = cmu_*tau(l, k)*yv*effgw(i)
2570  fe = cmu*cmu_*tau(l, k)*effgw(i)
2571  fpe = c(l)*cmu_*tau(l, k)*effgw(i)
2572  IF (k .EQ. kbot) THEN
2573  fpml = fpmx*xv + fpmy*yv
2574  fpel = fpe
2575  END IF
2576  IF (k .EQ. ktop) THEN
2577  fpmt = fpmx*xv + fpmy*yv
2578  fpet = fpe
2579  END IF
2580 ! Record outputs for GW fluxes
2581  taugwx(i, k) = taugwx(i, k) + fpmx
2582  taugwy(i, k) = taugwy(i, k) + fpmy
2583  fegw(i, k) = fegw(i, k) + fe
2584  fepgw(i, k) = fepgw(i, k) + fpe
2585  END IF
2586  END DO
2587  DO k=ktop+1,pver
2588  IF (k .GE. kbot + 1) THEN
2589 ! Define layer pressure and density
2590  pm = (pi(i, k-1)+pi(i, k))*0.5
2591  CALL get_ti(t_new, t(i, k))
2592  rhom = pm/(mapl_rgas*t_new)
2593  dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
2594  dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
2595  dtsrcl = -((fpel-fpet-ubm(k)*(fpml-fpmt))/(rhom*zlb*mapl_cp))
2596 ! Add sub-source wind and temperature tendencies
2597  dusrc(k) = dusrc(k) + dusrcl
2598  dvsrc(k) = dvsrc(k) + dvsrcl
2599  dtsrc(k) = dtsrc(k) + dtsrcl
2600  END IF
2601  END DO
2602  END DO
2603 ! For orographic waves, sub-source tendencies are set equal to zero.
2604  DO k=1,pver
2605  IF (ngwv .EQ. 0) THEN
2606  dusrc(k) = 0.0
2607  dvsrc(k) = 0.0
2608  dtsrc(k) = 0.0
2609  END IF
2610  END DO
2611 !-----------------------------------------------------------------------
2612 ! Adjust efficiency factor to prevent unrealistically strong forcing
2613 !-----------------------------------------------------------------------
2614  uhtmax = 0.0
2615  utfac = 1.0
2616  DO k=1,pver
2617  arg1 = ut(k)**2 + vt(k)**2
2618  x2 = sqrt(arg1)
2619  IF (x2 .LT. uhtmax) THEN
2620  uhtmax = uhtmax
2621  ELSE
2622  uhtmax = x2
2623  END IF
2624  END DO
2625  IF (uhtmax .GT. tndmax) utfac = tndmax/uhtmax
2626 !jkim effgw(i) = effgw(i)*utfac
2627  DO k=1,pver
2628  ut(k) = ut(k)*utfac
2629  vt(k) = vt(k)*utfac
2630  tt(k) = tt(k)*utfac
2631  dusrc(k) = dusrc(k)*utfac
2632  dvsrc(k) = dvsrc(k)*utfac
2633  dtsrc(k) = dtsrc(k)*utfac
2634  END DO
2635 !-----------------------------------------------------------------------
2636 ! Project the c=0 stress (scaled) in the direction of the source wind
2637 ! for recording on the output file.
2638 !-----------------------------------------------------------------------
2639  tau0x = tau(0, kbot)*xv*effgw(i)*utfac
2640  tau0y = tau(0, kbot)*yv*effgw(i)*utfac
2641  RETURN
2642  END SUBROUTINE gw_drag_prof
2643 ! Differentiation of gw_drag_prof_bgnd in forward (tangent) mode:
2644 ! variations of output variables: tau ut vt dvsrc dusrc
2645 ! with respect to input variables: xv yv ubi
2646  SUBROUTINE gw_drag_prof_bgnd_d(i, pcols, pver, pgwv, ngwv, kbot, ktop&
2647 & , c, u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha&
2648 & , dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, &
2649 & ubm, xv, xvd, yv, yvd, ut, utd, vt, vtd, tt, taugwx, taugwy, fegw, &
2650 & fepgw, dusrc, dusrcd, dvsrc, dvsrcd, dtsrc, tau0x, tau0y, effgw)
2651  IMPLICIT NONE
2652 !-----------------------------------------------------------------------
2653 ! Solve for the drag profile from the multiple gravity wave drag
2654 ! parameterization.
2655 ! 1. scan up from the wave source to determine the stress profile
2656 ! 2. scan down the stress profile to determine the tendencies
2657 ! => apply bounds to the tendency
2658 ! a. from wkb solution
2659 ! b. from computational stability constraints
2660 ! => adjust stress on interface below to reflect actual bounded tendency
2661 !-----------------------------------------------------------------------
2662 !------------------------------Arguments--------------------------------
2663 ! number of current column
2664  INTEGER, INTENT(IN) :: i
2665 ! number of atmospheric columns
2666  INTEGER, INTENT(IN) :: pcols
2667 ! number of atmospheric columns
2668  INTEGER, INTENT(IN) :: pver
2669 ! index of bottom (source) interface
2670  INTEGER, INTENT(IN) :: kbot
2671 ! index of top interface of gwd region
2672  INTEGER, INTENT(IN) :: ktop
2673 ! number of gravity waves possible
2674  INTEGER, INTENT(IN) :: pgwv
2675 ! number of gravity waves to use
2676  INTEGER, INTENT(IN) :: ngwv
2677 ! top interface of low level stress divergence region
2678  INTEGER, INTENT(IN) :: kldv
2679 ! min value of kldv
2680  INTEGER, INTENT(IN) :: kldvmn
2681 ! index of top interface of source region
2682  INTEGER, INTENT(IN) :: ksrc
2683 ! min value of ksrc
2684  INTEGER, INTENT(IN) :: ksrcmn
2685 ! wave phase speeds
2686  REAL :: c(-pgwv:pgwv)
2687 ! midpoint zonal wind
2688  REAL :: u(pcols, pver)
2689 ! midpoint meridional wind
2690  REAL :: v(pcols, pver)
2691 ! midpoint temperatures
2692  REAL :: t(pcols, pver)
2693 ! interface pressures
2694  REAL :: pi(pcols, 0:pver)
2695 ! midpoint delta p (pi(k)-pi(k-1))
2696  REAL :: dpm(pcols, pver)
2697 ! 1. / (pi(k)-pi(k-1))
2698  REAL :: rdpm(pcols, pver)
2699 ! ln(interface pressures)
2700  REAL :: piln(pcols, 0:pver)
2701  REAL :: rlat(pcols)
2702 ! interface density
2703  REAL :: rhoi(0:pver)
2704 ! interface Brunt-Vaisalla frequency
2705  REAL :: ni(0:pver)
2706 ! interface temperature
2707  REAL :: ti(0:pver)
2708 ! midpoint Brunt-Vaisalla frequency
2709  REAL :: nm(pver)
2710 ! time step
2711  REAL :: dt
2712 ! newtonian cooling coefficients
2713  REAL :: alpha(0:pver)
2714 ! newtonian cooling coefficients
2715  REAL :: dback(0:pver)
2716 ! 1/dp across low level divergence region
2717  REAL :: rdpldv
2718 ! projection of wind at interfaces
2719  REAL :: ubi(0:pver)
2720  REAL :: ubid(0:pver)
2721 ! projection of wind at midpoints
2722  REAL :: ubm(pver)
2723 ! unit vectors of source wind (x)
2724  REAL :: xv
2725  REAL :: xvd
2726 ! unit vectors of source wind (y)
2727  REAL :: yv
2728  REAL :: yvd
2729 ! tendency efficiency for gwd
2730  REAL :: effgw(pcols)
2731 ! wave Reynolds stress
2732  REAL :: tau(-pgwv:pgwv, 0:pver)
2733  REAL :: taud(-pgwv:pgwv, 0:pver)
2734 ! zonal wind tendency
2735  REAL :: ut(pver)
2736  REAL :: utd(pver)
2737 ! meridional wind tendency
2738  REAL :: vt(pver)
2739  REAL :: vtd(pver)
2740 ! temperature tendency
2741  REAL :: tt(pver)
2742 ! Total zonal GW momentum flux
2743  REAL :: taugwx(pcols, 0:pver)
2744 ! Total meridional GW momentum flux
2745  REAL :: taugwy(pcols, 0:pver)
2746 ! Total GW energy flux
2747  REAL :: fegw(pcols, 0:pver)
2748 ! Total GW pseudo energy flux
2749  REAL :: fepgw(pcols, 0:pver)
2750 ! Total U tendency below launch level
2751  REAL :: dusrc(pver)
2752  REAL :: dusrcd(pver)
2753 ! Total V tendency below launch level
2754  REAL :: dvsrc(pver)
2755  REAL :: dvsrcd(pver)
2756 ! Total V tendency below launch level
2757  REAL :: dtsrc(pver)
2758 ! c=0 sfc. stress (zonal)
2759  REAL :: tau0x
2760 ! c=0 sfc. stress (meridional)
2761  REAL :: tau0y
2762 !---------------------------Local storage-------------------------------
2763 ! loop indexes
2764  INTEGER :: k, l
2765 ! fraction of dsat to use
2766  REAL :: dscal
2767 ! imaginary part of vertical wavenumber
2768  REAL :: mi
2769 ! stress after damping
2770  REAL :: taudmp
2771 !real :: taumax(pcols) ! max(tau) for any l
2772 ! saturation stress
2773  REAL :: tausat
2774  REAL :: tausatd
2775 ! (ub-c)
2776  REAL :: ubmc
2777  REAL :: ubmcd
2778 ! (ub-c)**2
2779  REAL :: ubmc2
2780 ! ubar tendency
2781  REAL :: ubt
2782  REAL :: ubtd
2783 ! tbar tendency
2784  REAL :: tbt
2785 ! ubar tendency from wave l
2786  REAL :: ubtl
2787  REAL :: ubtld
2788 ! saturation tendency
2789  REAL :: ubtlsat
2790 ! layer pressure
2791  REAL :: pm
2792 ! layer density
2793  REAL :: rhom
2794 ! launch level height
2795  REAL :: zlb
2796 ! c-u
2797  REAL :: cmu
2798  REAL :: dzm, hscal, tautmp
2799  REAL :: utl
2800  REAL :: utld
2801  REAL :: ttl
2802 ! zonal pseudomomentum flux spectrum
2803  REAL :: fpmx
2804  REAL :: fpmxd
2805 ! meridional pseudomomentum flux spectrum
2806  REAL :: fpmy
2807  REAL :: fpmyd
2808 ! energy flux (p'w') spectrum
2809  REAL :: fe
2810 ! pseudoenergy flux (p'w'+U rho u'w') spectrum
2811  REAL :: fpe
2812  REAL :: fpml
2813  REAL :: fpmld
2814  REAL :: fpmt
2815  REAL :: fpmtd
2816  REAL :: fpel
2817  REAL :: fpet
2818  REAL :: dusrcl
2819  REAL :: dusrcld
2820  REAL :: dvsrcl
2821  REAL :: dvsrcld
2822  REAL :: dtsrcl
2823  REAL :: zi
2824  REAL :: effkwvmap
2825  REAL :: zfac
2826  REAL :: uhtmax
2827  REAL :: uhtmaxd
2828  REAL :: utfac
2829  REAL :: utfacd
2830  REAL :: tau_min, cmu_, t_new
2831  REAL :: tau_mind
2832  REAL :: arg1
2833  REAL :: arg1d
2834  INTRINSIC max
2835  INTRINSIC sign
2836  INTRINSIC abs
2837  REAL :: x2
2838  REAL :: x2d
2839  REAL :: x1
2840  REAL :: x1d
2841  INTRINSIC sqrt
2842 ! Initialize gravity wave drag tendencies to zero
2843  DO l=-ngwv,ngwv
2844  DO k=pver-1,ktop,-1
2845  IF (k .LE. kbot - 1) THEN
2846  taud(l, k) = 0.0_8
2847  tau(l, k) = 0.
2848  END IF
2849  END DO
2850  END DO
2851  DO k=1,pver
2852  utd(k) = 0.0_8
2853  ut(k) = 0.
2854  vtd(k) = 0.0_8
2855  vt(k) = 0.
2856  tt(k) = 0.
2857  dusrcd(k) = 0.0_8
2858  dusrc(k) = 0.
2859  dvsrcd(k) = 0.0_8
2860  dvsrc(k) = 0.
2861  dtsrc(k) = 0.
2862  END DO
2863 ! Initialize total momentum and energy fluxes to zero
2864  DO k=0,pver
2865  taugwx(i, k) = 0.
2866  taugwy(i, k) = 0.
2867  fegw(i, k) = 0.
2868  fepgw(i, k) = 0.
2869  END DO
2870 ! Initialize surface wave stress at c = 0 to zero
2871  tau0x = 0.
2872  tau0y = 0.
2873  taud(:, :) = 0.0_8
2874 !---------------------------------------------------------------------------
2875 ! Compute the stress profiles and diffusivities
2876 !---------------------------------------------------------------------------
2877 ! Determine the absolute value of the saturation stress and the diffusivity
2878 ! for each wave.
2879 ! Define critical levels where the sign of (u-c) changes between interfaces.
2880 ! Loop from bottom to top to get stress profiles
2881  DO l=-ngwv,ngwv
2882  DO k=pver-1,ktop,-1
2883  IF (k .LE. kbot - 1) THEN
2884 !!!!!!!!jk d = dback(k)
2885  ubmcd = ubid(k)
2886  ubmc = ubi(k) - c(l)
2887  IF (ngwv .GT. 0) THEN
2888  CALL get_effkwvmap_1(effkwvmap, rlat(i))
2889  IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
2890 & .LT. 15.0) effkwvmap = fcrit2*kwvbeq
2891  ELSE
2892  IF (pi(i, k) .LT. 1000.0) THEN
2893  zfac = (pi(i, k)/1000.0)**3
2894  ELSE
2895  zfac = 1.0
2896  END IF
2897  IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
2898  CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
2899  ELSE
2900  CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
2901  END IF
2902  END IF
2903  x1d = effkwvmap*rhoi(k)*3*ubmc**2*ubmcd/(2.*ni(k))
2904  x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
2905  IF (x1 .GE. 0.) THEN
2906  tausatd = x1d
2907  tausat = x1
2908  ELSE
2909  tausatd = -x1d
2910  tausat = -x1
2911  END IF
2912  IF (tausat .LE. taumin) THEN
2913  tausat = 0.0
2914  tausatd = 0.0_8
2915  END IF
2916  IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) THEN
2917  tausat = 0.0
2918  tausatd = 0.0_8
2919  END IF
2920  IF (k .EQ. ktop) THEN
2921  tausat = 0.
2922  tausatd = 0.0_8
2923  END IF
2924 !
2925  IF (k .EQ. ktop + 1) THEN
2926  tausatd = 0.02*tausatd
2927  tausat = tausat*0.02
2928  END IF
2929  IF (k .EQ. ktop + 2) THEN
2930  tausatd = 0.05*tausatd
2931  tausat = tausat*0.05
2932  END IF
2933  IF (k .EQ. ktop + 3) THEN
2934  tausatd = 0.10*tausatd
2935  tausat = tausat*0.10
2936  END IF
2937  IF (k .EQ. ktop + 4) THEN
2938  tausatd = 0.20*tausatd
2939  tausat = tausat*0.20
2940  END IF
2941  IF (k .EQ. ktop + 5) THEN
2942  tausatd = 0.50*tausatd
2943  tausat = tausat*0.50
2944  END IF
2945 !
2946  tau_mind = tausatd
2947  tau_min = tausat
2948  IF (tau_min .GT. tau(l, k+1)) THEN
2949  tau_mind = taud(l, k+1)
2950  tau_min = tau(l, k+1)
2951  END IF
2952  taud(l, k) = tau_mind
2953  tau(l, k) = tau_min
2954  END IF
2955  END DO
2956  END DO
2957  utd(:) = 0.0_8
2958  vtd(:) = 0.0_8
2959 !---------------------------------------------------------------------------
2960 ! Compute the tendencies from the stress divergence.
2961 !---------------------------------------------------------------------------
2962 ! Accumulate the mean wind tendency over wavenumber.
2963 ! Loop over levels from top to bottom
2964  DO k=ktop+1,pver
2965  ubt = 0.0
2966  tbt = 0.0
2967  ubtd = 0.0_8
2968  DO l=-ngwv,ngwv
2969  IF (k .LE. kbot) THEN
2970 ! Determine the wind tendency including excess stress carried down from above.
2971  ubtld = mapl_grav*rdpm(i, k)*(taud(l, k)-taud(l, k-1))
2972  ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
2973 ! Calculate the sign of wind tendency
2974  utld = ubtld*sign(1.d0, ubtl*(c(l)-ubi(k)))
2975  utl = sign(ubtl, c(l) - ubi(k))
2976 ! Accumulate the mean wind tendency over wavenumber.
2977  ubtd = ubtd + utld
2978  ubt = ubt + utl
2979 ! Calculate irreversible temperature tendency associated with gravity wave breaking.
2980  ttl = (c(l)-ubm(k))*utl/mapl_cp
2981 ! Adding frictional heating associated with the GW momentum forcing
2982  tbt = tbt + ttl
2983  END IF
2984  END DO
2985 ! Project the mean wind tendency onto the components and scale by "efficiency".
2986  IF (k .LE. kbot) THEN
2987  utd(k) = effgw(i)*(ubtd*xv+ubt*xvd)
2988  ut(k) = ubt*xv*effgw(i)
2989  vtd(k) = effgw(i)*(ubtd*yv+ubt*yvd)
2990  vt(k) = ubt*yv*effgw(i)
2991  tt(k) = tbt*effgw(i)
2992  END IF
2993  END DO
2994 !-----------------------------------------------------------------------
2995 ! Calculates wind and temperature tendencies below launch level for
2996 ! energy and momentum conservation (does nothing for orographic GWs).
2997 !-----------------------------------------------------------------------
2998 !----kimmmmm here
2999 ! Calculate launch level height
3000  zlb = 0.
3001  DO k=ktop+1,pver
3002  IF (k .GE. kbot + 1) THEN
3003 ! Define layer pressure and density
3004  pm = (pi(i, k-1)+pi(i, k))*0.5
3005  CALL get_ti(t_new, t(i, k))
3006  rhom = pm/(mapl_rgas*t_new)
3007  zlb = zlb + dpm(i, k)/mapl_grav/rhom
3008  END IF
3009  END DO
3010  dvsrcd(:) = 0.0_8
3011  dusrcd(:) = 0.0_8
3012  fpmld = 0.0_8
3013  fpmtd = 0.0_8
3014 !-----------------------------------------------------------------------
3015 ! Calculates energy and momentum flux profiles
3016 !-----------------------------------------------------------------------
3017  DO l=-ngwv,ngwv
3018  DO k=ktop,pver
3019  IF (k .LE. kbot) THEN
3020  cmu = c(l) - ubi(k)
3021  CALL get_cmu(cmu_, cmu)
3022  fpmxd = cmu_*effgw(i)*(taud(l, k)*xv+tau(l, k)*xvd)
3023  fpmx = cmu_*tau(l, k)*xv*effgw(i)
3024  fpmyd = cmu_*effgw(i)*(taud(l, k)*yv+tau(l, k)*yvd)
3025  fpmy = cmu_*tau(l, k)*yv*effgw(i)
3026  fe = cmu*cmu_*tau(l, k)*effgw(i)
3027  fpe = c(l)*cmu_*tau(l, k)*effgw(i)
3028  IF (k .EQ. kbot) THEN
3029  fpmld = fpmxd*xv + fpmx*xvd + fpmyd*yv + fpmy*yvd
3030  fpml = fpmx*xv + fpmy*yv
3031  fpel = fpe
3032  END IF
3033  IF (k .EQ. ktop) THEN
3034  fpmtd = fpmxd*xv + fpmx*xvd + fpmyd*yv + fpmy*yvd
3035  fpmt = fpmx*xv + fpmy*yv
3036  fpet = fpe
3037  END IF
3038 ! Record outputs for GW fluxes
3039  taugwx(i, k) = taugwx(i, k) + fpmx
3040  taugwy(i, k) = taugwy(i, k) + fpmy
3041  fegw(i, k) = fegw(i, k) + fe
3042  fepgw(i, k) = fepgw(i, k) + fpe
3043  END IF
3044  END DO
3045  DO k=ktop+1,pver
3046  IF (k .GE. kbot + 1) THEN
3047 ! Define layer pressure and density
3048  pm = (pi(i, k-1)+pi(i, k))*0.5
3049  CALL get_ti(t_new, t(i, k))
3050  rhom = pm/(mapl_rgas*t_new)
3051  dusrcld = -((fpmld-fpmtd)*xv/(rhom*zlb)+(fpml-fpmt)*xvd/(rhom*&
3052 & zlb))
3053  dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
3054  dvsrcld = -((fpmld-fpmtd)*yv/(rhom*zlb)+(fpml-fpmt)*yvd/(rhom*&
3055 & zlb))
3056  dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
3057  dtsrcl = -((fpel-fpet-ubm(k)*(fpml-fpmt))/(rhom*zlb*mapl_cp))
3058 ! Add sub-source wind and temperature tendencies
3059  dusrcd(k) = dusrcd(k) + dusrcld
3060  dusrc(k) = dusrc(k) + dusrcl
3061  dvsrcd(k) = dvsrcd(k) + dvsrcld
3062  dvsrc(k) = dvsrc(k) + dvsrcl
3063  dtsrc(k) = dtsrc(k) + dtsrcl
3064  END IF
3065  END DO
3066  END DO
3067 ! For orographic waves, sub-source tendencies are set equal to zero.
3068  DO k=1,pver
3069  IF (ngwv .EQ. 0) THEN
3070  dusrcd(k) = 0.0_8
3071  dusrc(k) = 0.0
3072  dvsrcd(k) = 0.0_8
3073  dvsrc(k) = 0.0
3074  dtsrc(k) = 0.0
3075  END IF
3076  END DO
3077 !-----------------------------------------------------------------------
3078 ! Adjust efficiency factor to prevent unrealistically strong forcing
3079 !-----------------------------------------------------------------------
3080  uhtmax = 0.0
3081  utfac = 1.0
3082  uhtmaxd = 0.0_8
3083  DO k=1,pver
3084  arg1d = 2*ut(k)*utd(k) + 2*vt(k)*vtd(k)
3085  arg1 = ut(k)**2 + vt(k)**2
3086  IF (arg1 .EQ. 0.0) THEN
3087  x2d = 0.0_8
3088  ELSE
3089  x2d = arg1d/(2.0*sqrt(arg1))
3090  END IF
3091  x2 = sqrt(arg1)
3092  IF (x2 .LT. uhtmax) THEN
3093  uhtmax = uhtmax
3094  ELSE
3095  uhtmaxd = x2d
3096  uhtmax = x2
3097  END IF
3098  END DO
3099  IF (uhtmax .GT. tndmax) THEN
3100  utfacd = -(tndmax*uhtmaxd/uhtmax**2)
3101  utfac = tndmax/uhtmax
3102  ELSE
3103  utfacd = 0.0_8
3104  END IF
3105 !jkim effgw(i) = effgw(i)*utfac
3106  DO k=1,pver
3107  utd(k) = utd(k)*utfac + ut(k)*utfacd
3108  ut(k) = ut(k)*utfac
3109  vtd(k) = vtd(k)*utfac + vt(k)*utfacd
3110  vt(k) = vt(k)*utfac
3111  tt(k) = tt(k)*utfac
3112  dusrcd(k) = dusrcd(k)*utfac + dusrc(k)*utfacd
3113  dusrc(k) = dusrc(k)*utfac
3114  dvsrcd(k) = dvsrcd(k)*utfac + dvsrc(k)*utfacd
3115  dvsrc(k) = dvsrc(k)*utfac
3116  dtsrc(k) = dtsrc(k)*utfac
3117  END DO
3118 !-----------------------------------------------------------------------
3119 ! Project the c=0 stress (scaled) in the direction of the source wind
3120 ! for recording on the output file.
3121 !-----------------------------------------------------------------------
3122  tau0x = tau(0, kbot)*xv*effgw(i)*utfac
3123  tau0y = tau(0, kbot)*yv*effgw(i)*utfac
3124  RETURN
3125  END SUBROUTINE gw_drag_prof_bgnd_d
3126  SUBROUTINE gw_drag_prof_bgnd(i, pcols, pver, pgwv, ngwv, kbot, ktop, c&
3127 & , u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, &
3128 & dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut&
3129 & , vt, tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, &
3130 & tau0y, effgw)
3131  IMPLICIT NONE
3132 !-----------------------------------------------------------------------
3133 ! Solve for the drag profile from the multiple gravity wave drag
3134 ! parameterization.
3135 ! 1. scan up from the wave source to determine the stress profile
3136 ! 2. scan down the stress profile to determine the tendencies
3137 ! => apply bounds to the tendency
3138 ! a. from wkb solution
3139 ! b. from computational stability constraints
3140 ! => adjust stress on interface below to reflect actual bounded tendency
3141 !-----------------------------------------------------------------------
3142 !------------------------------Arguments--------------------------------
3143 ! number of current column
3144  INTEGER, INTENT(IN) :: i
3145 ! number of atmospheric columns
3146  INTEGER, INTENT(IN) :: pcols
3147 ! number of atmospheric columns
3148  INTEGER, INTENT(IN) :: pver
3149 ! index of bottom (source) interface
3150  INTEGER, INTENT(IN) :: kbot
3151 ! index of top interface of gwd region
3152  INTEGER, INTENT(IN) :: ktop
3153 ! number of gravity waves possible
3154  INTEGER, INTENT(IN) :: pgwv
3155 ! number of gravity waves to use
3156  INTEGER, INTENT(IN) :: ngwv
3157 ! top interface of low level stress divergence region
3158  INTEGER, INTENT(IN) :: kldv
3159 ! min value of kldv
3160  INTEGER, INTENT(IN) :: kldvmn
3161 ! index of top interface of source region
3162  INTEGER, INTENT(IN) :: ksrc
3163 ! min value of ksrc
3164  INTEGER, INTENT(IN) :: ksrcmn
3165 ! wave phase speeds
3166  REAL :: c(-pgwv:pgwv)
3167 ! midpoint zonal wind
3168  REAL :: u(pcols, pver)
3169 ! midpoint meridional wind
3170  REAL :: v(pcols, pver)
3171 ! midpoint temperatures
3172  REAL :: t(pcols, pver)
3173 ! interface pressures
3174  REAL :: pi(pcols, 0:pver)
3175 ! midpoint delta p (pi(k)-pi(k-1))
3176  REAL :: dpm(pcols, pver)
3177 ! 1. / (pi(k)-pi(k-1))
3178  REAL :: rdpm(pcols, pver)
3179 ! ln(interface pressures)
3180  REAL :: piln(pcols, 0:pver)
3181  REAL :: rlat(pcols)
3182 ! interface density
3183  REAL :: rhoi(0:pver)
3184 ! interface Brunt-Vaisalla frequency
3185  REAL :: ni(0:pver)
3186 ! interface temperature
3187  REAL :: ti(0:pver)
3188 ! midpoint Brunt-Vaisalla frequency
3189  REAL :: nm(pver)
3190 ! time step
3191  REAL :: dt
3192 ! newtonian cooling coefficients
3193  REAL :: alpha(0:pver)
3194 ! newtonian cooling coefficients
3195  REAL :: dback(0:pver)
3196 ! 1/dp across low level divergence region
3197  REAL :: rdpldv
3198 ! projection of wind at interfaces
3199  REAL :: ubi(0:pver)
3200 ! projection of wind at midpoints
3201  REAL :: ubm(pver)
3202 ! unit vectors of source wind (x)
3203  REAL :: xv
3204 ! unit vectors of source wind (y)
3205  REAL :: yv
3206 ! tendency efficiency for gwd
3207  REAL :: effgw(pcols)
3208 ! wave Reynolds stress
3209  REAL :: tau(-pgwv:pgwv, 0:pver)
3210 ! zonal wind tendency
3211  REAL :: ut(pver)
3212 ! meridional wind tendency
3213  REAL :: vt(pver)
3214 ! temperature tendency
3215  REAL :: tt(pver)
3216 ! Total zonal GW momentum flux
3217  REAL :: taugwx(pcols, 0:pver)
3218 ! Total meridional GW momentum flux
3219  REAL :: taugwy(pcols, 0:pver)
3220 ! Total GW energy flux
3221  REAL :: fegw(pcols, 0:pver)
3222 ! Total GW pseudo energy flux
3223  REAL :: fepgw(pcols, 0:pver)
3224 ! Total U tendency below launch level
3225  REAL :: dusrc(pver)
3226 ! Total V tendency below launch level
3227  REAL :: dvsrc(pver)
3228 ! Total V tendency below launch level
3229  REAL :: dtsrc(pver)
3230 ! c=0 sfc. stress (zonal)
3231  REAL :: tau0x
3232 ! c=0 sfc. stress (meridional)
3233  REAL :: tau0y
3234 !---------------------------Local storage-------------------------------
3235 ! loop indexes
3236  INTEGER :: k, l
3237 ! fraction of dsat to use
3238  REAL :: dscal
3239 ! imaginary part of vertical wavenumber
3240  REAL :: mi
3241 ! stress after damping
3242  REAL :: taudmp
3243 !real :: taumax(pcols) ! max(tau) for any l
3244 ! saturation stress
3245  REAL :: tausat
3246 ! (ub-c)
3247  REAL :: ubmc
3248 ! (ub-c)**2
3249  REAL :: ubmc2
3250 ! ubar tendency
3251  REAL :: ubt
3252 ! tbar tendency
3253  REAL :: tbt
3254 ! ubar tendency from wave l
3255  REAL :: ubtl
3256 ! saturation tendency
3257  REAL :: ubtlsat
3258 ! layer pressure
3259  REAL :: pm
3260 ! layer density
3261  REAL :: rhom
3262 ! launch level height
3263  REAL :: zlb
3264 ! c-u
3265  REAL :: cmu
3266  REAL :: dzm, hscal, tautmp
3267  REAL :: utl
3268  REAL :: ttl
3269 ! zonal pseudomomentum flux spectrum
3270  REAL :: fpmx
3271 ! meridional pseudomomentum flux spectrum
3272  REAL :: fpmy
3273 ! energy flux (p'w') spectrum
3274  REAL :: fe
3275 ! pseudoenergy flux (p'w'+U rho u'w') spectrum
3276  REAL :: fpe
3277  REAL :: fpml
3278  REAL :: fpmt
3279  REAL :: fpel
3280  REAL :: fpet
3281  REAL :: dusrcl
3282  REAL :: dvsrcl
3283  REAL :: dtsrcl
3284  REAL :: zi
3285  REAL :: effkwvmap
3286  REAL :: zfac
3287  REAL :: uhtmax
3288  REAL :: utfac
3289  REAL :: tau_min, cmu_, t_new
3290  REAL :: arg1
3291  INTRINSIC max
3292  INTRINSIC sign
3293  INTRINSIC abs
3294  REAL :: x2
3295  REAL :: x1
3296  INTRINSIC sqrt
3297 ! Initialize gravity wave drag tendencies to zero
3298  DO l=-ngwv,ngwv
3299  DO k=pver-1,ktop,-1
3300  IF (k .LE. kbot - 1) tau(l, k) = 0.
3301  END DO
3302  END DO
3303  DO k=1,pver
3304  ut(k) = 0.
3305  vt(k) = 0.
3306  tt(k) = 0.
3307  dusrc(k) = 0.
3308  dvsrc(k) = 0.
3309  dtsrc(k) = 0.
3310  END DO
3311 ! Initialize total momentum and energy fluxes to zero
3312  DO k=0,pver
3313  taugwx(i, k) = 0.
3314  taugwy(i, k) = 0.
3315  fegw(i, k) = 0.
3316  fepgw(i, k) = 0.
3317  END DO
3318 ! Initialize surface wave stress at c = 0 to zero
3319  tau0x = 0.
3320  tau0y = 0.
3321 !---------------------------------------------------------------------------
3322 ! Compute the stress profiles and diffusivities
3323 !---------------------------------------------------------------------------
3324 ! Determine the absolute value of the saturation stress and the diffusivity
3325 ! for each wave.
3326 ! Define critical levels where the sign of (u-c) changes between interfaces.
3327 ! Loop from bottom to top to get stress profiles
3328  DO l=-ngwv,ngwv
3329  DO k=pver-1,ktop,-1
3330  IF (k .LE. kbot - 1) THEN
3331 !!!!!!!!jk d = dback(k)
3332  ubmc = ubi(k) - c(l)
3333  IF (ngwv .GT. 0) THEN
3334  CALL get_effkwvmap_1(effkwvmap, rlat(i))
3335  IF (-15.0 .LT. rlat(i)*180./pi_gwd .AND. rlat(i)*180./pi_gwd&
3336 & .LT. 15.0) effkwvmap = fcrit2*kwvbeq
3337  ELSE
3338  IF (pi(i, k) .LT. 1000.0) THEN
3339  zfac = (pi(i, k)/1000.0)**3
3340  ELSE
3341  zfac = 1.0
3342  END IF
3343  IF (rlat(i)*180./pi_gwd .LT. -20.0) THEN
3344  CALL get_effkwvmap_2(effkwvmap, rlat(i), zfac)
3345  ELSE
3346  CALL get_effkwvmap_3(effkwvmap, rlat(i), zfac)
3347  END IF
3348  END IF
3349  x1 = effkwvmap*rhoi(k)*ubmc**3/(2.*ni(k))
3350  IF (x1 .GE. 0.) THEN
3351  tausat = x1
3352  ELSE
3353  tausat = -x1
3354  END IF
3355  IF (tausat .LE. taumin) tausat = 0.0
3356  IF (ubmc*(ubi(k+1)-c(l)) .LE. 0.0) tausat = 0.0
3357  IF (k .EQ. ktop) tausat = 0.
3358 !
3359  IF (k .EQ. ktop + 1) tausat = tausat*0.02
3360  IF (k .EQ. ktop + 2) tausat = tausat*0.05
3361  IF (k .EQ. ktop + 3) tausat = tausat*0.10
3362  IF (k .EQ. ktop + 4) tausat = tausat*0.20
3363  IF (k .EQ. ktop + 5) tausat = tausat*0.50
3364 !
3365  tau_min = tausat
3366  IF (tau_min .GT. tau(l, k+1)) tau_min = tau(l, k+1)
3367  tau(l, k) = tau_min
3368  END IF
3369  END DO
3370  END DO
3371 !---------------------------------------------------------------------------
3372 ! Compute the tendencies from the stress divergence.
3373 !---------------------------------------------------------------------------
3374 ! Accumulate the mean wind tendency over wavenumber.
3375 ! Loop over levels from top to bottom
3376  DO k=ktop+1,pver
3377  ubt = 0.0
3378  tbt = 0.0
3379  DO l=-ngwv,ngwv
3380  IF (k .LE. kbot) THEN
3381 ! Determine the wind tendency including excess stress carried down from above.
3382  ubtl = mapl_grav*(tau(l, k)-tau(l, k-1))*rdpm(i, k)
3383 ! Calculate the sign of wind tendency
3384  utl = sign(ubtl, c(l) - ubi(k))
3385 ! Accumulate the mean wind tendency over wavenumber.
3386  ubt = ubt + utl
3387 ! Calculate irreversible temperature tendency associated with gravity wave breaking.
3388  ttl = (c(l)-ubm(k))*utl/mapl_cp
3389 ! Adding frictional heating associated with the GW momentum forcing
3390  tbt = tbt + ttl
3391  END IF
3392  END DO
3393 ! Project the mean wind tendency onto the components and scale by "efficiency".
3394  IF (k .LE. kbot) THEN
3395  ut(k) = ubt*xv*effgw(i)
3396  vt(k) = ubt*yv*effgw(i)
3397  tt(k) = tbt*effgw(i)
3398  END IF
3399  END DO
3400 !-----------------------------------------------------------------------
3401 ! Calculates wind and temperature tendencies below launch level for
3402 ! energy and momentum conservation (does nothing for orographic GWs).
3403 !-----------------------------------------------------------------------
3404 !----kimmmmm here
3405 ! Calculate launch level height
3406  zlb = 0.
3407  DO k=ktop+1,pver
3408  IF (k .GE. kbot + 1) THEN
3409 ! Define layer pressure and density
3410  pm = (pi(i, k-1)+pi(i, k))*0.5
3411  CALL get_ti(t_new, t(i, k))
3412  rhom = pm/(mapl_rgas*t_new)
3413  zlb = zlb + dpm(i, k)/mapl_grav/rhom
3414  END IF
3415  END DO
3416 !-----------------------------------------------------------------------
3417 ! Calculates energy and momentum flux profiles
3418 !-----------------------------------------------------------------------
3419  DO l=-ngwv,ngwv
3420  DO k=ktop,pver
3421  IF (k .LE. kbot) THEN
3422  cmu = c(l) - ubi(k)
3423  CALL get_cmu(cmu_, cmu)
3424  fpmx = cmu_*tau(l, k)*xv*effgw(i)
3425  fpmy = cmu_*tau(l, k)*yv*effgw(i)
3426  fe = cmu*cmu_*tau(l, k)*effgw(i)
3427  fpe = c(l)*cmu_*tau(l, k)*effgw(i)
3428  IF (k .EQ. kbot) THEN
3429  fpml = fpmx*xv + fpmy*yv
3430  fpel = fpe
3431  END IF
3432  IF (k .EQ. ktop) THEN
3433  fpmt = fpmx*xv + fpmy*yv
3434  fpet = fpe
3435  END IF
3436 ! Record outputs for GW fluxes
3437  taugwx(i, k) = taugwx(i, k) + fpmx
3438  taugwy(i, k) = taugwy(i, k) + fpmy
3439  fegw(i, k) = fegw(i, k) + fe
3440  fepgw(i, k) = fepgw(i, k) + fpe
3441  END IF
3442  END DO
3443  DO k=ktop+1,pver
3444  IF (k .GE. kbot + 1) THEN
3445 ! Define layer pressure and density
3446  pm = (pi(i, k-1)+pi(i, k))*0.5
3447  CALL get_ti(t_new, t(i, k))
3448  rhom = pm/(mapl_rgas*t_new)
3449  dusrcl = -((fpml-fpmt)/(rhom*zlb)*xv)
3450  dvsrcl = -((fpml-fpmt)/(rhom*zlb)*yv)
3451  dtsrcl = -((fpel-fpet-ubm(k)*(fpml-fpmt))/(rhom*zlb*mapl_cp))
3452 ! Add sub-source wind and temperature tendencies
3453  dusrc(k) = dusrc(k) + dusrcl
3454  dvsrc(k) = dvsrc(k) + dvsrcl
3455  dtsrc(k) = dtsrc(k) + dtsrcl
3456  END IF
3457  END DO
3458  END DO
3459 ! For orographic waves, sub-source tendencies are set equal to zero.
3460  DO k=1,pver
3461  IF (ngwv .EQ. 0) THEN
3462  dusrc(k) = 0.0
3463  dvsrc(k) = 0.0
3464  dtsrc(k) = 0.0
3465  END IF
3466  END DO
3467 !-----------------------------------------------------------------------
3468 ! Adjust efficiency factor to prevent unrealistically strong forcing
3469 !-----------------------------------------------------------------------
3470  uhtmax = 0.0
3471  utfac = 1.0
3472  DO k=1,pver
3473  arg1 = ut(k)**2 + vt(k)**2
3474  x2 = sqrt(arg1)
3475  IF (x2 .LT. uhtmax) THEN
3476  uhtmax = uhtmax
3477  ELSE
3478  uhtmax = x2
3479  END IF
3480  END DO
3481  IF (uhtmax .GT. tndmax) utfac = tndmax/uhtmax
3482 !jkim effgw(i) = effgw(i)*utfac
3483  DO k=1,pver
3484  ut(k) = ut(k)*utfac
3485  vt(k) = vt(k)*utfac
3486  tt(k) = tt(k)*utfac
3487  dusrc(k) = dusrc(k)*utfac
3488  dvsrc(k) = dvsrc(k)*utfac
3489  dtsrc(k) = dtsrc(k)*utfac
3490  END DO
3491 !-----------------------------------------------------------------------
3492 ! Project the c=0 stress (scaled) in the direction of the source wind
3493 ! for recording on the output file.
3494 !-----------------------------------------------------------------------
3495  tau0x = tau(0, kbot)*xv*effgw(i)*utfac
3496  tau0y = tau(0, kbot)*yv*effgw(i)*utfac
3497  RETURN
3498  END SUBROUTINE gw_drag_prof_bgnd
3499 ! Differentiation of get_uv in forward (tangent) mode:
3500 ! variations of output variables: uv_out
3501 ! with respect to input variables: uv_in
3502  SUBROUTINE get_uv_d(uv_out, uv_outd, uv_in, uv_ind)
3503  IMPLICIT NONE
3504  REAL :: uv_out, uv_in
3505  REAL :: uv_outd, uv_ind
3506  LOGICAL :: src_bypass
3507  IF (src_bypass) THEN
3508  uv_outd = uv_ind
3509  uv_out = uv_in
3510  ELSE
3511  uv_outd = uv_ind
3512  uv_out = uv_in
3513  END IF
3514  END SUBROUTINE get_uv_d
3515  SUBROUTINE get_uv(uv_out, uv_in)
3516  IMPLICIT NONE
3517  REAL :: uv_out, uv_in
3518  LOGICAL :: src_bypass
3519  IF (src_bypass) THEN
3520  uv_out = uv_in
3521  ELSE
3522  uv_out = uv_in
3523  END IF
3524  END SUBROUTINE get_uv
3525  SUBROUTINE get_effkwvmap_1(effkwvmap_, rlat_)
3526  IMPLICIT NONE
3527  REAL :: tanh1, tanh2, effkwvmap_, rlat_
3528  INTRINSIC tanh
3529  tanh1 = (rlat_*180./pi_gwd-20.0)/6.0
3530  tanh2 = -((rlat_*180./pi_gwd+20.0)/6.0)
3531  tanh1 = tanh(tanh1)
3532  tanh2 = tanh(tanh2)
3533  effkwvmap_ = fcrit2*kwvb*(0.5*(1.0+tanh1)+0.5*(1.0+tanh2))
3534  END SUBROUTINE get_effkwvmap_1
3535  SUBROUTINE get_effkwvmap_2(effkwvmap_, rlat_, zfac)
3536  IMPLICIT NONE
3537  REAL :: effkwvmap_, rlat_, zfac
3538  INTRINSIC exp
3539  INTRINSIC abs
3540  REAL :: abs1
3541  IF (rlat_*180./pi_gwd + 20. .GE. 0.) THEN
3542  abs1 = rlat_*180./pi_gwd + 20.
3543  ELSE
3544  abs1 = -(rlat_*180./pi_gwd+20.)
3545  END IF
3546  effkwvmap_ = fcrit2*kwvo*zfac*(2.0-0.65*exp(-(abs1/5.)))
3547  END SUBROUTINE get_effkwvmap_2
3548  SUBROUTINE get_effkwvmap_3(effkwvmap_, rlat_, zfac)
3549  IMPLICIT NONE
3550  REAL :: effkwvmap_, rlat_, zfac
3551  INTRINSIC exp
3552  INTRINSIC abs
3553  REAL :: abs1
3554  IF (rlat_*180./pi_gwd + 20. .GE. 0.) THEN
3555  abs1 = rlat_*180./pi_gwd + 20.
3556  ELSE
3557  abs1 = -(rlat_*180./pi_gwd+20.)
3558  END IF
3559  effkwvmap_ = fcrit2*kwvo*zfac*(0.7+0.65*exp(-(abs1/5.)))
3560  END SUBROUTINE get_effkwvmap_3
3561  SUBROUTINE get_cmu(cmu_, cmu)
3562  IMPLICIT NONE
3563  REAL :: cmu_, cmu
3564  INTRINSIC sign
3565  cmu_ = sign(1.0, cmu)
3566  END SUBROUTINE get_cmu
3567  SUBROUTINE get_ti(t_out, t_in)
3568  IMPLICIT NONE
3569  REAL :: t_out, t_in
3570  LOGICAL :: src_bypass
3571  t_out = t_in
3572  END SUBROUTINE get_ti
3573 END MODULE gw_drag_d
gw_drag_d::get_effkwvmap_1
subroutine get_effkwvmap_1(effkwvmap_, rlat_)
Definition: gw_drag_d.F90:3526
gw_drag_d::gw_intr
subroutine gw_intr(i, pcols, pver, dt, pgwv, effgworo_dev, effgwbkg_dev, pint_dev, t_dev, u_dev, v_dev, sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, rlat_dev, dudt_gwd_dev, dvdt_gwd_dev, dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev, dtdt_org_dev, taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
Definition: gw_drag_d.F90:590
gw_drag_init::isize_s1_tau
integer, public isize_s1_tau
Definition: gw_drag_init.F90:16
gw_drag_d::fracldv
real, parameter, public fracldv
Definition: gw_drag_d.F90:45
gw_drag_init::do_oro
logical, public do_oro
Definition: gw_drag_init.F90:28
gw_drag_d::gw_bgnd_d
subroutine gw_bgnd_d(i, pcols, pver, c, u, ud, v, vd, t, pm, pi, dpm, rdpm, piln, rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubid, ubm, ubmd, xv, xvd, yv, yvd, ngwv, kbot)
Definition: gw_drag_d.F90:1373
gw_drag_d::gw_drag_prof_bgnd_d
subroutine gw_drag_prof_bgnd_d(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, ubm, xv, xvd, yv, yvd, ut, utd, vt, vtd, tt, taugwx, taugwy, fegw, fepgw, dusrc, dusrcd, dvsrc, dvsrcd, dtsrc, tau0x, tau0y, effgw)
Definition: gw_drag_d.F90:2651
gw_drag_d::gw_drag_prof_bgnd
subroutine gw_drag_prof_bgnd(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut, vt, tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y, effgw)
Definition: gw_drag_d.F90:3131
gw_drag_init::p_src
real *8, parameter, public p_src
Definition: gw_drag_init.F90:23
gw_drag_d::mxrange
real, parameter, public mxrange
Definition: gw_drag_d.F90:49
gw_drag_d::gw_drag_prof
subroutine gw_drag_prof(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ut, vt, tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y, effgw)
Definition: gw_drag_d.F90:2275
gw_drag_d::taubgnd
real, parameter, public taubgnd
Definition: gw_drag_d.F90:59
gw_drag_d::gw_prof
subroutine gw_prof(i, k, pcols, pver, u, v, t, pm, pi, rhoi, ni, ti, nm)
Definition: gw_drag_d.F90:821
gw_drag_d::umcfac
real, parameter, public umcfac
Definition: gw_drag_d.F90:67
mapl_constantsmod::mapl_cp
real, parameter, public mapl_cp
Definition: MAPL_Constants.F90:46
mapl_constantsmod::mapl_vireps
real, parameter, public mapl_vireps
Definition: MAPL_Constants.F90:47
gw_drag_d::get_cmu
subroutine get_cmu(cmu_, cmu)
Definition: gw_drag_d.F90:3562
gw_drag_d::get_uv_d
subroutine get_uv_d(uv_out, uv_outd, uv_in, uv_ind)
Definition: gw_drag_d.F90:3503
gw_drag_init::isize_e1_tau
integer, public isize_e1_tau
Definition: gw_drag_init.F90:16
gw_drag_d::taumin
real, parameter, public taumin
Definition: gw_drag_d.F90:61
gw_drag_d::kwvbeq
real, parameter, public kwvbeq
Definition: gw_drag_d.F90:41
gw_drag_d::gw_bgnd
subroutine gw_bgnd(i, pcols, pver, c, u, v, t, pm, pi, dpm, rdpm, piln, rlat, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, ngwv, kbot)
Definition: gw_drag_d.F90:1612
gw_drag_d::kwvb
real, parameter, public kwvb
Definition: gw_drag_d.F90:39
mapl_constantsmod::mapl_kappa
real, parameter, public mapl_kappa
Definition: MAPL_Constants.F90:38
gw_drag_d::gw_oro
subroutine gw_oro(i, pcols, pver, pgwv, u, v, t, sgh, pm, pi, dpm, zm, nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, ubi, ubm, xv, yv, kbot, rlat)
Definition: gw_drag_d.F90:1162
gw_drag_d::pi_gwd
real, parameter, public pi_gwd
Definition: gw_drag_d.F90:76
mapl_constantsmod::mapl_p00
real, parameter, public mapl_p00
Definition: MAPL_Constants.F90:49
gw_drag_init::isize_s2_tau
integer, public isize_s2_tau
Definition: gw_drag_init.F90:16
gw_drag_d::tauscal
real, parameter, public tauscal
Definition: gw_drag_d.F90:63
gw_drag_init::do_bgnd
logical, public do_bgnd
Definition: gw_drag_init.F90:27
gw_drag_d::oroko2
real, parameter, public oroko2
Definition: gw_drag_d.F90:74
gw_drag_d::gw_oro_d
subroutine gw_oro_d(i, pcols, pver, pgwv, u, ud, v, vd, t, sgh, pm, pi, dpm, zm, nm, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, ubm, ubmd, xv, xvd, yv, yvd, kbot, rlat)
Definition: gw_drag_d.F90:912
gw_drag_d::get_uv
subroutine get_uv(uv_out, uv_in)
Definition: gw_drag_d.F90:3516
gw_drag_d::rog
real, parameter, public rog
Definition: gw_drag_d.F90:72
gw_drag_d::gw_main_d
subroutine, public gw_main_d(pcols, pver, dt, pgwv, effgworo_dev, effgwbkg_dev, pint_dev, t_dev, u_dev, u_devd, v_dev, v_devd, sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, qvt_dev, rog, mapl_vireps_, rlat_dev)
Definition: gw_drag_d.F90:83
gw_drag_d::mxasym
real, parameter, public mxasym
Definition: gw_drag_d.F90:47
gw_drag_d::get_ti
subroutine get_ti(t_out, t_in)
Definition: gw_drag_d.F90:3568
gw_drag_d::n2min
real, parameter, public n2min
Definition: gw_drag_d.F90:51
gw_drag_init::isize_e2_tau
integer, public isize_e2_tau
Definition: gw_drag_init.F90:16
max
#define max(a, b)
Definition: mosaic_util.h:33
gw_drag_d::kwvo
real, parameter, public kwvo
Definition: gw_drag_d.F90:43
gw_drag_d::fcrit2
real, parameter, public fcrit2
Definition: gw_drag_d.F90:53
gw_drag_d::orovmin
real, parameter, public orovmin
Definition: gw_drag_d.F90:57
gw_drag_d
Definition: gw_drag_d.F90:5
mapl_constantsmod::mapl_rgas
real, parameter, public mapl_rgas
Definition: MAPL_Constants.F90:45
gw_drag_init::bypass_bgnd
logical, parameter, public bypass_bgnd
Definition: gw_drag_init.F90:25
mapl_constantsmod
Definition: MAPL_Constants.F90:1
mapl_constantsmod::mapl_grav
real, parameter, public mapl_grav
Definition: MAPL_Constants.F90:18
gw_drag_d::zldvcon
real, parameter, public zldvcon
Definition: gw_drag_d.F90:71
gw_drag_d::gw_drag_prof_d
subroutine gw_drag_prof_d(i, pcols, pver, pgwv, ngwv, kbot, ktop, c, u, v, t, pi, dpm, rdpm, piln, rlat, rhoi, ni, ti, nm, dt, alpha, dback, kldv, kldvmn, ksrc, ksrcmn, rdpldv, tau, taud, ubi, ubid, ubm, xv, xvd, yv, yvd, ut, utd, vt, vtd, tt, taugwx, taugwy, fegw, fepgw, dusrc, dvsrc, dtsrc, tau0x, tau0y, effgw)
Definition: gw_drag_d.F90:1825
min
#define min(a, b)
Definition: mosaic_util.h:32
gw_drag_init::bypass_oro
logical, parameter, public bypass_oro
Definition: gw_drag_init.F90:26
gw_drag_d::get_effkwvmap_2
subroutine get_effkwvmap_2(effkwvmap_, rlat_, zfac)
Definition: gw_drag_d.F90:3536
gw_drag_d::tndmax
real, parameter, public tndmax
Definition: gw_drag_d.F90:65
gw_drag_d::orohmin
real, parameter, public orohmin
Definition: gw_drag_d.F90:55
gw_drag_d::gw_intr_d
subroutine gw_intr_d(i, pcols, pver, dt, pgwv, effgworo_dev, effgwbkg_dev, pint_dev, t_dev, u_dev, u_devd, v_dev, v_devd, sgh_dev, pref_dev, pmid_dev, pdel_dev, rpdel_dev, lnpint_dev, zm_dev, rlat_dev, dudt_gwd_dev, dudt_gwd_devd, dvdt_gwd_dev, dvdt_gwd_devd, dtdt_gwd_dev, dudt_org_dev, dvdt_org_dev, dtdt_org_dev, taugwdx_dev, taugwdy_dev, tauox_dev, tauoy_dev, feo_dev, taubkgx_dev, taubkgy_dev, taubx_dev, tauby_dev, feb_dev, fepo_dev, fepb_dev, utbsrc_dev, vtbsrc_dev, ttbsrc_dev)
Definition: gw_drag_d.F90:325
gw_drag_d::ubmc2mn
real, parameter, public ubmc2mn
Definition: gw_drag_d.F90:69
gw_drag_d::get_effkwvmap_3
subroutine get_effkwvmap_3(effkwvmap_, rlat_, zfac)
Definition: gw_drag_d.F90:3549
gw_drag_init
Definition: gw_drag_init.F90:4