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_init
20 : private :: compute_frontogenesis
21 :
22 : type (EdgeBuffer_t) :: edge3
23 : type (derivative_t) :: deriv
24 : real(r8) :: psurf_ref
25 :
26 : !----------------------------------------------------------------------
27 : CONTAINS
28 : !----------------------------------------------------------------------
29 :
30 372480 : subroutine gws_init(elem)
31 : use parallel_mod, only : par
32 : use edge_mod, only : initEdgeBuffer
33 : use hycoef, only : hypi
34 : use pmgrid, only : plev
35 : use thread_mod, only : horz_num_threads
36 : implicit none
37 :
38 : ! Elem will be needed for future updates to edge code
39 : type(element_t), pointer :: elem(:)
40 :
41 : ! Set up variables similar to dyn_comp and prim_driver_mod initializations
42 1536 : call initEdgeBuffer(par, edge3, elem, 3*nlev,nthreads=1)
43 :
44 1536 : psurf_ref = hypi(plev+1)
45 :
46 1536 : end subroutine gws_init
47 :
48 370944 : subroutine gws_src_fnct(elem, tl, tlq, frontgf, frontga,nphys)
49 1536 : use derivative_mod, only : derivinit
50 : use dimensions_mod, only : npsq, nelemd
51 : use dof_mod, only : UniquePoints
52 : use hybrid_mod, only : config_thread_region, get_loop_ranges
53 : use parallel_mod, only : par
54 : use ppgrid, only : pver
55 : use thread_mod, only : horz_num_threads
56 : use dimensions_mod, only : fv_nphys
57 : implicit none
58 : type (element_t), intent(inout), dimension(:) :: elem
59 : integer, intent(in) :: tl, nphys, tlq
60 : real (kind=r8), intent(out) :: frontgf(nphys*nphys,pver,nelemd)
61 : real (kind=r8), intent(out) :: frontga(nphys*nphys,pver,nelemd)
62 :
63 : ! Local variables
64 : type (hybrid_t) :: hybrid
65 : integer :: nets, nete, ithr, ncols, ie
66 370944 : real(kind=r8), allocatable :: frontgf_thr(:,:,:,:)
67 370944 : real(kind=r8), allocatable :: frontga_thr(:,:,:,:)
68 :
69 : ! This does not need to be a thread private data-structure
70 370944 : call derivinit(deriv)
71 : !!$OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(nets,nete,hybrid,ie,ncols,frontgf_thr,frontga_thr)
72 : ! hybrid = config_thread_region(par,'horizontal')
73 370944 : hybrid = config_thread_region(par,'serial')
74 370944 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
75 :
76 2225664 : allocate(frontgf_thr(nphys,nphys,nlev,nets:nete))
77 1483776 : allocate(frontga_thr(nphys,nphys,nlev,nets:nete))
78 370944 : call compute_frontogenesis(frontgf_thr,frontga_thr,tl,tlq,elem,deriv,hybrid,nets,nete,nphys)
79 370944 : if (fv_nphys>0) then
80 2979144 : do ie=nets,nete
81 7824600 : frontgf(:,:,ie) = RESHAPE(frontgf_thr(:,:,:,ie),(/nphys*nphys,nlev/))
82 8195544 : frontga(:,:,ie) = RESHAPE(frontga_thr(:,:,:,ie),(/nphys*nphys,nlev/))
83 : end do
84 : else
85 0 : do ie=nets,nete
86 0 : ncols = elem(ie)%idxP%NumUniquePts
87 0 : call UniquePoints(elem(ie)%idxP, nlev, frontgf_thr(:,:,:,ie), frontgf(1:ncols,:,ie))
88 0 : call UniquePoints(elem(ie)%idxP, nlev, frontga_thr(:,:,:,ie), frontga(1:ncols,:,ie))
89 : end do
90 : end if
91 370944 : deallocate(frontga_thr)
92 370944 : deallocate(frontgf_thr)
93 : !!$OMP END PARALLEL
94 :
95 741888 : end subroutine gws_src_fnct
96 :
97 370944 : subroutine compute_frontogenesis(frontgf,frontga,tl,tlq,elem,ederiv,hybrid,nets,nete,nphys)
98 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
99 : ! compute frontogenesis function F
100 : ! F = -gradth dot C
101 : ! with:
102 : ! theta = potential temperature
103 : ! gradth = grad(theta)
104 : ! C = ( gradth dot grad ) U
105 : !
106 : ! Original by Mark Taylor, July 2011
107 : ! Change by Santos, 10 Aug 2011:
108 : ! Integrated into gravity_waves_sources module, several arguments made global
109 : ! to prevent repeated allocation/initialization
110 : !
111 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
112 370944 : use physconst, only: cappa
113 : use air_composition,only: dry_air_species_num, thermodynamic_active_species_num
114 : use air_composition,only: thermodynamic_active_species_idx_dycore
115 : use derivative_mod, only: gradient_sphere, ugradv_sphere
116 : use edge_mod, only: edgevpack, edgevunpack
117 : use bndry_mod, only: bndry_exchange
118 : use dyn_grid, only: hvcoord
119 : use dimensions_mod, only: fv_nphys,ntrac
120 : use fvm_mapping, only: dyn2phys_vector,dyn2phys
121 :
122 : type(hybrid_t), intent(in) :: hybrid
123 : type(element_t), intent(inout), target :: elem(:)
124 : type(derivative_t), intent(in) :: ederiv
125 : integer, intent(in) :: nets,nete,nphys
126 : integer, intent(in) :: tl,tlq
127 : real(r8), intent(out) :: frontgf(nphys,nphys,nlev,nets:nete)
128 : real(r8), intent(out) :: frontga(nphys,nphys,nlev,nets:nete)
129 :
130 : ! local
131 741888 : real(r8) :: area_inv(fv_nphys,fv_nphys), tmp(np,np)
132 741888 : real(r8) :: uv_tmp(fv_nphys*fv_nphys,2,nlev)
133 741888 : real(r8) :: frontgf_gll(np,np,nlev,nets:nete)
134 : real(r8) :: frontga_gll(np,np,nlev,nets:nete)
135 : integer :: k,kptr,i,j,ie,component,h,nq,m_cnst
136 741888 : real(r8) :: gradth(np,np,2,nlev,nets:nete) ! grad(theta)
137 : real(r8) :: p(np,np) ! pressure at mid points
138 : real(r8) :: pint(np,np) ! pressure at interface points
139 : real(r8) :: theta(np,np) ! potential temperature at mid points
140 : real(r8) :: C(np,np,2), sum_water(np,np)
141 :
142 2979144 : do ie=nets,nete
143 : ! pressure at model top
144 54772200 : pint(:,:) = hvcoord%hyai(1)*hvcoord%ps0
145 245170800 : do k=1,nlev
146 : ! moist pressure at mid points
147 5093814600 : sum_water(:,:) = 1.0_r8
148 1697938200 : do nq=dry_air_species_num+1,thermodynamic_active_species_num
149 1455375600 : m_cnst = thermodynamic_active_species_idx_dycore(nq)
150 : !
151 : ! make sure Q is updated
152 : !
153 30805450200 : sum_water(:,:) = sum_water(:,:) + elem(ie)%state%Qdp(:,:,k,m_cnst,tlq)/elem(ie)%state%dp3d(:,:,k,tl)
154 : end do
155 5093814600 : p(:,:) = pint(:,:) + 0.5_r8*sum_water(:,:)*elem(ie)%state%dp3d(:,:,k,tl)
156 : ! moist pressure at interface for next iteration
157 5093814600 : pint(:,:) = pint(:,:)+elem(ie)%state%dp3d(:,:,k,tl)
158 : !
159 5093814600 : theta(:,:) = elem(ie)%state%T(:,:,k,tl)*(psurf_ref / p(:,:))**cappa
160 : ! gradth(:,:,:,k,ie) = gradient_sphere(theta,ederiv,elem(ie)%Dinv)
161 242562600 : call gradient_sphere(theta,ederiv,elem(ie)%Dinv,gradth(:,:,:,k,ie))
162 : ! compute C = (grad(theta) dot grad ) u
163 242562600 : C(:,:,:) = ugradv_sphere(gradth(:,:,:,k,ie), elem(ie)%state%v(:,:,:,k,tl),ederiv,elem(ie))
164 : ! gradth dot C
165 5093814600 : frontgf_gll(:,:,k,ie) = -( C(:,:,1)*gradth(:,:,1,k,ie) + C(:,:,2)*gradth(:,:,2,k,ie) )
166 : ! apply mass matrix
167 5093814600 : gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%spheremp(:,:)
168 5093814600 : gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%spheremp(:,:)
169 5096422800 : frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%spheremp(:,:)
170 : enddo
171 : ! pack
172 2608200 : call edgeVpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
173 2979144 : call edgeVpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
174 : enddo
175 370944 : call bndry_exchange(hybrid,edge3,location='compute_frontogenesis')
176 2979144 : do ie=nets,nete
177 2608200 : call edgeVunpack(edge3, frontgf_gll(:,:,:,ie),nlev,0,ie)
178 2608200 : call edgeVunpack(edge3, gradth(:,:,:,:,ie),2*nlev,nlev,ie)
179 : ! apply inverse mass matrix,
180 245170800 : do k=1,nlev
181 5093814600 : gradth(:,:,1,k,ie)=gradth(:,:,1,k,ie)*elem(ie)%rspheremp(:,:)
182 5093814600 : gradth(:,:,2,k,ie)=gradth(:,:,2,k,ie)*elem(ie)%rspheremp(:,:)
183 5096422800 : frontgf_gll(:,:,k,ie)=frontgf_gll(:,:,k,ie)*elem(ie)%rspheremp(:,:)
184 : end do
185 2979144 : if (fv_nphys>0) then
186 2608200 : uv_tmp(:,:,:) = dyn2phys_vector(gradth(:,:,:,:,ie),elem(ie))
187 245170800 : do k=1,nlev
188 242562600 : h=0
189 972858600 : do j=1,fv_nphys
190 3153313800 : do i=1,fv_nphys
191 2183063400 : h=h+1
192 2910751200 : frontga(i,j,k,ie) = atan2 ( uv_tmp(h,2,k) , uv_tmp(h,1,k) + 1.e-10_r8 )
193 : end do
194 : end do
195 : end do
196 : !
197 : ! compute inverse physgrid area for mapping of scaler
198 : !
199 54772200 : tmp = 1.0_r8
200 2608200 : area_inv = dyn2phys(tmp,elem(ie)%metdet)
201 33906600 : area_inv = 1.0_r8/area_inv
202 245170800 : do k=1,nlev
203 245170800 : frontgf(:,:,k,ie) = dyn2phys(frontgf_gll(:,:,k,ie),elem(ie)%metdet,area_inv)
204 : end do
205 : else
206 0 : do k=1,nlev
207 0 : frontgf(:,:,k,ie)=frontgf_gll(:,:,k,ie)
208 : ! Frontogenesis angle
209 0 : frontga(:,:,k,ie) = atan2 ( gradth(:,:,2,k,ie) , gradth(:,:,1,k,ie) + 1.e-10_r8 )
210 : end do
211 : end if
212 : enddo
213 370944 : end subroutine compute_frontogenesis
214 :
215 :
216 : end module gravity_waves_sources
|
