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 0 : 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 0 : call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1)
45 0 : call initEdgeBuffer(par, edge1, elem, nlev,nthreads=1)
46 :
47 0 : psurf_ref = hypi(plev+1)
48 :
49 0 : end subroutine gws_init
50 :
51 0 : subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga, nphys)
52 0 : 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 0 : real(kind=r8), allocatable :: frontgf_thr(:,:,:,:)
73 0 : real(kind=r8), allocatable :: frontga_thr(:,:,:,:)
74 :
75 :
76 : ! This does not need to be a thread private data-structure
77 0 : call derivinit(deriv)
78 : !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,frontgf_thr,frontga_thr)
79 0 : hybrid = config_thread_region(par,'serial')
80 0 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
81 :
82 0 : allocate(frontgf_thr(nphys,nphys,nlev,nets:nete), stat=ierr)
83 0 : call handle_allocate_error(ierr, 'gws_src_fnct', 'frontgf_thr')
84 :
85 0 : allocate(frontga_thr(nphys,nphys,nlev,nets:nete), stat=ierr)
86 0 : call handle_allocate_error(ierr, 'gws_src_fnct', 'frontga_thr')
87 :
88 :
89 0 : call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
90 :
91 0 : if (fv_nphys>0) then
92 0 : do ie=nets,nete
93 0 : frontgf(:,:,ie) = RESHAPE(frontgf_thr(:,:,:,ie),(/nphys*nphys,nlev/))
94 0 : 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 0 : deallocate(frontga_thr)
104 0 : deallocate(frontgf_thr)
105 :
106 : !!$OMP END PARALLEL
107 :
108 0 : end subroutine gws_src_fnct
109 :
110 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
111 0 : subroutine gws_src_vort(elem, tl, tlq, vort4gw, nphys)
112 0 : 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 0 : real(kind=r8), allocatable :: vort4gw_thr(:,:,:,:)
135 :
136 : ! This does not need to be a thread private data-structure
137 0 : call derivinit(deriv)
138 : !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,vort4gw_thr)
139 0 : hybrid = config_thread_region(par,'serial')
140 0 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
141 :
142 0 : allocate(vort4gw_thr(nphys,nphys,nlev,nets:nete), stat=ierr)
143 0 : call handle_allocate_error(ierr, 'gws_src_vort', 'vort4gw_thr')
144 :
145 0 : call compute_vorticity_4gw(vort4gw_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
146 :
147 0 : if (fv_nphys>0) then
148 0 : do ie=nets,nete
149 0 : 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 0 : deallocate(vort4gw_thr)
158 :
159 : !!$OMP END PARALLEL
160 :
161 0 : end subroutine gws_src_vort
162 :
163 0 : 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 0 : 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 0 : real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
188 0 : 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 0 : n0=tl
196 0 : do ie=nets,nete
197 0 : do k=1,nlev
198 0 : call vorticity_sphere(elem(ie)%state%v(:,:,:,k,n0),ederiv,elem(ie),vort_gll(:,:,k,ie))
199 : end do
200 0 : do k=1,nlev
201 0 : vort_gll(:,:,k,ie) = vort_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
202 : end do
203 : ! pack
204 0 : call edgeVpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie)
205 : enddo
206 0 : call bndry_exchange(hybrid,edge1,location='compute_vorticity_4gw')
207 0 : do ie=nets,nete
208 0 : call edgeVunpack(edge1, vort_gll(:,:,:,ie),nlev,0,ie)
209 : ! apply inverse mass matrix,
210 0 : do k=1,nlev
211 0 : 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 0 : if (fv_nphys>0) then
219 0 : tmp = 1.0_r8
220 0 : area_inv = dyn2phys(tmp,elem(ie)%metdet)
221 0 : area_inv = 1.0_r8/area_inv
222 0 : do k=1,nlev
223 0 : 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 0 : end subroutine compute_vorticity_4gw
234 :
235 :
236 0 : 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 0 : 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 0 : real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
271 0 : real(r8) :: uv_tmp(fv_nphys*fv_nphys,2,nlev)
272 0 : 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 0 : 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 0 : do ie=nets,nete
282 : ! pressure at model top
283 0 : pint(:,:) = hvcoord%hyai(1)*hvcoord%ps0
284 0 : do k=1,nlev
285 : ! moist pressure at mid points
286 0 : sum_water(:,:) = 1.0_r8
287 0 : do nq=dry_air_species_num+1,thermodynamic_active_species_num
288 0 : m_cnst = thermodynamic_active_species_idx_dycore(nq)
289 : !
290 : ! make sure Q is updated
291 : !
292 0 : sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl)
293 : end do
294 0 : p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
295 : ! moist pressure at interface for next iteration
296 0 : pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl)
297 : !
298 0 : theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa
299 : ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv)
300 0 : call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
301 : ! compute C = (grad(theta) dot grad ) u
302 0 : C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie))
303 : ! gradth dot C
304 0 : frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) + C(:,:,2)*gradth(:,:,2,k,ie) )
305 : ! apply mass matrix
306 0 : gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
307 0 : gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
308 0 : frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
309 : enddo
310 : ! pack
311 0 : call edgeVpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
312 0 : call edgeVpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
313 : enddo
314 0 : call bndry_exchange(hybrid,edge3,location='compute_frontogenesis')
315 0 : do ie=nets,nete
316 0 : call edgeVunpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
317 0 : call edgeVunpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
318 : ! apply inverse mass matrix,
319 0 : do k=1,nlev
320 0 : gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
321 0 : gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
322 0 : frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:)
323 : end do
324 0 : if (fv_nphys>0) then
325 0 : uv_tmp(:,:,:) = dyn2phys_vector(gradth(:,:,:,:,ie),elem(ie))
326 0 : do k=1,nlev
327 0 : h=0
328 0 : do j=1,fv_nphys
329 0 : do i=1,fv_nphys
330 0 : h=h+1
331 0 : 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 0 : tmp = 1.0_r8
339 0 : area_inv = dyn2phys(tmp,elem(ie)%metdet)
340 0 : area_inv = 1.0_r8/area_inv
341 0 : do k=1,nlev
342 0 : 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 0 : end subroutine compute_frontogenesis
353 :
354 :
355 : end module gravity_waves_sources
|
