Line data Source code
1 : #ifdef HAVE_CONFIG_H
2 : #include "config.h"
3 : #endif
4 :
5 : !MODULE FVM_CONTROL_VOLUME_MOD---------------------------------------------CE-for FVM
6 : ! AUTHOR: Christoph Erath, 11.June 2011 !
7 : ! This module contains everything to initialize the arrival. It also provides the !
8 : ! interpolation points for the reconstruction (projection from one face to another !
9 : ! when the element is on the cube edge) !
10 : ! It also intialize the start values, see also fvm_analytic !
11 : !-----------------------------------------------------------------------------------!
12 : module fvm_control_volume_mod
13 : use shr_kind_mod, only: r8=>shr_kind_r8
14 : use coordinate_systems_mod, only: spherical_polar_t
15 : use element_mod, only: element_t
16 : use dimensions_mod, only: nc, nhe, nlev, ntrac_d, qsize_d,ne, np, nhr, ns, nhc
17 : use dimensions_mod, only: fv_nphys, nhe_phys, nhr_phys, ns_phys, nhc_phys,fv_nphys
18 : use dimensions_mod, only: irecons_tracer
19 : use cam_abortutils, only: endrun
20 :
21 : implicit none
22 : private
23 : integer, parameter, private:: nh = nhr+(nhe-1) ! = 2 (nhr=2; nhe=1)
24 : ! = 3 (nhr=2; nhe=2)
25 :
26 : type, public :: fvm_struct
27 : ! fvm tracer mixing ratio: (kg/kg)
28 : real (kind=r8) :: c(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev,ntrac_d)
29 : real (kind=r8) :: se_flux(1-nhe:nc+nhe,1-nhe:nc+nhe,4,nlev)
30 :
31 : real (kind=r8) :: dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev)
32 : real (kind=r8) :: dp_ref(nlev)
33 : real (kind=r8) :: dp_ref_inverse(nlev)
34 : real (kind=r8) :: psc(nc,nc)
35 :
36 : real (kind=r8) :: inv_area_sphere(nc,nc) ! inverse area_sphere
37 : real (kind=r8) :: inv_se_area_sphere(nc,nc) ! inverse area_sphere
38 :
39 : integer :: faceno !face number
40 : ! number of south,....,swest and 0 for interior element
41 : integer :: cubeboundary
42 :
43 : #ifdef waccm_debug
44 : real (kind=r8) :: CSLAM_gamma(nc,nc,nlev,4)
45 : #endif
46 : real (kind=r8) :: displ_max(1-nhc:nc+nhc,1-nhc:nc+nhc,4)
47 : integer :: flux_vec (2,1-nhc:nc+nhc,1-nhc:nc+nhc,4)
48 : !
49 : !
50 : ! cartesian location of vertices for flux sides
51 : !
52 : ! x-coordinate of vertex 1: vtx_cart(1,1i,j,1,1) = fvm%acartx(i)
53 : ! y-coordinate of vertex 1: vtx_cart(1,2,i,j,2,1) = fvm%acarty(j)
54 : !
55 : ! x-coordinate of vertex 2: vtx_cart(2,1,i,j) = fvm%acartx(i+1)
56 : ! y-coordinate of vertex 2: vtx_cart(2,2,i,j) = fvm%acarty(j )
57 : !
58 : ! x-coordinate of vertex 3: vtx_cart(3,1,i,j) = fvm%acartx(i+1)
59 : ! y-coordinate of vertex 3: vtx_cart(3,2,i,j) = fvm%acarty(j+1)
60 : !
61 : ! x-coordinate of vertex 4: vtx_cart(4,1,i,j) = fvm%acartx(i )
62 : ! y-coordinate of vertex 4: vtx_cart(4,2,i,j) = fvm%acarty(j+1)
63 : !
64 : real (kind=r8) :: vtx_cart (4,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
65 : !
66 : ! flux_orient(1,i,j) = panel on which control volume (i,j) is located
67 : ! flux_orient(2,i,j) = cshift value for vertex permutation
68 : !
69 : real (kind=r8) :: flux_orient(2 ,1-nhc:nc+nhc,1-nhc:nc+nhc)
70 : !
71 : ! i,j: indicator function for non-existent cells (0 for corner halo and 1 elsewhere)
72 : !
73 : integer :: ifct (1-nhc:nc+nhc,1-nhc:nc+nhc)
74 : integer :: rot_matrix(2,2,1-nhc:nc+nhc,1-nhc:nc+nhc)
75 : !
76 : real (kind=r8) :: dalpha, dbeta ! central-angle for gnomonic coordinates
77 : type (spherical_polar_t) :: center_cart(nc,nc) ! center of fvm cell in gnomonic coordinates
78 : real (kind=r8) :: area_sphere(nc,nc) ! spherical area of fvm cell
79 : real (kind=r8) :: spherecentroid(irecons_tracer-1,1-nhc:nc+nhc,1-nhc:nc+nhc) ! centroids
80 : !
81 : ! pre-computed metric terms (for efficiency)
82 : !
83 : ! recons_metrics(1,:,:) = spherecentroid(1,:,:)**2 -spherecentroid(3,:,:)
84 : ! recons_metrics(2,:,:) = spherecentroid(2,:,:)**2 -spherecentroid(4,:,:)
85 : ! recons_metrics(3,:,:) = spherecentroid(1,:,:)*spherecentroid(2,:,:)-spherecentroid(5,:,:)
86 : !
87 : real (kind=r8) :: recons_metrics(3,1-nhe:nc+nhe,1-nhe:nc+nhe)
88 : !
89 : ! recons_metrics_integral(1,:,:) = 2.0_r8*spherecentroid(1,:,:)**2 -spherecentroid(3,:,:)
90 : ! recons_metrics_integral(2,:,:) = 2.0_r8*spherecentroid(2,:,:)**2 -spherecentroid(4,:,:)
91 : ! recons_metrics_integral(3,:,:) = 2.0_r8*spherecentroid(1,:,:)*spherecentroid(2,:,:)-spherecentroid(5,:,:)
92 : !
93 : real (kind=r8) :: recons_metrics_integral(3,1-nhe:nc+nhe,1-nhe:nc+nhe)
94 : !
95 : integer :: jx_min(3), jx_max(3), jy_min(3), jy_max(3) !bounds for computation
96 :
97 : ! provide fixed interpolation points with respect to the arrival grid for
98 : ! reconstruction
99 : integer :: ibase(1-nh:nc+nh,1:nhr,2)
100 : real (kind=r8) :: halo_interp_weight(1:ns,1-nh:nc+nh,1:nhr,2)
101 : real (kind=r8) :: centroid_stretch(7,1-nhe:nc+nhe,1-nhe:nc+nhe) !for finite-difference reconstruction
102 : !
103 : ! pre-compute weights for reconstruction at cell vertices
104 : !
105 : ! ! Evaluate constant order terms
106 : ! value = fcube(a,b) + &
107 : ! ! Evaluate linear order terms
108 : ! recons(1,a,b) * (cartx - centroid(1,a,b)) + &
109 : ! recons(2,a,b) * (carty - centroid(2,a,b)) + &
110 : ! ! Evaluate second order terms
111 : ! recons(3,a,b) * (centroid(1,a,b)**2 - centroid(3,a,b)) + &
112 : ! recons(4,a,b) * (centroid(2,a,b)**2 - centroid(4,a,b)) + &
113 : ! recons(5,a,b) * (centroid(1,a,b) * centroid(2,a,b) - centroid(5,a,b)) + &
114 : !
115 : ! recons(3,a,b) * (cartx - centroid(1,a,b))**2 + &
116 : ! recons(4,a,b) * (carty - centroid(2,a,b))**2 + &
117 : ! recons(5,a,b) * (cartx - centroid(1,a,b)) * (carty - centroid(2,a,b))
118 : !
119 : real (kind=r8) :: vertex_recons_weights(4,1:irecons_tracer-1,1-nhe:nc+nhe,1-nhe:nc+nhe)
120 : !
121 : ! for mapping fvm2dyn
122 : !
123 : real (kind=r8) :: norm_elem_coord(2,1-nhc:nc+nhc,1-nhc:nc+nhc)
124 : !
125 : !******************************************
126 : !
127 : ! separate physics grid variables
128 : !
129 : !******************************************
130 : !
131 : real (kind=r8) , allocatable :: vtx_cart_physgrid(:,:,:,:)
132 : real (kind=r8) , allocatable :: flux_orient_physgrid(:,:,:)
133 : integer , allocatable :: ifct_physgrid(:,:)
134 : integer , allocatable :: rot_matrix_physgrid(:,:,:,:)
135 : real (kind=r8) , allocatable :: spherecentroid_physgrid(:,:,:)
136 : real (kind=r8) , allocatable :: recons_metrics_physgrid(:,:,:)
137 : real (kind=r8) , allocatable :: recons_metrics_integral_physgrid(:,:,:)
138 : ! centroid_stretch_physgrid for finite-difference reconstruction
139 : real (kind=r8) , allocatable :: centroid_stretch_physgrid (:,:,:)
140 : real (kind=r8) :: dalpha_physgrid, dbeta_physgrid ! central-angle for gnomonic coordinates
141 : type (spherical_polar_t) , allocatable :: center_cart_physgrid(:,:) ! center of fvm cell in gnomonic coordinates
142 : real (kind=r8) , allocatable :: area_sphere_physgrid(:,:) ! spherical area of fvm cell
143 : integer :: jx_min_physgrid(3), jx_max_physgrid(3) !bounds for computation
144 : integer :: jy_min_physgrid(3), jy_max_physgrid(3) !bounds for computation
145 : integer , allocatable :: ibase_physgrid(:,:,:)
146 : real (kind=r8) , allocatable :: halo_interp_weight_physgrid(:,:,:,:)
147 : real (kind=r8) , allocatable :: vertex_recons_weights_physgrid(:,:,:,:)
148 :
149 : real (kind=r8) , allocatable :: norm_elem_coord_physgrid(:,:,:)
150 : real (kind=r8) , allocatable :: Dinv_physgrid(:,:,:,:)
151 :
152 : real (kind=r8) , allocatable :: fc(:,:,:,:)
153 : real (kind=r8) , allocatable :: fc_phys(:,:,:,:)
154 : real (kind=r8) , allocatable :: ft(:,:,:)
155 : real (kind=r8) , allocatable :: fm(:,:,:,:)
156 : real (kind=r8) , allocatable :: dp_phys(:,:,:)
157 : end type fvm_struct
158 :
159 : public :: fvm_mesh, fvm_set_cubeboundary, allocate_physgrid_vars
160 :
161 :
162 : real (kind=r8),parameter, public :: bignum = 1.0E20_r8
163 :
164 : contains
165 10800 : subroutine fvm_set_cubeboundary(elem, fvm)
166 : implicit none
167 : type (element_t) , intent(in) :: elem
168 : type (fvm_struct), intent(inout) :: fvm
169 :
170 : logical :: corner
171 : integer :: j, mynbr_cnt, mystart
172 : integer :: nbrsface(8)! store the neighbours in north, south
173 :
174 10800 : fvm%faceno=elem%FaceNum
175 : ! write the neighbors in the structure
176 10800 : fvm%cubeboundary=0
177 10800 : corner=.FALSE.
178 97200 : do j=1,8
179 86400 : mynbr_cnt = elem%vertex%nbrs_ptr(j+1) - elem%vertex%nbrs_ptr(j) !length of neighbor location
180 86400 : mystart = elem%vertex%nbrs_ptr(j)
181 : !NOTE: assuming that we do not have multiple corner neighbors (so not a refined mesh)
182 97200 : if (mynbr_cnt > 0 ) then
183 86352 : nbrsface(j)=elem%vertex%nbrs_face(mystart)
184 : ! note that if the element lies on a corner, it will be at j=5,6,7,8
185 86352 : if ((nbrsface(j) /= fvm%faceno) .AND. (j<5)) then
186 1440 : fvm%cubeboundary=j
187 : endif
188 : else ! corner on the cube
189 48 : if (.NOT. corner) then
190 48 : nbrsface(j)=-1
191 48 : fvm%cubeboundary=j
192 48 : corner=.TRUE.
193 : else
194 0 : if ( ne == 0 ) then
195 : ! dont check this condition. note that we call this code
196 : ! generate phys grid template files, so we need to be able
197 : ! to call create_ari() to create the subcells even though
198 : ! cslam cant run with the unstructed ne=0 case
199 : else
200 0 : print *,'Error in fvm_CONTROL_VOLUME_MOD - Subroutine fvm_MESH_ARI: '
201 0 : call endrun('Do not allow one element per face for fvm, please increase ne!')
202 : endif
203 : endif
204 : end if
205 : end do
206 10800 : end subroutine fvm_set_cubeboundary
207 :
208 10800 : subroutine fvm_mesh(elem, fvm)
209 : use fvm_analytic_mod, only : compute_halo_vars
210 : use fvm_analytic_mod, only : create_interpolation_points
211 : use derivative_mod , only : subcell_integration
212 :
213 : implicit none
214 : type (element_t), intent(in) :: elem
215 : type (fvm_struct), intent(inout) :: fvm
216 : integer :: i,j
217 : real (kind=r8) :: tmp(np,np)
218 : !
219 : ! initialize metric and related terms on panel
220 : !
221 : call compute_halo_vars(& !input
222 : fvm%faceno,fvm%cubeboundary,nc,nhc,nhe, & !input
223 : fvm%jx_min,fvm%jx_max,fvm%jy_min,fvm%jy_max,&!output
224 10800 : fvm%flux_orient,fvm%ifct,fvm%rot_matrix) !output
225 43200 : do j=1,nc
226 140400 : do i=1,nc
227 97200 : fvm%norm_elem_coord(1,i,j) = elem%corners(1)%x+(i-0.5_r8)*fvm%dalpha
228 129600 : fvm%norm_elem_coord(2,i,j) = elem%corners(1)%y+(j-0.5_r8)*fvm%dalpha
229 : end do
230 : end do
231 :
232 : !
233 : ! overwrite areas for consistency with SE areas (that are O(10E-5) incorrect)
234 : !
235 : ! tmp = 1.0_r8
236 : ! call subcell_integration(tmp, np, nc, elem%metdet,fvm%area_sphere)
237 : !
238 : ! do the same for physics grid
239 : !
240 : call compute_halo_vars(&
241 : fvm%faceno,fvm%cubeboundary,fv_nphys,nhc_phys,nhe_phys,&
242 : fvm%jx_min_physgrid,fvm%jx_max_physgrid,fvm%jy_min_physgrid,fvm%jy_max_physgrid,&
243 10800 : fvm%flux_orient_physgrid,fvm%ifct_physgrid,fvm%rot_matrix_physgrid)
244 43200 : do j=1,fv_nphys
245 140400 : do i=1,fv_nphys
246 97200 : fvm%norm_elem_coord_physgrid(1,i,j) = elem%corners(1)%x+(i-0.5_r8)*fvm%dalpha_physgrid
247 129600 : fvm%norm_elem_coord_physgrid(2,i,j) = elem%corners(1)%y+(j-0.5_r8)*fvm%dalpha_physgrid
248 : end do
249 : end do
250 : !
251 : ! initialize halo interpolation variables
252 : !
253 : call create_interpolation_points(elem,&
254 : nc,nhc,nhr,ns,nh,fvm%cubeboundary,&
255 10800 : fvm%dalpha,fvm%dbeta,fvm%ibase,fvm%halo_interp_weight)
256 : call create_interpolation_points(elem,&
257 : fv_nphys,nhc_phys,nhr_phys,ns_phys,nhr_phys,fvm%cubeboundary,&
258 10800 : fvm%dalpha_physgrid,fvm%dbeta_physgrid,fvm%ibase_physgrid,fvm%halo_interp_weight_physgrid)
259 10800 : end subroutine fvm_mesh
260 :
261 :
262 1536 : subroutine allocate_physgrid_vars(fvm,par)
263 10800 : use cam_logfile , only : iulog
264 : use parallel_mod , only : parallel_t
265 : use dimensions_mod, only : nelemd
266 : type (fvm_struct), intent(inout) :: fvm(:)
267 : type (parallel_t), intent(in) :: par
268 : integer :: ie
269 :
270 1536 : nhc_phys = fv_nphys
271 1536 : nhe_phys = 0
272 1536 : nhr_phys = 2
273 1536 : ns_phys = MAX(fv_nphys,2)
274 :
275 1536 : if(par%masterproc) then
276 2 : write(iulog,*)"allocating physgrid grid vars"
277 2 : write(iulog,*)"fv_nphys,nhc_phys,nhe_phys,nhr_phys,ns_phys = ",&
278 4 : fv_nphys,nhc_phys,nhe_phys,nhr_phys,ns_phys
279 : end if
280 :
281 12336 : do ie=1,nelemd
282 43200 : allocate(fvm(ie)%vtx_cart_physgrid (4,2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
283 43200 : allocate(fvm(ie)%flux_orient_physgrid (2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
284 43200 : allocate(fvm(ie)%ifct_physgrid (1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
285 43200 : allocate(fvm(ie)%rot_matrix_physgrid (2,2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
286 :
287 0 : allocate(fvm(ie)%spherecentroid_physgrid(irecons_tracer-1,&
288 43200 : 1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
289 43200 : allocate(fvm(ie)%recons_metrics_physgrid (3,1-nhe_phys:fv_nphys+nhe_phys,1-nhe_phys:fv_nphys+nhe_phys))
290 32400 : allocate(fvm(ie)%recons_metrics_integral_physgrid(3,1-nhe_phys:fv_nphys+nhe_phys,1-nhe_phys:fv_nphys+nhe_phys))
291 43200 : allocate(fvm(ie)%centroid_stretch_physgrid (7,1-nhe_phys:fv_nphys+nhe_phys,1-nhe_phys:fv_nphys+nhe_phys))
292 43200 : allocate(fvm(ie)%center_cart_physgrid(fv_nphys,fv_nphys))
293 43200 : allocate(fvm(ie)%area_sphere_physgrid(fv_nphys,fv_nphys))
294 54000 : allocate(fvm(ie)%ibase_physgrid(1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,2))
295 64800 : allocate(fvm(ie)%halo_interp_weight_physgrid(1:ns_phys,1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,2))
296 0 : allocate(fvm(ie)%vertex_recons_weights_physgrid(4,1:irecons_tracer-1,1-nhe_phys:fv_nphys+nhe_phys,&
297 43200 : 1-nhe_phys:fv_nphys+nhe_phys))
298 :
299 32400 : allocate(fvm(ie)%norm_elem_coord_physgrid(2,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys ))
300 64800 : allocate(fvm(ie)%Dinv_physgrid ( 1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,2,2))
301 :
302 10800 : allocate(fvm(ie)%fc(nc,nc,nlev,max(ntrac_d,qsize_d)))
303 64800 : allocate(fvm(ie)%fc_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev,max(ntrac_d,qsize_d)))
304 43200 : allocate(fvm(ie)%ft(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev))
305 54000 : allocate(fvm(ie)%fm(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,2,nlev))
306 33936 : allocate(fvm(ie)%dp_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev))
307 : end do
308 1536 : end subroutine allocate_physgrid_vars
309 0 : end module fvm_control_volume_mod
|
