Line data Source code
1 : module gravity_waves_sources
2 :
3 : use shr_kind_mod, only: r8 => shr_kind_r8
4 : use pmgrid, only: plev, plevp
5 : use hycoef, only : hypi
6 :
7 : implicit none
8 : save
9 : private
10 :
11 : public :: gws_src_fnct
12 :
13 : ! ghosting added by Francis Vitt -- 7 July 2008
14 : !
15 : ! moved from waccm to fv, changed source of psurf_ref
16 : ! -- S Santos -- 10 Aug 2011
17 :
18 : contains
19 :
20 :
21 : !=================================================================
22 :
23 8064 : subroutine gws_src_fnct(grid, u3,v3,pt, q3, pe, frontgf, frontga)
24 :
25 : use dynamics_vars, only: t_fvdycore_grid
26 : use physconst, only: zvir, cappa, aearth => rearth
27 :
28 : implicit none
29 :
30 : ! Input/Output arguments
31 : type (t_fvdycore_grid), intent(in) :: grid ! grid for XY decomp
32 : real(r8), intent(in) :: u3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! zonal velocity
33 : real(r8), intent(in) :: v3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! meridional velocity
34 : real(r8), intent(in) :: pt(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,plev) ! virtual temperature
35 : real(r8), intent(in) :: q3(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,plev) ! water constituent
36 : real(r8), intent(in) :: pe(grid%ifirstxy:grid%ilastxy,plevp,grid%jfirstxy:grid%jlastxy) ! interface pressure
37 :
38 : real(r8), intent(out) :: frontgf(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! Frontogenesis function
39 : real(r8), intent(out) :: frontga(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! Frontogenesis angle
40 :
41 : ! Locals
42 : real(r8) :: psurf_ref ! surface reference pressure
43 :
44 16128 : real(r8) :: ptemp(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! temperature
45 16128 : real(r8) :: pm(grid%ifirstxy:grid%ilastxy ,plev,grid%jfirstxy:grid%jlastxy) ! mid-point pressure
46 : real(r8) :: pexf ! Exner function
47 :
48 16128 : real(r8) :: pty(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! temperature meridional gradient
49 16128 : real(r8) :: ptx(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! temperature zonal gradient
50 :
51 16128 : real(r8) :: uy(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! U-wind meridional gradient
52 16128 : real(r8) :: ux(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! U-wind zonal gradient
53 :
54 16128 : real(r8) :: vy(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! V-wind meridional gradient
55 16128 : real(r8) :: vx(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! V-wind zonal gradient
56 :
57 16128 : real(r8) :: ptg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! temperature ghosted
58 16128 : real(r8) :: ug(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! U-wind ghosted
59 16128 : real(r8) :: vg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! V-wind ghosted
60 :
61 : real(r8) :: tglat ! tangent-latitude
62 : integer :: i,j,k
63 : integer :: im, ip
64 : integer :: beglatxy, endlatxy, beglonxy, endlonxy
65 : !-----------------------------------------------------------------------------------------
66 :
67 8064 : beglatxy = grid%jfirstxy
68 8064 : endlatxy = grid%jlastxy
69 8064 : beglonxy = grid%ifirstxy
70 8064 : endlonxy = grid%ilastxy
71 :
72 40916736 : pty(:,:,:) = 0._r8
73 40916736 : ptx(:,:,:) = 0._r8
74 40916736 : uy(:,:,:) = 0._r8
75 40916736 : ux(:,:,:) = 0._r8
76 40916736 : vy(:,:,:) = 0._r8
77 40916736 : vx(:,:,:) = 0._r8
78 40916736 : frontgf(:,:,:) = 0._r8
79 40916736 : frontga(:,:,:) = 0._r8
80 :
81 8064 : psurf_ref = hypi(plev+1)
82 :
83 : !$omp parallel do private (i,j,k,pexf)
84 32256 : do j = beglatxy, endlatxy
85 3177216 : do k = 1, plev
86 40908672 : do i = beglonxy, endlonxy
87 :
88 : ! Calculate pressure and Exner function
89 37739520 : pm(i,k,j) = 0.5_r8 * ( pe(i,k,j) + pe(i,k+1,j) )
90 37739520 : pexf = (psurf_ref / pm(i,k,j))**cappa
91 :
92 : ! Convert virtual temperature to temperature and calculate potential temperature
93 40884480 : ptemp(i,k,j) = pt(i,j,k) / (1._r8 + zvir*q3(i,j,k)) * pexf
94 :
95 : end do
96 : end do
97 : end do
98 :
99 8064 : call ghost_array(grid, ptemp, ptg)
100 8064 : call ghost_array(grid, u3, ug)
101 8064 : call ghost_array(grid, v3, vg)
102 :
103 : !$omp parallel do private (i,j,k)
104 1056384 : do k=1, plev
105 13636224 : do i=beglonxy, endlonxy
106 51367680 : do j=beglatxy, endlatxy
107 :
108 : ! Pot. Temperature
109 37739520 : pty(i,k,j) = ( ptg(i,k,j+1) - ptg(i,k,j-1) ) / (2._r8 * grid%dp)
110 37739520 : pty(i,k,j) = pty(i,k,j) / aearth
111 :
112 : ! U-wind
113 37739520 : uy(i,k,j) = ( ug(i,k,j+1) - ug(i,k,j-1) ) / (2._r8 * grid%dp)
114 37739520 : uy(i,k,j) = uy(i,k,j) / aearth
115 :
116 : ! V-wind
117 37739520 : vy(i,k,j) = ( vg(i,k,j+1) - vg(i,k,j-1) ) / (2._r8 * grid%dp)
118 50319360 : vy(i,k,j) = vy(i,k,j) / aearth
119 :
120 : end do
121 : end do
122 : end do
123 :
124 : !++rrg use 1.e-3 floor on cosine terms in the denominator of frontgf
125 :
126 : !$omp parallel do private (i,j,k,im,ip)
127 1056384 : do k=1, plev
128 4201344 : do j=beglatxy, endlatxy
129 41932800 : do i=beglonxy, endlonxy
130 :
131 37739520 : im = i-1
132 37739520 : ip = i+1
133 :
134 : ! Pot. Temperature
135 37739520 : ptx(i,k,j) = ( ptg(ip,k,j) - ptg(im,k,j) ) / (2._r8 * grid%dl)
136 37739520 : ptx(i,k,j) = ptx(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))
137 :
138 : ! U-wind
139 37739520 : ux(i,k,j) = ( ug(ip,k,j) - ug(im,k,j) ) / (2._r8 *grid%dl)
140 37739520 : ux(i,k,j) = ux(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))
141 :
142 : ! V-wind
143 37739520 : vx(i,k,j) = ( vg(ip,k,j) - vg(im,k,j) ) / (2._r8 *grid%dl)
144 40884480 : vx(i,k,j) = vx(i,k,j) / (aearth * (grid%cosp(j)+1.e-3_r8))
145 :
146 : end do
147 : end do
148 : end do
149 :
150 : !$omp parallel do private (i,j,k, tglat)
151 32256 : do j=beglatxy, endlatxy
152 :
153 24192 : tglat = grid%sinp(j) / (grid%cosp(j)+1.e-3_r8)
154 :
155 3177216 : do k=1, plev
156 40908672 : do i=beglonxy, endlonxy
157 :
158 113218560 : frontgf(i,k,j) = &
159 : - ptx(i,k,j)**2._r8 * (ux(i,k,j) - v3(i,k,j) * tglat / aearth) &
160 : - pty(i,k,j)**2._r8 * vy(i,k,j) &
161 116363520 : - ptx(i,k,j) * pty(i,k,j) * ( vx(i,k,j) + uy(i,k,j) + u3(i,k,j) * tglat / aearth )
162 :
163 : end do
164 : end do
165 :
166 : end do
167 :
168 : !--rrg use 1.e-3 floor on cosine terms in the denominator of frontgf
169 :
170 : !$omp parallel do private (i,j,k)
171 32256 : do j=beglatxy, endlatxy
172 3177216 : do k=1, plev
173 40908672 : do i=beglonxy, endlonxy
174 40884480 : frontga(i,k,j) = atan2 ( pty(i,k,j) , ptx(i,k,j) + 1.e-10_r8 )
175 : end do
176 : end do
177 : end do
178 :
179 8064 : return
180 :
181 : end subroutine gws_src_fnct
182 :
183 24192 : subroutine ghost_array(grid, x, xg)
184 :
185 : ! subroutine added by Francis Vitt -- 7 July 2008
186 :
187 : #if defined( SPMD )
188 : use mod_comm, only: mp_send3d, mp_recv3d
189 : #endif
190 : use dynamics_vars, only: T_FVDYCORE_GRID
191 :
192 : implicit none
193 :
194 : ! Input/Output arguments
195 : type (T_FVDYCORE_GRID), intent(in) :: grid ! grid for XY decomp
196 : real(r8), intent(in) :: x(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! zonal velocity
197 : real(r8), intent(out) :: xg(grid%ifirstxy-1:grid%ilastxy+1,plev,grid%jfirstxy-1:grid%jlastxy+1) ! zonal velocity
198 :
199 : ! local variables
200 48384 : real(r8) :: north(grid%ifirstxy:grid%ilastxy,plev)
201 48384 : real(r8) :: south(grid%ifirstxy:grid%ilastxy,plev)
202 48384 : real(r8) :: east(plev,grid%jfirstxy:grid%jlastxy)
203 48384 : real(r8) :: west(plev,grid%jfirstxy:grid%jlastxy)
204 : integer :: im, jm, km, ifirstxy, ilastxy, jfirstxy, jlastxy, iam, myidxy_y, nprxy_x
205 : integer :: itot, dest, src, j, k
206 :
207 24192 : im = grid%im
208 24192 : jm = grid%jm
209 24192 : km = grid%km
210 :
211 24192 : ifirstxy = grid%ifirstxy
212 24192 : ilastxy = grid%ilastxy
213 24192 : jfirstxy = grid%jfirstxy
214 24192 : jlastxy = grid%jlastxy
215 :
216 24192 : iam = grid%iam
217 24192 : myidxy_y = grid%myidxy_y
218 24192 : nprxy_x = grid%nprxy_x
219 24192 : itot = ilastxy-ifirstxy+1
220 :
221 122774400 : xg(ifirstxy:ilastxy,:,jfirstxy:jlastxy) = x(ifirstxy:ilastxy,:,jfirstxy:jlastxy)
222 :
223 : #if defined( SPMD )
224 :
225 : ! north
226 : call mp_send3d( grid%commxy, iam-nprxy_x, iam+nprxy_x, im, km, jm, &
227 : ifirstxy, ilastxy, 1, km, jfirstxy, jlastxy, &
228 24192 : ifirstxy, ilastxy, 1, km, jfirstxy, jfirstxy, x )
229 : call mp_recv3d( grid%commxy, iam+nprxy_x, im, jm, km, &
230 : ifirstxy, ilastxy, 1, km, jlastxy+1, jlastxy+1, &
231 24192 : ifirstxy, ilastxy, 1, km, jlastxy+1, jlastxy+1, north )
232 :
233 : ! south
234 : call mp_send3d( grid%commxy, iam+nprxy_x, iam-nprxy_x, im, km, jm, &
235 : ifirstxy, ilastxy, 1, km, jfirstxy, jlastxy, &
236 24192 : ifirstxy, ilastxy, 1, km, jlastxy, jlastxy, x )
237 : call mp_recv3d( grid%commxy, iam-nprxy_x, im, jm, km, &
238 : ifirstxy, ilastxy, 1, km, jfirstxy-1, jfirstxy-1, &
239 24192 : ifirstxy, ilastxy, 1, km, jfirstxy-1, jfirstxy-1, south )
240 :
241 : #endif
242 :
243 24192 : if (itot .ne. im) then
244 : #if defined( SPMD )
245 :
246 : ! east
247 :
248 24192 : dest = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
249 24192 : src = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
250 : call mp_send3d( grid%commxy, dest, src, im, km, jm, &
251 : ifirstxy, ilastxy, 1, km, jfirstxy, jlastxy, &
252 24192 : ifirstxy, ifirstxy, 1, km, jfirstxy, jlastxy, x )
253 : call mp_recv3d( grid%commxy, src, im, km, jm, &
254 : ilastxy+1, ilastxy+1, 1, km, jfirstxy, jlastxy, &
255 24192 : ilastxy+1, ilastxy+1, 1, km, jfirstxy, jlastxy, east )
256 :
257 : ! west
258 :
259 24192 : dest = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
260 24192 : src = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
261 : call mp_send3d( grid%commxy, dest, src, im, km, jm, &
262 : ifirstxy, ilastxy, 1, km, jfirstxy, jlastxy, &
263 24192 : ilastxy, ilastxy, 1, km, jfirstxy, jlastxy, x )
264 : call mp_recv3d( grid%commxy, src, im, km, jm, &
265 : ifirstxy-1, ifirstxy-1, 1, km, jfirstxy, jlastxy, &
266 24192 : ifirstxy-1, ifirstxy-1, 1, km, jfirstxy, jlastxy, west )
267 : #endif
268 :
269 : else
270 : !$omp parallel do private(j, k)
271 0 : do k = 1,km
272 0 : do j=jfirstxy,jlastxy
273 0 : east(k,j) = x(1, k,j)
274 0 : west(k,j) = x(im,k,j)
275 : enddo
276 : enddo
277 : endif
278 :
279 24192 : if ( jfirstxy == 1 ) then
280 1278396 : xg(ifirstxy:ilastxy,:,jfirstxy-1) = xg(ifirstxy:ilastxy,:,jfirstxy)
281 : else
282 39630276 : xg(ifirstxy:ilastxy,:,jfirstxy-1) = south
283 : endif
284 :
285 24192 : if ( jlastxy == jm ) then
286 1278396 : xg(ifirstxy:ilastxy,:,jlastxy+1) = xg(ifirstxy:ilastxy,:,jlastxy)
287 : else
288 39630276 : xg(ifirstxy:ilastxy,:,jlastxy+1) = north
289 : endif
290 :
291 9531648 : xg(ifirstxy-1,:,jfirstxy:jlastxy) = west
292 9531648 : xg( ilastxy+1,:,jfirstxy:jlastxy) = east
293 :
294 24192 : end subroutine ghost_array
295 : !=================================================================
296 :
297 : end module gravity_waves_sources
|
