Line data Source code
1 : module gravity_waves_sources
2 : use derivative_mod, only: derivative_t
3 : use dimensions_mod, only: np,nlev
4 : use edgetype_mod, only: EdgeBuffer_t
5 : use element_mod, only: element_t
6 : use hybrid_mod, only: hybrid_t
7 : use shr_kind_mod, only: r8 => shr_kind_r8
8 :
9 : implicit none
10 : private
11 : save
12 :
13 : !! gravity_waves_sources created by S Santos, 10 Aug 2011
14 : !!
15 : !! gws_src_fnct starts parallel environment and computes frontogenesis
16 : !! for use by WACCM (via dp_coupling)
17 :
18 : public :: gws_src_fnct
19 : public :: gws_src_vort
20 : public :: gws_init
21 : private :: compute_frontogenesis
22 : private :: compute_vorticity_4gw
23 :
24 : type (EdgeBuffer_t) :: edge3,edge1
25 : type (derivative_t) :: deriv
26 : real(r8) :: psurf_ref
27 :
28 : !----------------------------------------------------------------------
29 : CONTAINS
30 : !----------------------------------------------------------------------
31 :
32 26880 : subroutine gws_init(elem)
33 : use parallel_mod, only : par
34 : use edge_mod, only : initEdgeBuffer
35 : use hycoef, only : hypi
36 : use pmgrid, only : plev
37 : use thread_mod, only : horz_num_threads
38 : implicit none
39 :
40 : ! Elem will be needed for future updates to edge code
41 : type(element_t), pointer :: elem(:)
42 :
43 : ! Set up variables similar to dyn_comp and prim_driver_mod initializations
44 2304 : call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1)
45 2304 : call initEdgeBuffer(par, edge1, elem, nlev,nthreads=1)
46 :
47 2304 : psurf_ref = hypi(plev+1)
48 :
49 2304 : end subroutine gws_init
50 :
51 24576 : subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga, nphys)
52 2304 : use derivative_mod, only : derivinit
53 : use dimensions_mod, only : nelemd
54 : use dof_mod, only : UniquePoints
55 : use hybrid_mod, only : config_thread_region, get_loop_ranges
56 : use parallel_mod, only : par
57 : use ppgrid, only : pver
58 : use thread_mod, only : horz_num_threads
59 : use dimensions_mod, only : fv_nphys
60 : use cam_abortutils, only : handle_allocate_error
61 :
62 : implicit none
63 : type (element_t), intent(inout), dimension(:) :: elem
64 : integer, intent(in) :: tl, nphys, tlq
65 : real (kind=r8), intent(out) :: frontgf(nphys*nphys,pver,nelemd)
66 : real (kind=r8), intent(out) :: frontga(nphys*nphys,pver,nelemd)
67 :
68 :
69 : ! Local variables
70 : type (hybrid_t) :: hybrid
71 : integer :: nets, nete, ithr, ncols, ie, ierr
72 24576 : real(kind=r8), allocatable :: frontgf_thr(:,:,:,:)
73 24576 : real(kind=r8), allocatable :: frontga_thr(:,:,:,:)
74 :
75 :
76 : ! This does not need to be a thread private data-structure
77 24576 : call derivinit(deriv)
78 : !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,frontgf_thr,frontga_thr)
79 24576 : hybrid = config_thread_region(par,'serial')
80 24576 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
81 :
82 147456 : allocate(frontgf_thr(nphys,nphys,nlev,nets:nete), stat=ierr)
83 24576 : call handle_allocate_error(ierr, 'gws_src_fnct', 'frontgf_thr')
84 :
85 122880 : allocate(frontga_thr(nphys,nphys,nlev,nets:nete), stat=ierr)
86 24576 : call handle_allocate_error(ierr, 'gws_src_fnct', 'frontga_thr')
87 :
88 :
89 24576 : call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
90 :
91 24576 : if (fv_nphys>0) then
92 197376 : do ie=nets,nete
93 518400 : frontgf(:,:,ie) = RESHAPE(frontgf_thr(:,:,:,ie),(/nphys*nphys,nlev/))
94 542976 : frontga(:,:,ie) = RESHAPE(frontga_thr(:,:,:,ie),(/nphys*nphys,nlev/))
95 : end do
96 : else
97 0 : do ie=nets,nete
98 0 : ncols = elem(ie)%idxP%NumUniquePts
99 0 : call UniquePoints(elem(ie)%idxP, nlev, frontgf_thr(:,:,:,ie), frontgf(1:ncols,:,ie))
100 0 : call UniquePoints(elem(ie)%idxP, nlev, frontga_thr(:,:,:,ie), frontga(1:ncols,:,ie))
101 : end do
102 : end if
103 24576 : deallocate(frontga_thr)
104 24576 : deallocate(frontgf_thr)
105 :
106 : !!$OMP END PARALLEL
107 :
108 49152 : end subroutine gws_src_fnct
109 :
110 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
111 24576 : subroutine gws_src_vort(elem, tl, tlq, vort4gw, nphys)
112 24576 : use derivative_mod, only : derivinit
113 : use dimensions_mod, only : nelemd
114 : use dof_mod, only : UniquePoints
115 : use hybrid_mod, only : config_thread_region, get_loop_ranges
116 : use parallel_mod, only : par
117 : use ppgrid, only : pver
118 : use thread_mod, only : horz_num_threads
119 : use dimensions_mod, only : fv_nphys
120 : use cam_abortutils, only : handle_allocate_error
121 :
122 : implicit none
123 : type (element_t), intent(in), dimension(:) :: elem
124 : integer, intent(in) :: tl, nphys, tlq
125 :
126 : !
127 : real (kind=r8), intent(out) :: vort4gw(nphys*nphys,pver,nelemd)
128 :
129 : ! Local variables
130 : type (hybrid_t) :: hybrid
131 : integer :: nets, nete, ithr, ncols, ie, ierr
132 :
133 : !
134 24576 : real(kind=r8), allocatable :: vort4gw_thr(:,:,:,:)
135 :
136 : ! This does not need to be a thread private data-structure
137 24576 : call derivinit(deriv)
138 : !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,vort4gw_thr)
139 24576 : hybrid = config_thread_region(par,'serial')
140 24576 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
141 :
142 147456 : allocate(vort4gw_thr(nphys,nphys,nlev,nets:nete), stat=ierr)
143 24576 : call handle_allocate_error(ierr, 'gws_src_vort', 'vort4gw_thr')
144 :
145 24576 : call compute_vorticity_4gw(vort4gw_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
146 :
147 24576 : if (fv_nphys>0) then
148 197376 : do ie=nets,nete
149 542976 : vort4gw(:,:,ie) = RESHAPE(vort4gw_thr(:,:,:,ie),(/nphys*nphys,nlev/))
150 : end do
151 : else
152 0 : do ie=nets,nete
153 0 : ncols = elem(ie)%idxP%NumUniquePts
154 0 : call UniquePoints(elem(ie)%idxP, nlev, vort4gw_thr(:,:,:,ie), vort4gw(1:ncols,:,ie))
155 : end do
156 : end if
157 24576 : deallocate(vort4gw_thr)
158 :
159 : !!$OMP END PARALLEL
160 :
161 49152 : end subroutine gws_src_vort
162 :
163 24576 : subroutine compute_vorticity_4gw(vort4gw,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys)
164 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165 : ! compute vorticity for use in gw params
166 : ! F = ( curl ) [U,V]
167 : !
168 : ! Original by Peter Lauritzen, Julio Bacmeister*, Dec 2024
169 : ! Patterned on 'compute_frontogenesis'
170 : !
171 : ! * corresponding/blame-able
172 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173 24576 : use derivative_mod, only: vorticity_sphere
174 : use edge_mod, only: edgevpack, edgevunpack
175 : use bndry_mod, only: bndry_exchange
176 : use dimensions_mod, only: fv_nphys
177 : use fvm_mapping, only: dyn2phys
178 :
179 : type(hybrid_t), intent(in) :: hybrid
180 : type(element_t), intent(in) :: elem(:)
181 : type(derivative_t), intent(in) :: ederiv
182 : integer, intent(in) :: nets,nete,nphys
183 : integer, intent(in) :: tl,tlq
184 : real(r8), intent(out) :: vort4gw(nphys,nphys,nlev,nets:nete)
185 :
186 : ! local
187 49152 : real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
188 49152 : real(r8) :: vort_gll(np,np,nlev,nets:nete)
189 : integer :: k,kptr,i,j,ie,component,h,nq,m_cnst,n0
190 :
191 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
192 : ! First calculate vorticity on GLL grid
193 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
194 : ! set timelevel=1 fro velocities
195 24576 : n0=tl
196 197376 : do ie=nets,nete
197 16243200 : do k=1,nlev
198 16243200 : call vorticity_sphere(elem(ie)%state%v(:,:,:,k,n0),ederiv,elem(ie),vort_gll(:,:,k,ie))
199 : end do
200 16243200 : do k=1,nlev
201 337651200 : vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
202 : end do
203 : ! pack
204 197376 : call edgeVpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie)
205 : enddo
206 24576 : call bndry_exchange(hybrid,edge1,location='compute_vorticity_4gw')
207 197376 : do ie=nets,nete
208 172800 : call edgeVunpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie)
209 : ! apply inverse mass matrix,
210 16243200 : do k=1,nlev
211 337651200 : vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:)
212 : end do
213 :
214 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
215 : ! Now regrid from GLL to PhysGrid if necessary
216 : ! otherwise just return vorticity on GLL grid
217 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
218 197376 : if (fv_nphys>0) then
219 3628800 : tmp = 1.0_r8
220 172800 : area_inv = dyn2phys(tmp,elem(ie)%metdet)
221 2246400 : area_inv = 1.0_r8/area_inv
222 16243200 : do k=1,nlev
223 16243200 : vort4gw(:,:,k,ie) = dyn2phys( vort_gll(:,:,k,ie) , elem(ie)%metdet , area_inv )
224 : end do
225 : else
226 0 : do k=1,nlev
227 0 : vort4gw(:,:,k,ie) = vort_gll(:,:,k,ie)
228 : end do
229 : end if
230 : enddo
231 :
232 :
233 24576 : end subroutine compute_vorticity_4gw
234 :
235 :
236 24576 : subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys)
237 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
238 : ! compute frontogenesis function F
239 : ! F = -gradth dot C
240 : ! with:
241 : ! theta = potential temperature
242 : ! gradth = grad(theta)
243 : ! C = ( gradth dot grad ) U
244 : !
245 : ! Original by Mark Taylor, July 2011
246 : ! Change by Santos, 10 Aug 2011:
247 : ! Integrated into gravity_waves_sources module, several arguments made global
248 : ! to prevent repeated allocation/initialization
249 : !
250 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
251 24576 : use physconst, only: cappa
252 : use air_composition,only: dry_air_species_num, thermodynamic_active_species_num
253 : use air_composition,only: thermodynamic_active_species_idx_dycore
254 : use derivative_mod, only: gradient_sphere, ugradv_sphere
255 : use edge_mod, only: edgevpack, edgevunpack
256 : use bndry_mod, only: bndry_exchange
257 : use dyn_grid, only: hvcoord
258 : use dimensions_mod, only: fv_nphys,ntrac
259 : use fvm_mapping, only: dyn2phys_vector,dyn2phys
260 :
261 : type(hybrid_t), intent(in) :: hybrid
262 : type(element_t), intent(inout), target :: elem(:)
263 : type(derivative_t), intent(in) :: ederiv
264 : integer, intent(in) :: nets,nete,nphys
265 : integer, intent(in) :: tl,tlq
266 : real(r8), intent(out) :: frontgf(nphys,nphys,nlev,nets:nete)
267 : real(r8), intent(out) :: frontga(nphys,nphys,nlev,nets:nete)
268 :
269 : ! local
270 49152 : real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
271 49152 : real(r8) :: uv_tmp(fv_nphys*fv_nphys,2,nlev)
272 49152 : real(r8) :: frontgf_gll(np,np,nlev,nets:nete)
273 : real(r8) :: frontga_gll(np,np,nlev,nets:nete)
274 : integer :: k,kptr,i,j,ie,component,h,nq,m_cnst
275 49152 : real(r8) :: gradth(np,np,2,nlev,nets:nete) ! grad(theta)
276 : real(r8) :: p(np,np) ! pressure at mid points
277 : real(r8) :: pint(np,np) ! pressure at interface points
278 : real(r8) :: theta(np,np) ! potential temperature at mid points
279 : real(r8) :: C(np,np,2), sum_water(np,np)
280 :
281 197376 : do ie=nets,nete
282 : ! pressure at model top
283 3628800 : pint(:,:) = hvcoord%hyai(1)*hvcoord%ps0
284 16243200 : do k=1,nlev
285 : ! moist pressure at mid points
286 337478400 : sum_water(:,:) = 1.0_r8
287 112492800 : do nq=dry_air_species_num+1,thermodynamic_active_species_num
288 96422400 : m_cnst = thermodynamic_active_species_idx_dycore(nq)
289 : !
290 : ! make sure Q is updated
291 : !
292 2040940800 : sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl)
293 : end do
294 337478400 : p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
295 : ! moist pressure at interface for next iteration
296 337478400 : pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl)
297 : !
298 337478400 : theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa
299 : ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv)
300 16070400 : call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
301 : ! compute C = (grad(theta) dot grad ) u
302 16070400 : C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie))
303 : ! gradth dot C
304 337478400 : frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) + C(:,:,2)*gradth(:,:,2,k,ie) )
305 : ! apply mass matrix
306 337478400 : gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
307 337478400 : gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
308 337651200 : frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
309 : enddo
310 : ! pack
311 172800 : call edgeVpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
312 197376 : call edgeVpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
313 : enddo
314 24576 : call bndry_exchange(hybrid,edge3,location='compute_frontogenesis')
315 197376 : do ie=nets,nete
316 172800 : call edgeVunpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
317 172800 : call edgeVunpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
318 : ! apply inverse mass matrix,
319 16243200 : do k=1,nlev
320 337478400 : gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
321 337478400 : gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
322 337651200 : frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:)
323 : end do
324 197376 : if (fv_nphys>0) then
325 172800 : uv_tmp(:,:,:) = dyn2phys_vector(gradth(:,:,:,:,ie),elem(ie))
326 16243200 : do k=1,nlev
327 16070400 : h=0
328 64454400 : do j=1,fv_nphys
329 208915200 : do i=1,fv_nphys
330 144633600 : h=h+1
331 192844800 : frontga(i,j,k,ie) = atan2 ( uv_tmp(h,2,k) , uv_tmp(h,1,k) + 1.e-10_r8 )
332 : end do
333 : end do
334 : end do
335 : !
336 : ! compute inverse physgrid area for mapping of scaler
337 : !
338 3628800 : tmp = 1.0_r8
339 172800 : area_inv = dyn2phys(tmp,elem(ie)%metdet)
340 2246400 : area_inv = 1.0_r8/area_inv
341 16243200 : do k=1,nlev
342 16243200 : frontgf(:,:,k,ie) = dyn2phys(frontgf_gll(:,:,k,ie),elem(ie)%metdet,area_inv)
343 : end do
344 : else
345 0 : do k=1,nlev
346 0 : frontgf(:,:,k,ie)=frontgf_gll(:,:,k,ie)
347 : ! Frontogenesis angle
348 0 : frontga(:,:,k,ie) = atan2 ( gradth(:,:,2,k,ie) , gradth(:,:,1,k,ie) + 1.e-10_r8 )
349 : end do
350 : end if
351 : enddo
352 24576 : end subroutine compute_frontogenesis
353 :
354 :
355 : end module gravity_waves_sources
|
