Line data Source code
1 : module gw_diffusion
2 :
3 : !
4 : ! This module contains code computing the effective diffusion of
5 : ! constituents and dry static energy due to gravity wave breaking.
6 : !
7 :
8 : use gw_utils, only: r8
9 : use linear_1d_operators, only: TriDiagDecomp
10 :
11 : implicit none
12 : private
13 : save
14 :
15 : public :: gw_ediff
16 : public :: gw_diff_tend
17 :
18 : contains
19 :
20 : !==========================================================================
21 :
22 1489176 : subroutine gw_ediff(ncol, pver, ngwv, kbot, ktop, tend_level, &
23 2978352 : gwut, ubm, nm, rho, dt, prndl, gravit, p, c, vramp, &
24 0 : egwdffi, decomp, ro_adjust)
25 : !
26 : ! Calculate effective diffusivity associated with GW forcing.
27 : !
28 : ! Author: F. Sassi, Jan 31, 2001
29 : !
30 : use gw_utils, only: midpoint_interp
31 : use coords_1d, only: Coords1D
32 : use vdiff_lu_solver, only: fin_vol_lu_decomp
33 :
34 : !-------------------------------Input Arguments----------------------------
35 :
36 : ! Column, level, and gravity wave spectrum dimensions.
37 : integer, intent(in) :: ncol, pver, ngwv
38 : ! Bottom and top levels to operate on.
39 : integer, intent(in) :: kbot, ktop
40 : ! Per-column bottom index where tendencies are applied.
41 : integer, intent(in) :: tend_level(ncol)
42 : ! GW zonal wind tendencies at midpoint.
43 : real(r8), intent(in) :: gwut(ncol,pver,-ngwv:ngwv)
44 : ! Projection of wind at midpoints.
45 : real(r8), intent(in) :: ubm(ncol,pver)
46 : ! Brunt-Vaisalla frequency.
47 : real(r8), intent(in) :: nm(ncol,pver)
48 :
49 : ! Density at interfaces.
50 : real(r8), intent(in) :: rho(ncol,pver+1)
51 : ! Time step.
52 : real(r8), intent(in) :: dt
53 : ! Inverse Prandtl number.
54 : real(r8), intent(in) :: prndl
55 : ! Acceleration due to gravity.
56 : real(r8), intent(in) :: gravit
57 : ! Pressure coordinates.
58 : type(Coords1D), intent(in) :: p
59 : ! Wave phase speeds for each column.
60 : real(r8), intent(in) :: c(ncol,-ngwv:ngwv)
61 :
62 : ! Coefficient to ramp down diffusion coeff.
63 : real(r8), pointer, intent(in) :: vramp(:)
64 :
65 : ! Adjustment parameter for IGWs.
66 : real(r8), intent(in), optional :: &
67 : ro_adjust(ncol,-ngwv:ngwv,pver+1)
68 :
69 : !-----------------------------Output Arguments-----------------------------
70 : ! Effective gw diffusivity at interfaces.
71 : real(r8), intent(out) :: egwdffi(ncol,pver+1)
72 : ! LU decomposition.
73 : type(TriDiagDecomp), intent(out) :: decomp
74 :
75 : !-----------------------------Local Workspace------------------------------
76 :
77 : ! Effective gw diffusivity at midpoints.
78 2978352 : real(r8) :: egwdffm(ncol,pver)
79 : ! Temporary used to hold gw_diffusivity for one level and wavenumber.
80 2978352 : real(r8) :: egwdff_lev(ncol)
81 : ! (dp/dz)^2 == (gravit*rho)^2
82 1489176 : real(r8) :: dpidz_sq(ncol,pver+1)
83 : ! Level and wave indices.
84 : integer :: k, l
85 :
86 : ! Density scale height.
87 : real(r8), parameter :: dscale=7000._r8
88 :
89 : !--------------------------------------------------------------------------
90 :
91 672865128 : egwdffi = 0._r8
92 647999352 : egwdffm = 0._r8
93 :
94 : ! Calculate effective diffusivity at midpoints.
95 2978352 : do l = -ngwv, ngwv
96 41696928 : do k = ktop, kbot
97 :
98 : egwdff_lev = &
99 646510176 : prndl * 0.5_r8 * gwut(:,k,l) * (c(:,l)-ubm(:,k)) / nm(:,k)**2
100 :
101 : ! IGWs need ro_adjust factor.
102 38718576 : if (present(ro_adjust)) then
103 0 : egwdff_lev = egwdff_lev * ro_adjust(:,l,k)**2
104 : end if
105 :
106 647999352 : egwdffm(:,k) = egwdffm(:,k) + egwdff_lev
107 :
108 : end do
109 : end do
110 :
111 1489176 : if (associated(vramp)) then
112 0 : do k = ktop,kbot
113 0 : egwdffm(:,k) = egwdffm(:,k) * vramp(k)
114 : end do
115 : endif
116 :
117 : ! Interpolate effective diffusivity to interfaces.
118 : ! Assume zero at top and bottom interfaces.
119 1489176 : egwdffi(:,ktop+1:kbot) = midpoint_interp(egwdffm(:,ktop:kbot))
120 :
121 : ! Do not calculate diffusivities below level where tendencies are
122 : ! actually allowed.
123 38718576 : do k = ktop+1, kbot
124 623133576 : where (k > tend_level) egwdffi(:,k) = 0.0_r8
125 : enddo
126 :
127 : ! Calculate (dp/dz)^2.
128 672865128 : dpidz_sq = rho*gravit
129 672865128 : dpidz_sq = dpidz_sq*dpidz_sq
130 :
131 : ! Decompose the diffusion matrix.
132 : decomp = fin_vol_lu_decomp(dt, p%section([1,ncol],[ktop,kbot]), &
133 678821832 : coef_q_diff=egwdffi(:,ktop:kbot+1)*dpidz_sq(:,ktop:kbot+1))
134 :
135 1489176 : end subroutine gw_ediff
136 :
137 : !==========================================================================
138 :
139 5956704 : subroutine gw_diff_tend(ncol, pver, kbot, ktop, q, dt, decomp, dq)
140 :
141 : !
142 : ! Calculates tendencies from effective diffusion due to gravity wave
143 : ! breaking.
144 : !
145 : ! Method:
146 : ! A constituent flux on interfaces is given by:
147 : !
148 : ! rho * (w'q') = rho * Deff qz
149 : !
150 : ! where (all evaluated on interfaces):
151 : !
152 : ! rho = density
153 : ! qz = constituent vertical gradient
154 : ! Deff = effective diffusivity
155 : !
156 : ! An effective diffusivity is calculated by adding up the diffusivities
157 : ! from all waves (see gw_ediff). The tendency is calculated by invoking LU
158 : ! decomposition and solving as for a regular diffusion equation.
159 : !
160 : ! Author: Sassi - Jan 2001
161 : !--------------------------------------------------------------------------
162 :
163 : !---------------------------Input Arguments--------------------------------
164 :
165 : ! Column and level dimensions.
166 : integer, intent(in) :: ncol, pver
167 : ! Bottom and top levels to operate on.
168 : integer, intent(in) :: kbot, ktop
169 :
170 : ! Constituent to diffuse.
171 : real(r8), intent(in) :: q(ncol,pver)
172 : ! Time step.
173 : real(r8), intent(in) :: dt
174 :
175 : ! LU decomposition.
176 : type(TriDiagDecomp), intent(in) :: decomp
177 :
178 : !--------------------------Output Arguments--------------------------------
179 :
180 : ! Constituent tendencies.
181 : real(r8), intent(out) :: dq(ncol,pver)
182 :
183 : !--------------------------Local Workspace---------------------------------
184 :
185 : ! Temporary storage for constituent.
186 5956704 : real(r8) :: qnew(ncol,pver)
187 :
188 : !--------------------------------------------------------------------------
189 :
190 2591997408 : dq = 0.0_r8
191 2591997408 : qnew = q
192 :
193 5956704 : call decomp%left_div(qnew(:,ktop:kbot))
194 :
195 : ! Evaluate tendency to be reported back.
196 2591997408 : dq = (qnew-q) / dt
197 :
198 5956704 : end subroutine gw_diff_tend
199 :
200 : end module gw_diffusion
|
