Line data Source code
1 : #define FVM_TIMERS .FALSE.
2 : !-----------------------------------------------------------------------------!
3 : !MODULE FVM_MOD-----------------------------------------------------CE-for FVM!
4 : ! FVM_MOD File for the fvm project in HOMME !
5 : ! Author: Christoph Erath !
6 : ! Date: 25.January 2011 !
7 : ! MAIN module to run fvm on HOMME !
8 : ! 14.November 2011: reorganisation done !
9 : ! 7.Februar 2012: cslam_run and cslam_runair !
10 : !-----------------------------------------------------------------------------!
11 :
12 : module fvm_mod
13 : use shr_kind_mod, only: r8=>shr_kind_r8
14 : use edge_mod, only: initghostbuffer, freeghostbuffer, ghostpack, ghostunpack
15 : use edgetype_mod, only: edgebuffer_t
16 : use bndry_mod, only: ghost_exchange
17 : use thread_mod, only: horz_num_threads, vert_num_threads
18 :
19 : use element_mod, only: element_t
20 : use fvm_control_volume_mod, only: fvm_struct
21 : use hybrid_mod, only: hybrid_t
22 :
23 : implicit none
24 : private
25 : save
26 :
27 : type (EdgeBuffer_t) :: edgeveloc
28 : type (EdgeBuffer_t), public :: ghostBufQnhc_s
29 : type (EdgeBuffer_t), public :: ghostBufQnhc_t1
30 : type (EdgeBuffer_t), public :: ghostBufQnhc_vh
31 : type (EdgeBuffer_t), public :: ghostBufQnhc_h
32 : type (EdgeBuffer_t), public :: ghostBufQ1_h
33 : type (EdgeBuffer_t), public :: ghostBufQ1_vh
34 : ! type (EdgeBuffer_t), private :: ghostBufFlux_h
35 : type (EdgeBuffer_t), public :: ghostBufFlux_vh
36 : type (EdgeBuffer_t), public :: ghostBufQnhcJet_h
37 : type (EdgeBuffer_t), public :: ghostBufFluxJet_h
38 : type (EdgeBuffer_t), public :: ghostBufPG_s
39 : type (EdgeBuffer_t), public :: ghostBuf_cslam2gll
40 :
41 : interface fill_halo_fvm
42 : module procedure fill_halo_fvm_noprealloc
43 : module procedure fill_halo_fvm_prealloc
44 : end interface
45 :
46 :
47 : public :: edgeveloc, fvm_init1,fvm_init2, fill_halo_fvm, fvm_pg_init,fvm_init3,fill_halo_and_extend_panel
48 :
49 : contains
50 :
51 0 : subroutine fill_halo_fvm_noprealloc(elem,fvm,hybrid,nets,nete,ndepth,kmin,kmax,ksize)
52 : use perf_mod, only : t_startf, t_stopf ! _EXTERNAL
53 : use dimensions_mod, only: nc, ntrac, nlev
54 : implicit none
55 : type (element_t),intent(inout) :: elem(:)
56 : type (fvm_struct),intent(inout) :: fvm(:)
57 : type (hybrid_t),intent(in) :: hybrid
58 :
59 0 : type (edgeBuffer_t) :: cellghostbuf
60 :
61 : integer,intent(in) :: nets,nete
62 : integer,intent(in) :: ndepth ! depth of halo
63 : integer,intent(in) :: kmin,kmax ! min and max vertical level
64 : integer,intent(in) :: ksize ! the total number of vertical
65 :
66 : integer :: ie,i1,i2,kblk,kptr,q
67 : !
68 : !
69 :
70 0 : if(kmin .ne. 1 .or. kmax .ne. nlev) then
71 0 : print *,'WARNING: fill_halo_fvm_noprealloc does not support the passing of non-contigous arrays'
72 0 : print *,'WARNING: incorrect answers are likely'
73 : endif
74 : if(FVM_TIMERS) call t_startf('FVM:initbuf')
75 0 : i1=1-ndepth
76 0 : i2=nc+ndepth
77 0 : kblk = kmax-kmin+1
78 0 : call initghostbuffer(hybrid%par,cellghostbuf,elem,kblk*(ntrac+1),ndepth,nc)
79 : if(FVM_TIMERS) call t_stopf('FVM:initbuf')
80 : if(FVM_TIMERS) call t_startf('FVM:pack')
81 0 : do ie=nets,nete
82 0 : kptr = kmin-1
83 0 : call ghostpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk, kptr,ie)
84 0 : do q=1,ntrac
85 0 : kptr = kptr + ksize
86 0 : call ghostpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,q) ,kblk,kptr,ie)
87 : enddo
88 : end do
89 : if(FVM_TIMERS) call t_stopf('FVM:pack')
90 : if(FVM_TIMERS) call t_startf('FVM:Communication')
91 0 : call ghost_exchange(hybrid,cellghostbuf,location='fill_halo_fvm_noprealloc')
92 : if(FVM_TIMERS) call t_stopf('FVM:Communication')
93 : !-----------------------------------------------------------------------------------!
94 : if(FVM_TIMERS) call t_startf('FVM:Unpack')
95 0 : do ie=nets,nete
96 0 : kptr = kmin-1
97 0 : call ghostunpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk ,kptr,ie)
98 0 : do q=1,ntrac
99 0 : kptr = kptr + ksize
100 0 : call ghostunpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,:), kblk,kptr,ie)
101 : enddo
102 : enddo
103 : if(FVM_TIMERS) call t_stopf('FVM:Unpack')
104 : if(FVM_TIMERS) call t_startf('FVM:freebuf')
105 0 : call freeghostbuffer(cellghostbuf)
106 : if(FVM_TIMERS) call t_stopf('FVM:freebuf')
107 0 : end subroutine fill_halo_fvm_noprealloc
108 :
109 738816 : subroutine fill_halo_fvm_prealloc(cellghostbuf,elem,fvm,hybrid,nets,nete,ndepth,kmin,kmax,ksize,active)
110 0 : use perf_mod, only : t_startf, t_stopf ! _EXTERNAL
111 : use dimensions_mod, only: nc, ntrac, nlev
112 : implicit none
113 : type (EdgeBuffer_t), intent(inout) :: cellghostbuf
114 : type (element_t),intent(inout) :: elem(:)
115 : type (fvm_struct),intent(inout) :: fvm(:)
116 : type (hybrid_t),intent(in) :: hybrid
117 :
118 :
119 : integer,intent(in) :: nets,nete
120 : integer,intent(in) :: ndepth ! depth of halo
121 : integer,intent(in) :: kmin,kmax ! min and max vertical level
122 : integer,intent(in) :: ksize ! the total number of vertical
123 : logical, optional :: active ! indicates if te current thread is active
124 : integer :: ie,i1,i2,kblk,q,kptr
125 : !
126 : !
127 : logical :: lactive
128 :
129 738816 : if(present(active)) then
130 738816 : lactive = active
131 : else
132 : lactive = .true.
133 : endif
134 : ! call t_startf('FVM:initbuf')
135 738816 : i1=1-ndepth
136 738816 : i2=nc+ndepth
137 738816 : kblk = kmax-kmin+1
138 : if(FVM_TIMERS) call t_startf('FVM:pack')
139 738816 : if(lactive) then
140 5933616 : do ie=nets,nete
141 5194800 : kptr = kmin-1
142 4192203600 : call ghostpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk, kptr,ie)
143 21518016 : do q=1, ntrac
144 15584400 : kptr = kptr + ksize
145 12581805600 : call ghostpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,q) ,kblk,kptr,ie)
146 : enddo
147 : end do
148 : endif
149 : if(FVM_TIMERS) call t_stopf('FVM:pack')
150 : if(FVM_TIMERS) call t_startf('FVM:Communication')
151 738816 : call ghost_exchange(hybrid,cellghostbuf,location='fill_halo_fvm_prealloc')
152 : if(FVM_TIMERS) call t_stopf('FVM:Communication')
153 : !-----------------------------------------------------------------------------------!
154 : if(FVM_TIMERS) call t_startf('FVM:Unpack')
155 738816 : if(lactive) then
156 5933616 : do ie=nets,nete
157 5194800 : kptr = kmin-1
158 8379212400 : call ghostunpack(cellghostbuf, fvm(ie)%dp_fvm(i1:i2,i1:i2,kmin:kmax),kblk, kptr,ie)
159 21518016 : do q=1, ntrac
160 15584400 : kptr = kptr + ksize
161 25142832000 : call ghostunpack(cellghostbuf, fvm(ie)%c(i1:i2,i1:i2,kmin:kmax,q), kblk,kptr,ie)
162 : enddo
163 : enddo
164 : endif
165 : if(FVM_TIMERS) call t_stopf('FVM:Unpack')
166 :
167 738816 : end subroutine fill_halo_fvm_prealloc
168 :
169 : subroutine PrintArray(i1,i2,array)
170 : ! debug routine potentially called from any MPI rank
171 : integer :: i1,i2
172 : real(kind=r8) :: array(i1:i2,i1:i2)
173 : integer :: sz,i,ub
174 :
175 : sz = size(array,dim=1)
176 :
177 : if (sz == 9) then
178 : do i=i2,i1,-1
179 : write(6,9) array(-2,i),array(-1,i), array(0,i), &
180 : array( 1,i), array(2,i), array(3,i), &
181 : array( 4,i), array(5,i), array(6,i)
182 : enddo
183 : endif
184 :
185 : 9 format('|',9(f10.1,'|'))
186 :
187 :
188 738816 : end subroutine
189 :
190 :
191 0 : subroutine fill_halo_and_extend_panel(elem,fvm,fld,hybrid,nets,nete,nphys,nhcc, ndepth,numlev,num_flds,lfill_halo,lextend_panel)
192 : use hybrid_mod, only: hybrid_t
193 : use edge_mod, only: initghostbuffer, freeghostbuffer, ghostpack, ghostunpack
194 :
195 : use fvm_reconstruction_mod, only: extend_panel_interpolate
196 : use cam_abortutils, only: endrun
197 : use dimensions_mod, only: fv_nphys,nhr,nhr_phys,nhc,nhc_phys,ns,ns_phys,nhe_phys,nc
198 : use perf_mod, only : t_startf, t_stopf ! _EXTERNAL
199 :
200 : integer , intent(in) :: nets,nete,nphys,ndepth,numlev,num_flds,nhcc
201 : real (kind=r8) , intent(inout) :: fld(1-nhcc:nphys+nhcc,1-nhcc:nphys+nhcc,numlev,num_flds,nets:nete)
202 : type (hybrid_t) , intent(in) :: hybrid ! distributed parallel structure (shared)
203 : type (element_t) , intent(inout) :: elem(:)
204 : type(fvm_struct) , intent(in) :: fvm(:)
205 : logical , intent(in) :: lfill_halo,lextend_panel
206 : ! real (kind=r8) , allocatable :: ftmp(:,:)
207 : ! real (kind=r8) :: ftmp(1-nhcc:nphys+nhcc,1-nhcc:nphys+nhcc,numlev,num_flds,nets:nete)
208 0 : real (kind=r8), allocatable :: fld_tmp(:,:)
209 :
210 : integer :: ie,k,itr,nht_phys,nh_phys
211 0 : type (edgeBuffer_t) :: cellghostbuf
212 :
213 0 : if (lfill_halo) then
214 : !
215 : !*********************************************
216 : !
217 : ! halo exchange
218 : !
219 : !*********************************************
220 : !
221 0 : call t_startf('fill_halo_and_extend_panel initbuffer')
222 0 : call initghostbuffer(hybrid%par,cellghostbuf,elem,numlev*num_flds,nhcc,nphys)
223 0 : call t_stopf('fill_halo_and_extend_panel initbuffer')
224 0 : do ie=nets,nete
225 0 : call ghostpack(cellghostbuf, fld(:,:,:,:,ie),numlev*num_flds,0,ie)
226 : end do
227 0 : call ghost_exchange(hybrid,cellghostbuf,location='fill_halo_and_extend_panel')
228 0 : do ie=nets,nete
229 0 : call ghostunpack(cellghostbuf, fld(:,:,:,:,ie),numlev*num_flds,0,ie)
230 : end do
231 0 : call freeghostbuffer(cellghostbuf)
232 : end if
233 0 : if (lextend_panel) then
234 : !
235 : !*********************************************
236 : !
237 : ! extend panel
238 : !
239 : !*********************************************
240 : !
241 0 : if (nphys==fv_nphys) then
242 0 : if (ndepth>nhr_phys) &
243 0 : call endrun("fill_halo_and_extend_panel: ndepth>nhr_phys")
244 0 : nht_phys = nhe_phys+nhr_phys
245 0 : nh_phys = nhr_phys
246 0 : allocate(fld_tmp(1-nht_phys:nphys+nht_phys,1-nht_phys:nphys+nht_phys))
247 0 : do ie=nets,nete
248 0 : do itr=1,num_flds
249 0 : do k=1,numlev
250 : call extend_panel_interpolate(fv_nphys,nhc_phys,nhr_phys,nht_phys,ns_phys,nh_phys,&
251 0 : fld(:,:,k,itr,ie),fvm(ie)%cubeboundary,&
252 0 : fvm(ie)%halo_interp_weight_physgrid(1:ns_phys,1-nh_phys:fv_nphys+nh_phys,1:nhr_phys,:),&
253 0 : fvm(ie)%ibase_physgrid(1-nh_phys:fv_nphys+nh_phys,1:nhr_phys,:),&
254 0 : fld_tmp)
255 0 : fld(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth,k,itr,ie) = fld_tmp(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth)
256 : end do
257 : end do
258 : end do
259 0 : deallocate(fld_tmp)
260 0 : else if (nphys==nc) then
261 0 : if (ndepth>nhr) &
262 0 : call endrun("fill_halo_and_extend_panel: ndepth>nhr")
263 0 : nhe_phys= 0
264 0 : nht_phys= nhe_phys+nhr
265 0 : nh_phys = nhr
266 0 : allocate(fld_tmp(1-nht_phys:nphys+nht_phys,1-nht_phys:nphys+nht_phys))
267 0 : do ie=nets,nete
268 0 : do itr=1,num_flds
269 0 : do k=1,numlev
270 : call extend_panel_interpolate(nc,nhc,nhr,nht_phys,ns,nh_phys,&
271 0 : fld(:,:,k,itr,ie),fvm(ie)%cubeboundary,&
272 0 : fvm(ie)%halo_interp_weight(1:ns,1-nh_phys:nc+nh_phys,1:nhr,:),&
273 : fvm(ie)%ibase(1-nh_phys:nc+nh_phys,1:nhr,:),&
274 0 : fld_tmp)
275 0 : fld(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth,k,itr,ie) = fld_tmp(1-ndepth:nphys+ndepth,1-ndepth:nphys+ndepth)
276 : end do
277 : end do
278 : end do
279 0 : deallocate(fld_tmp)
280 : else
281 0 : call endrun("fill_halo_and_extend_panel: resolution not supported")
282 : end if
283 : end if
284 0 : end subroutine fill_halo_and_extend_panel
285 :
286 :
287 : ! initialize global buffers shared by all threads
288 1536 : subroutine fvm_init1(par,elem)
289 0 : use parallel_mod, only: parallel_t
290 : use cam_abortutils, only: endrun
291 : use cam_logfile, only: iulog
292 : use control_mod, only: rsplit
293 : use dimensions_mod, only: qsize, qsize_d
294 : use dimensions_mod, only: fvm_supercycling, fvm_supercycling_jet
295 : use dimensions_mod, only: nc,nhe, nhc, nlev,ntrac, ntrac_d,ns, nhr, use_cslam
296 : use dimensions_mod, only: large_Courant_incr
297 : use dimensions_mod, only: kmin_jet,kmax_jet
298 :
299 : type (parallel_t) :: par
300 : type (element_t),intent(inout) :: elem(:)
301 : !
302 1536 : if (use_cslam) then
303 1536 : if (par%masterproc) then
304 2 : write(iulog,*) " "
305 2 : write(iulog,*) "|-----------------------------------------|"
306 2 : write(iulog,*) "| FVM tracer transport scheme information |"
307 2 : write(iulog,*) "|-----------------------------------------|"
308 2 : write(iulog,*) " "
309 : end if
310 1536 : if (use_cslam) then
311 1536 : if (par%masterproc) then
312 2 : write(iulog,*) "Running consistent SE-CSLAM, Lauritzen et al. (2017, MWR)."
313 2 : write(iulog,*) "CSLAM = Conservative Semi-LAgrangian Multi-tracer scheme"
314 2 : write(iulog,*) "Lauritzen et al., (2010), J. Comput. Phys."
315 2 : write(iulog,*) " "
316 : end if
317 : end if
318 : !
319 : ! PARAMETER ERROR CHECKING
320 : !
321 1536 : if (kmin_jet>kmax_jet) &
322 0 : call endrun("PARAMETER ERROR for fvm: kmin_jet must be < kmax_jet")
323 1536 : if (ntrac>ntrac_d) &
324 0 : call endrun("PARAMETER ERROR for fvm: ntrac > ntrac_d")
325 :
326 1536 : if (qsize>0.and.mod(rsplit,fvm_supercycling).ne.0) then
327 0 : if (par%masterproc) then
328 0 : write(iulog,*)'cannot supercycle fvm tracers with respect to se tracers'
329 0 : write(iulog,*)'with this choice of rsplit =',rsplit
330 0 : write(iulog,*)'rsplit must be a multiple of fvm_supercycling=',fvm_supercycling
331 : end if
332 0 : call endrun("PARAMETER ERROR for fvm: mod(rsplit,fvm_supercycling)<>0")
333 : endif
334 :
335 1536 : if (qsize>0.and.mod(rsplit,fvm_supercycling_jet).ne.0) then
336 0 : if (par%masterproc) then
337 0 : write(iulog,*)'cannot supercycle fvm tracers with respect to se tracers'
338 0 : write(iulog,*)'with this choice of rsplit =',rsplit
339 0 : write(iulog,*)'rsplit must be a multiple of fvm_supercycling_jet=',fvm_supercycling_jet
340 : end if
341 0 : call endrun("PARAMETER ERROR for fvm: mod(rsplit,fvm_supercycling_jet)<>0")
342 : endif
343 :
344 1536 : if (large_Courant_incr.and.(fvm_supercycling.ne.fvm_supercycling_jet)) then
345 0 : if (par%masterproc) then
346 0 : write(iulog,*)'Large Courant number increment requires no level dependent supercycling'
347 0 : write(iulog,*)'i.e. fvm_supercycling must be equal to fvm_supercycling_jet'
348 : end if
349 0 : call endrun("PARAMETER ERROR for fvm: large_courant_incr requires fvm_supercycling=fvm_supercycling_jet")
350 : endif
351 :
352 1536 : if (par%masterproc) then
353 2 : write(iulog,*) " "
354 2 : write(iulog,*) "Done Tracer transport scheme information "
355 2 : write(iulog,*) " "
356 : end if
357 :
358 :
359 1536 : if (par%masterproc) write(iulog,*) "fvm resolution is nc*nc in each element: nc = ",nc
360 1536 : if (par%masterproc) write(iulog,*)'ntrac,ntrac_d=',ntrac,ntrac_d
361 1536 : if (par%masterproc) write(iulog,*)'qsize,qsize_d=',qsize,qsize_d
362 :
363 : if (nc.ne.3) then
364 : if (par%masterproc) then
365 : write(iulog,*) "Only nc==3 is supported for CSLAM"
366 : endif
367 : call endrun("PARAMETER ERRROR for fvm: only nc=3 supported for CSLAM")
368 : end if
369 :
370 1536 : if (par%masterproc) then
371 2 : write(iulog,*) " "
372 : if (ns==1) then
373 : write(iulog,*) "ns==1: using no interpolation for mapping cell averages values across edges"
374 : write(iulog,*) "Note: this is not a recommended setting - large errors at panel edges!"
375 : else if (ns==2) then
376 : write(iulog,*) "ns==2: using linear interpolation for mapping cell averages values across edges"
377 : write(iulog,*) "Note that ns=4 is default CSLAM setting used in Lauritzen et al. (2010)"
378 : write(iulog,*) "so this option is slightly less accurate (but the stencil is smaller near panel edges!)"
379 :
380 : else if (ns==3) then
381 2 : write(iulog,*) "ns==3: using quadratic interpolation for mapping cell averages values across edges"
382 2 : write(iulog,*) "Note that ns=4 is default CSLAM setting used in Lauritzen et al. (2010)"
383 2 : write(iulog,*) "so this option is slightly less accurate (but the stencil is smaller near panel edges!)"
384 : else if (ns==4) then
385 : write(iulog,*) "ns==4: using cubic interpolation for mapping cell averages values across edges"
386 : write(iulog,*) "This is default CSLAM setting used in Lauritzen et al. (2010)"
387 : else
388 : write(iulog,*) "Not a tested value for ns but it should work! You choose ns = ",ns
389 : end if
390 :
391 : ! if (ns.NE.3) then
392 : ! write(*,*) "In fvm_reconstruction_mod function matmul_w has been hard-coded for ns=3 for performance"
393 : ! write(*,*) "Revert to general code - outcommented above"
394 : ! call endrun("stopping")
395 : ! end if
396 : end if
397 :
398 : if (MOD(ns,2)==0.and.nhr+(nhe-1)+ns/2>nc+nc) then
399 : write(iulog,*) "to run this combination of ns and nhr you need to increase nc to ",nhr+ns/2+nhe-1
400 : write(iulog,*) "You choose (ns,nhr,nc,nhe)=",ns,nhr,nc,nhe
401 : call endrun("stopping")
402 : end if
403 : if (MOD(ns,2)==1.and.nhr+(ns-1)/2+(nhe-1)>nc+nc) then
404 : write(iulog,*) "to run this combination of ns and nhr you need to increase nc to ",nhr+(ns-1)/2+nhe-1
405 : write(iulog,*) "You choose (ns,nhr,nc,nhe)=",ns,nhr,nc,nhe
406 : call endrun("stopping")
407 : end if
408 :
409 : if (nc==3.and.ns.ne.3) then
410 : if (par%masterproc) then
411 : write(iulog,*) "Recommended setting for nc=3 is ns=3 (linear interpolation in halo)"
412 : write(iulog,*) "You choose ns=",ns
413 : write(iulog,*) "Goto dimensions_mod to change value of ns"
414 : write(iulog,*) "or outcomment call haltmop below (i.e. you know what you are doing!)"
415 : endif
416 : call endrun("stopping")
417 : end if
418 : end if
419 :
420 : if (nc==4.and.ns.ne.4) then
421 : if (par%masterproc) then
422 : write(iulog,*) "Recommended setting for nc=4 is ns=4 (cubic interpolation in halo)"
423 : write(iulog,*) "You choose ns=",ns
424 : write(iulog,*) "Goto dimensions_mod to change value of ns"
425 : write(iulog,*) "or outcomment call haltmop below (i.e. you know what you are doing!)"
426 : endif
427 : call endrun("stopping")
428 : end if
429 :
430 : if (nhe .ne. 1) then
431 : if (par%masterproc) then
432 : write(iulog,*) "PARAMETER ERROR for fvm: Number of halo zone for the extended"
433 : write(iulog,*) "element nhe has to be 1, only this is available now! STOP!"
434 : endif
435 : call endrun("stopping")
436 : end if
437 1536 : end subroutine fvm_init1
438 :
439 :
440 :
441 :
442 :
443 : ! initialization that can be done in threaded regions
444 1536 : subroutine fvm_init2(elem,fvm,hybrid,nets,nete)
445 : use fvm_control_volume_mod, only: fvm_mesh,fvm_set_cubeboundary
446 : use bndry_mod, only: compute_ghost_corner_orientation
447 : use dimensions_mod, only: nlev, nc, nhc, nhe, ntrac, ntrac_d, np
448 : use dimensions_mod, only: nhc_phys, fv_nphys
449 : use dimensions_mod, only: fvm_supercycling, fvm_supercycling_jet
450 : use dimensions_mod, only: kmin_jet,kmax_jet
451 : use hycoef, only: hyai, hybi, ps0
452 : use derivative_mod, only: subcell_integration
453 : use air_composition, only: thermodynamic_active_species_num
454 :
455 : type (fvm_struct) :: fvm(:)
456 : type (element_t) :: elem(:)
457 : type (hybrid_t) :: hybrid
458 : integer :: ie,nets,nete,k,klev
459 : real(kind=r8) :: one(np,np)
460 :
461 32256 : one = 1.0_r8
462 12336 : do ie=nets,nete
463 293136 : do k = 1, nlev
464 280800 : fvm(ie)%dp_ref(k) = ( hyai(k+1) - hyai(k) )*ps0 + ( hybi(k+1) - hybi(k) )*ps0
465 291600 : fvm(ie)%dp_ref_inverse(k) = 1.0_r8/fvm(ie)%dp_ref(k)
466 : end do
467 : end do
468 :
469 1536 : call compute_ghost_corner_orientation(hybrid,elem,nets,nete)
470 : ! run some tests:
471 : ! call test_ghost(hybrid,elem,nets,nete)
472 :
473 12336 : do ie=nets,nete
474 10800 : call fvm_set_cubeboundary(elem(ie),fvm(ie))
475 : call fvm_mesh(elem(ie),fvm(ie))
476 140400 : fvm(ie)%inv_area_sphere = 1.0_r8/fvm(ie)%area_sphere
477 : !
478 : ! compute CSLAM areas consistent with SE area (at 1 degree they can be up to
479 : ! 1E-6 different than the correct spherical areas used in CSLAM)
480 : !
481 10800 : call subcell_integration(one, np, nc, elem(ie)%metdet,fvm(ie)%inv_se_area_sphere)
482 140400 : fvm(ie)%inv_se_area_sphere = 1.0_r8/fvm(ie)%inv_se_area_sphere
483 :
484 36622800 : fvm(ie)%fc(:,:,:,:) = 0.0_r8
485 51397200 : fvm(ie)%fm(:,:,:,:) = 0.0_r8
486 25575936 : fvm(ie)%ft(:,:,: ) = 0.0_r8
487 : enddo
488 : ! Need to allocate ghostBufQnhc after compute_ghost_corner_orientation because it
489 : ! changes the values for reverse
490 :
491 1536 : call initghostbuffer(hybrid%par,ghostBufQnhc_s,elem,nlev*(ntrac+1),nhc,nc,nthreads=1)
492 1536 : call initghostbuffer(hybrid%par,ghostBufQnhc_t1,elem,nlev, nhc,nc,nthreads=1)
493 1536 : call initghostbuffer(hybrid%par,ghostBufQnhc_h,elem,nlev*(ntrac+1),nhc,nc,nthreads=horz_num_threads)
494 1536 : call initghostbuffer(hybrid%par,ghostBufQnhc_vh,elem,nlev*(ntrac+1),nhc,nc,nthreads=vert_num_threads*horz_num_threads)
495 1536 : klev = kmax_jet-kmin_jet+1
496 1536 : call initghostbuffer(hybrid%par,ghostBufQ1_h,elem,klev*(ntrac+1),1,nc,nthreads=horz_num_threads)
497 1536 : call initghostbuffer(hybrid%par,ghostBufQ1_vh,elem,klev*(ntrac+1),1,nc,nthreads=vert_num_threads*horz_num_threads)
498 : ! call initghostbuffer(hybrid%par,ghostBufFlux_h,elem,4*nlev,nhe,nc,nthreads=horz_num_threads)
499 1536 : call initghostbuffer(hybrid%par,ghostBufFlux_vh,elem,4*nlev,nhe,nc,nthreads=vert_num_threads*horz_num_threads)
500 1536 : call initghostbuffer(hybrid%par,ghostBuf_cslam2gll,elem,nlev*thermodynamic_active_species_num,nhc,nc,nthreads=1)
501 : !
502 : ! preallocate buffers for physics-dynamics coupling
503 : !
504 1536 : if (fv_nphys.ne.nc) then
505 0 : call initghostbuffer(hybrid%par,ghostBufPG_s,elem,nlev*(4+ntrac),nhc_phys,fv_nphys,nthreads=1)
506 : else
507 1536 : call initghostbuffer(hybrid%par,ghostBufPG_s,elem,nlev*3,nhc_phys,fv_nphys,nthreads=1)
508 : end if
509 :
510 1536 : if (fvm_supercycling.ne.fvm_supercycling_jet) then
511 : !
512 : ! buffers for running different fvm time-steps in the jet region
513 : !
514 0 : klev = kmax_jet-kmin_jet+1
515 0 : call initghostbuffer(hybrid%par,ghostBufQnhcJet_h,elem,klev*(ntrac+1),nhc,nc,nthreads=horz_num_threads)
516 0 : call initghostbuffer(hybrid%par,ghostBufFluxJet_h,elem,4*klev,nhe,nc,nthreads=horz_num_threads)
517 : end if
518 1536 : end subroutine fvm_init2
519 :
520 :
521 1536 : subroutine fvm_init3(elem,fvm,hybrid,nets,nete,irecons)
522 1536 : use control_mod , only: neast, nwest, seast, swest
523 : use fvm_analytic_mod, only: compute_reconstruct_matrix
524 : use dimensions_mod , only: fv_nphys, use_cslam
525 : use dimensions_mod, only: nlev, nc, nhe, nlev, nhc
526 : use coordinate_systems_mod, only: cartesian2D_t,cartesian3D_t
527 : use coordinate_systems_mod, only: cubedsphere2cart, cart2cubedsphere
528 : implicit none
529 : type (element_t) ,intent(inout) :: elem(:)
530 : type (fvm_struct),intent(inout) :: fvm(:)
531 : type (hybrid_t) ,intent(in) :: hybrid
532 : integer ,intent(in) :: nets,nete,irecons
533 : !
534 1536 : type (edgeBuffer_t) :: cellghostbuf
535 : integer :: ie, ixy, ivertex, i, j,istart,itot,ishft,imin,imax
536 : integer, dimension(2,4) :: unit_vec
537 : integer :: rot90_matrix(2,2), iside
538 :
539 : type (cartesian2D_t) :: tmpgnom
540 : type (cartesian2D_t) :: gnom
541 : type(cartesian3D_t) :: tmpcart3d
542 :
543 1536 : if (use_cslam.and.nc.ne.fv_nphys) then
544 : !
545 : ! fill the fvm halo for mapping in d_p_coupling if
546 : ! physics grid resolution is different than fvm resolution
547 : !
548 0 : call fill_halo_fvm(elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
549 : end if
550 :
551 :
552 1536 : imin=1-nhc
553 1536 : imax=nc+nhc
554 : !
555 : ! fill halo start
556 : !
557 1536 : itot=9+irecons-1+2
558 1536 : call initghostbuffer(hybrid%par,cellghostbuf,elem,itot,nhc,nc)
559 12336 : do ie=nets,nete
560 10800 : istart = 0
561 982800 : call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord(1,:,:),1,istart,ie)
562 10800 : istart = istart+1
563 982800 : call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord(2,:,:),1,istart,ie)
564 10800 : istart = istart+1
565 32400 : do ixy=1,2
566 118800 : do ivertex=1,4
567 7862400 : call ghostpack(cellghostbuf, fvm(ie)%vtx_cart(ivertex,ixy,:,:) ,1,istart,ie)
568 108000 : istart = istart+1
569 : end do
570 : end do
571 982800 : call ghostpack(cellghostbuf, fvm(ie)%flux_orient(1,:,:) ,1,istart,ie)
572 66336 : do ixy=1,irecons-1
573 54000 : istart=istart+1
574 4924800 : call ghostpack(cellghostbuf, fvm(ie)%spherecentroid(ixy,:,:) ,1,istart,ie)
575 : end do
576 : end do
577 1536 : call ghost_exchange(hybrid,cellghostbuf,location='fvm_init3')
578 12336 : do ie=nets,nete
579 10800 : istart = 0
580 1954800 : call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord(1,:,:),1,istart,ie)
581 10800 : istart = istart+1
582 1954800 : call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord(2,:,:),1,istart,ie)
583 10800 : istart = istart+1
584 32400 : do ixy=1,2
585 118800 : do ivertex=1,4
586 15638400 : call ghostunpack(cellghostbuf, fvm(ie)%vtx_cart(ivertex,ixy,:,:) ,1,istart,ie)
587 108000 : istart = istart+1
588 : end do
589 : end do
590 1954800 : call ghostunpack(cellghostbuf, fvm(ie)%flux_orient(1,:,:) ,1,istart,ie)
591 66336 : do ixy=1,irecons-1
592 54000 : istart=istart+1
593 9784800 : call ghostunpack(cellghostbuf, fvm(ie)%spherecentroid(ixy,:,:) ,1,istart,ie)
594 : end do
595 : enddo
596 1536 : call freeghostbuffer(cellghostbuf)
597 : !
598 : ! indicator for non-existing cells
599 : ! set vtx_cart to corner value in non-existent cells
600 : !
601 12336 : do ie=nets,nete
602 12336 : if (fvm(ie)%cubeboundary==nwest) then
603 372 : fvm(ie)%flux_orient (: ,1-nhc :0 ,nc +1 :nc +nhc ) = -1
604 696 : fvm(ie)%spherecentroid (:, 1-nhc :0 ,nc +1 :nc +nhc ) = -1e5_r8
605 588 : fvm(ie)%vtx_cart(:,1,1-nhc:0 ,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(4,1,1,nc)
606 588 : fvm(ie)%vtx_cart(:,2,1-nhc:0 ,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(4,2,1,nc)
607 10788 : else if (fvm(ie)%cubeboundary==swest) then
608 372 : fvm(ie)%flux_orient (:,1-nhc :0 ,1-nhc :0 ) = -1
609 696 : fvm(ie)%spherecentroid (:,1-nhc :0 ,1-nhc :0 ) = -1e5_r8
610 588 : fvm(ie)%vtx_cart(:,1,1-nhc:0 ,1-nhc:0 ) = fvm(ie)%vtx_cart(1,1,1,1)
611 588 : fvm(ie)%vtx_cart(:,2,1-nhc:0 ,1-nhc:0 ) = fvm(ie)%vtx_cart(1,2,1,1)
612 10776 : else if (fvm(ie)%cubeboundary==neast) then
613 372 : fvm(ie)%flux_orient (:,nc +1 :nc +nhc ,nc +1 :nc +nhc ) = -1
614 696 : fvm(ie)%spherecentroid (:,nc +1 :nc +nhc ,nc +1 :nc +nhc ) = -1e5_r8
615 588 : fvm(ie)%vtx_cart(:,1,nc+1 :nc+nhc,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(3,1,nc,nc)
616 588 : fvm(ie)%vtx_cart(:,2,nc+1 :nc+nhc,nc+1 :nc+nhc) = fvm(ie)%vtx_cart(3,2,nc,nc)
617 10764 : else if (fvm(ie)%cubeboundary==seast) then
618 372 : fvm(ie)%flux_orient (:,nc +1 :nc +nhc ,1-nhc :0 ) = -1
619 696 : fvm(ie)%spherecentroid (:,nc +1 :nc +nhc ,1-nhc :0 ) = -1e5_r8
620 588 : fvm(ie)%vtx_cart(:,1,nc+1 :nc+nhc,1-nhc:0 ) = fvm(ie)%vtx_cart(2,1,nc,1)
621 588 : fvm(ie)%vtx_cart(:,2,nc+1 :nc+nhc,1-nhc:0 ) = fvm(ie)%vtx_cart(2,2,nc,1)
622 : end if
623 : end do
624 :
625 : !
626 : ! set vectors for perpendicular flux vector
627 : !
628 1536 : rot90_matrix(1,1) = 0; rot90_matrix(2,1) = 1 !counter-clockwise rotation matrix
629 1536 : rot90_matrix(1,2) =-1; rot90_matrix(2,2) = 0 !counter-clockwise rotation matrix
630 :
631 1536 : iside = 1
632 1536 : unit_vec(1,iside) = 0 !x-component of displacement vector for side 1
633 1536 : unit_vec(2,iside) = 1 !y-component of displacement vector for side 1
634 :
635 6144 : do iside=2,4
636 6144 : unit_vec(:,iside) = MATMUL(rot90_matrix(:,:),unit_vec(:,iside-1))
637 : end do
638 :
639 : !
640 : ! fill halo done
641 : !
642 : !-------------------------------
643 :
644 12336 : do ie=nets,nete
645 3942000 : fvm(ie)%displ_max = 0.0_r8
646 109536 : do j=imin,imax
647 982800 : do i=imin,imax
648 : !
649 : ! rotate gnomonic coordinate vector
650 : !
651 : ! fvm(ie)%norm_elem_coord(:,i,j) = MATMUL(fvm(ie)%rot_matrix(:,:,i,j),fvm(ie)%norm_elem_coord(:,i,j))
652 : !
653 874800 : ishft = NINT(fvm(ie)%flux_orient(2,i,j))
654 2721600 : do ixy=1,2
655 : !
656 : ! rotate coordinates if needed through permutation
657 : !
658 8748000 : fvm(ie)%vtx_cart(1:4,ixy,i,j) = cshift(fvm(ie)%vtx_cart(1:4,ixy,i,j),shift=ishft)
659 8748000 : fvm(ie)%flux_vec (ixy,i,j,1:4) = cshift(unit_vec (ixy,1:4 ),shift=ishft)
660 : !
661 : ! set flux vector to zero in non-existent cells (corner halo)
662 : !
663 8748000 : fvm(ie)%flux_vec (ixy,i,j,1:4) = fvm(ie)%ifct(i,j)*fvm(ie)%flux_vec(ixy,i,j,1:4)
664 :
665 1749600 : iside=1
666 : fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
667 1749600 : ABS(fvm(ie)%vtx_cart(4,ixy,i,j)-fvm(ie)%vtx_cart(1,ixy,i,j))
668 1749600 : iside=2
669 : fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
670 1749600 : ABS(fvm(ie)%vtx_cart(1,ixy,i,j)-fvm(ie)%vtx_cart(2,ixy,i,j))
671 1749600 : iside=3
672 : fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
673 1749600 : ABS(fvm(ie)%vtx_cart(2,ixy,i,j)-fvm(ie)%vtx_cart(3,ixy,i,j))
674 1749600 : iside=4
675 : fvm(ie)%displ_max(i,j,iside) = fvm(ie)%displ_max(i,j,iside)+&
676 2624400 : ABS(fvm(ie)%vtx_cart(2,ixy,i,j)-fvm(ie)%vtx_cart(1,ixy,i,j))
677 : end do
678 : end do
679 : end do
680 : end do
681 : !
682 : ! pre-compute derived metric terms used for integration, polynomial
683 : ! evaluation at fvm cell vertices, etc.
684 : !
685 12336 : do ie=nets,nete
686 10800 : call compute_reconstruct_matrix(nc,nhe,nhc,irecons,fvm(ie)%dalpha,fvm(ie)%dbeta,&
687 : fvm(ie)%spherecentroid,fvm(ie)%vtx_cart,fvm(ie)%centroid_stretch,&
688 23136 : fvm(ie)%vertex_recons_weights,fvm(ie)%recons_metrics,fvm(ie)%recons_metrics_integral)
689 : end do
690 : !
691 : ! create a normalized element coordinate system with a halo
692 : !
693 12336 : do ie=nets,nete
694 109536 : do j=1-nhc,nc+nhc
695 982800 : do i=1-nhc,nc+nhc
696 : !
697 : ! only compute for physically existent cells
698 : !
699 972000 : if (fvm(ie)%ifct(i,j)>0) then
700 874368 : gnom%x = fvm(ie)%norm_elem_coord(1,i,j)
701 874368 : gnom%y = fvm(ie)%norm_elem_coord(2,i,j)
702 : !
703 : ! coordinate transform only necessary for points on another panel
704 : !
705 874368 : if (NINT(fvm(ie)%flux_orient(1,1,1)).NE.NINT(fvm(ie)%flux_orient(1,i,j))) then
706 38016 : tmpcart3d=cubedsphere2cart(gnom,NINT(fvm(ie)%flux_orient(1,i,j)))
707 38016 : tmpgnom=cart2cubedsphere(tmpcart3d,NINT(fvm(ie)%flux_orient(1,1,1)))
708 : else
709 836352 : tmpgnom%x = fvm(ie)%norm_elem_coord(1,i,j)
710 836352 : tmpgnom%y = fvm(ie)%norm_elem_coord(2,i,j)
711 : end if
712 : !
713 : ! convert to element normalized coordinates
714 : !
715 874368 : fvm(ie)%norm_elem_coord(1,i,j) =(tmpgnom%x-elem(ie)%corners(1)%x)/&
716 874368 : (0.5_r8*dble(nc)*fvm(ie)%dalpha)-1.0_r8
717 : fvm(ie)%norm_elem_coord(2,i,j) =(tmpgnom%y-elem(ie)%corners(1)%y)/&
718 874368 : (0.5_r8*dble(nc)*fvm(ie)%dalpha)-1.0_r8
719 : else
720 432 : fvm(ie)%norm_elem_coord(1,i,j) = 1D9
721 432 : fvm(ie)%norm_elem_coord(2,i,j) = 1D9
722 : end if
723 : end do
724 : end do
725 : end do
726 :
727 1536 : end subroutine fvm_init3
728 :
729 :
730 1536 : subroutine fvm_pg_init(elem, fvm, hybrid, nets, nete,irecons)
731 : use coordinate_systems_mod, only : cartesian2D_t,cartesian3D_t
732 : use control_mod, only : neast, nwest, seast, swest
733 : use coordinate_systems_mod, only : cubedsphere2cart, cart2cubedsphere
734 : use dimensions_mod, only: fv_nphys, nhe_phys,nhc_phys
735 : use cube_mod ,only: dmap
736 : use control_mod ,only: cubed_sphere_map
737 : use fvm_analytic_mod, only: compute_reconstruct_matrix
738 :
739 : type (element_t) , intent(in) :: elem(:)
740 : type (fvm_struct), intent(inout) :: fvm(:)
741 : type (hybrid_t) , intent(in) :: hybrid
742 :
743 : type (cartesian2D_t) :: gnom
744 : type(cartesian3D_t) :: tmpcart3d
745 : type (cartesian2D_t) :: tmpgnom
746 :
747 :
748 : integer, intent(in) :: nets ! starting thread element number (private)
749 : integer, intent(in) :: nete,irecons ! ending thread element number (private)
750 :
751 : ! ==================================
752 : ! Local variables
753 : ! ==================================
754 :
755 : integer :: ie, ixy, ivertex, i, j,istart,itot,ishft,imin,imax
756 : integer, dimension(2,4) :: unit_vec
757 : integer :: rot90_matrix(2,2), iside
758 :
759 1536 : type (edgeBuffer_t) :: cellghostbuf
760 :
761 : ! D is derivative of gnomonic mapping
762 3072 : real (kind=r8) :: D(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,2,2)
763 : real (kind=r8) :: detD,x1,x2
764 :
765 1536 : if (fv_nphys>0) then
766 : !
767 : ! do the same as fvm_init3 for the metric terms of physgrid
768 : !
769 1536 : imin=1-nhc_phys
770 1536 : imax=fv_nphys+nhc_phys
771 : !
772 : ! fill halo start
773 : !
774 1536 : itot=9+irecons-1+2
775 1536 : call initghostbuffer(hybrid%par,cellghostbuf,elem,itot,nhc_phys,fv_nphys)
776 12336 : do ie=nets,nete
777 10800 : istart = 0
778 982800 : call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(1,:,:),1,istart,ie)
779 10800 : istart = istart+1
780 982800 : call ghostpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(2,:,:),1,istart,ie)
781 10800 : istart = istart+1
782 32400 : do ixy=1,2
783 118800 : do ivertex=1,4
784 7862400 : call ghostpack(cellghostbuf, fvm(ie)%vtx_cart_physgrid(ivertex,ixy,:,:) ,1,istart,ie)
785 108000 : istart = istart+1
786 : end do
787 : end do
788 982800 : call ghostpack(cellghostbuf, fvm(ie)%flux_orient_physgrid(1,:,:) ,1,istart,ie)
789 66336 : do ixy=1,irecons-1
790 54000 : istart=istart+1
791 4924800 : call ghostpack(cellghostbuf, fvm(ie)%spherecentroid_physgrid(ixy,:,:) ,1,istart,ie)
792 : end do
793 : end do
794 1536 : call ghost_exchange(hybrid,cellghostbuf,location='fvm_pg_init')
795 12336 : do ie=nets,nete
796 10800 : istart = 0
797 1954800 : call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(1,:,:),1,istart,ie)
798 10800 : istart = istart+1
799 1954800 : call ghostunpack(cellghostbuf, fvm(ie)%norm_elem_coord_physgrid(2,:,:),1,istart,ie)
800 10800 : istart = istart+1
801 32400 : do ixy=1,2
802 118800 : do ivertex=1,4
803 15638400 : call ghostunpack(cellghostbuf, fvm(ie)%vtx_cart_physgrid(ivertex,ixy,:,:) ,1,istart,ie)
804 108000 : istart = istart+1
805 : end do
806 : end do
807 1954800 : call ghostunpack(cellghostbuf, fvm(ie)%flux_orient_physgrid(1,:,:) ,1,istart,ie)
808 66336 : do ixy=1,irecons-1
809 54000 : istart=istart+1
810 9784800 : call ghostunpack(cellghostbuf, fvm(ie)%spherecentroid_physgrid(ixy,:,:) ,1,istart,ie)
811 : end do
812 : enddo
813 1536 : call freeghostbuffer(cellghostbuf)
814 : !
815 : ! indicator for non-existing cells
816 : ! set vtx_cart to corner value in non-existent cells
817 : !
818 12336 : do ie=nets,nete
819 12336 : if (fvm(ie)%cubeboundary==nwest) then
820 372 : fvm(ie)%flux_orient_physgrid (: ,1-nhc_phys :0 ,fv_nphys +1 :fv_nphys +nhc_phys ) = -1
821 696 : fvm(ie)%spherecentroid_physgrid(:, 1-nhc_phys :0 ,fv_nphys +1 :fv_nphys +nhc_phys ) = -1e5_r8
822 0 : fvm(ie)%vtx_cart_physgrid(:,1,1-nhc_phys:0 ,fv_nphys+1 :fv_nphys+nhc_phys) = &
823 588 : fvm(ie)%vtx_cart_physgrid(4,1,1,fv_nphys)
824 0 : fvm(ie)%vtx_cart_physgrid(:,2,1-nhc_phys:0 ,fv_nphys+1 :fv_nphys+nhc_phys) = &
825 588 : fvm(ie)%vtx_cart_physgrid(4,2,1,fv_nphys)
826 10788 : else if (fvm(ie)%cubeboundary==swest) then
827 372 : fvm(ie)%flux_orient_physgrid (:,1-nhc_phys :0 ,1-nhc_phys :0 ) = -1
828 696 : fvm(ie)%spherecentroid_physgrid(:,1-nhc_phys :0 ,1-nhc_phys :0 ) = -1e5_r8
829 588 : fvm(ie)%vtx_cart_physgrid(:,1,1-nhc_phys:0 ,1-nhc_phys:0 ) = fvm(ie)%vtx_cart_physgrid(1,1,1,1)
830 588 : fvm(ie)%vtx_cart_physgrid(:,2,1-nhc_phys:0 ,1-nhc_phys:0 ) = fvm(ie)%vtx_cart_physgrid(1,2,1,1)
831 10776 : else if (fvm(ie)%cubeboundary==neast) then
832 0 : fvm(ie)%flux_orient_physgrid (:,fv_nphys +1 :fv_nphys +nhc_phys , &
833 372 : fv_nphys +1 :fv_nphys +nhc_phys ) = -1
834 0 : fvm(ie)%spherecentroid_physgrid(:,fv_nphys +1 :fv_nphys +nhc_phys , &
835 696 : fv_nphys +1 :fv_nphys +nhc_phys ) = -1e5_r8
836 0 : fvm(ie)%vtx_cart_physgrid(:,1,fv_nphys+1 :fv_nphys+nhc_phys,fv_nphys+1 :fv_nphys+nhc_phys) = &
837 588 : fvm(ie)%vtx_cart_physgrid(3,1,fv_nphys,fv_nphys)
838 0 : fvm(ie)%vtx_cart_physgrid(:,2,fv_nphys+1 :fv_nphys+nhc_phys,fv_nphys+1 :fv_nphys+nhc_phys) = &
839 588 : fvm(ie)%vtx_cart_physgrid(3,2,fv_nphys,fv_nphys)
840 10764 : else if (fvm(ie)%cubeboundary==seast) then
841 372 : fvm(ie)%flux_orient_physgrid (:,fv_nphys +1 :fv_nphys +nhc_phys ,1-nhc_phys :0 ) = -1
842 696 : fvm(ie)%spherecentroid_physgrid(:,fv_nphys +1 :fv_nphys +nhc_phys ,1-nhc_phys :0 ) = -1e5_r8
843 0 : fvm(ie)%vtx_cart_physgrid(:,1,fv_nphys+1 :fv_nphys+nhc_phys,1-nhc_phys:0 ) = &
844 588 : fvm(ie)%vtx_cart_physgrid(2,1,fv_nphys,1)
845 0 : fvm(ie)%vtx_cart_physgrid(:,2,fv_nphys+1 :fv_nphys+nhc_phys,1-nhc_phys:0 ) = &
846 588 : fvm(ie)%vtx_cart_physgrid(2,2,fv_nphys,1)
847 : end if
848 : end do
849 :
850 : !
851 : ! set vectors for perpendicular flux vector
852 : !
853 1536 : rot90_matrix(1,1) = 0; rot90_matrix(2,1) = 1 !counter-clockwise rotation matrix
854 1536 : rot90_matrix(1,2) =-1; rot90_matrix(2,2) = 0 !counter-clockwise rotation matrix
855 :
856 1536 : iside = 1
857 1536 : unit_vec(1,iside) = 0 !x-component of displacement vector for side 1
858 1536 : unit_vec(2,iside) = 1 !y-component of displacement vector for side 1
859 :
860 6144 : do iside=2,4
861 6144 : unit_vec(:,iside) = MATMUL(rot90_matrix(:,:),unit_vec(:,iside-1))
862 : end do
863 :
864 : !
865 : ! fill halo done
866 : !
867 : !-------------------------------
868 :
869 12336 : do ie=nets,nete
870 109536 : do j=imin,imax
871 982800 : do i=imin,imax
872 : !
873 : ! rotate gnomonic coordinate vector
874 : !
875 874800 : ishft = NINT(fvm(ie)%flux_orient_physgrid(2,i,j))
876 2721600 : do ixy=1,2
877 : !
878 : ! rotate coordinates if needed through permutation
879 : !
880 9622800 : fvm(ie)%vtx_cart_physgrid(1:4,ixy,i,j) = cshift(fvm(ie)%vtx_cart_physgrid(1:4,ixy,i,j),shift=ishft)
881 : end do
882 : end do
883 : end do
884 : end do
885 : !
886 : ! pre-compute derived metric terms used for integration, polynomial
887 : ! evaluation at fvm cell vertices, etc.
888 : !
889 12336 : do ie=nets,nete
890 10800 : call compute_reconstruct_matrix(fv_nphys,nhe_phys,nhc_phys,irecons,fvm(ie)%dalpha_physgrid,fvm(ie)%dbeta_physgrid,&
891 : fvm(ie)%spherecentroid_physgrid,fvm(ie)%vtx_cart_physgrid,fvm(ie)%centroid_stretch_physgrid,&
892 23136 : fvm(ie)%vertex_recons_weights_physgrid,fvm(ie)%recons_metrics_physgrid,fvm(ie)%recons_metrics_integral_physgrid)
893 : end do
894 : !
895 : ! code specific for physgrid
896 : !
897 : !
898 : ! create a normalized element coordinate system with a halo
899 : !
900 12336 : do ie=nets,nete
901 109536 : do j=1-nhc_phys,fv_nphys+nhc_phys
902 982800 : do i=1-nhc_phys,fv_nphys+nhc_phys
903 : !
904 : ! only compute for physically existent cells
905 : !
906 972000 : if (fvm(ie)%ifct_physgrid(i,j)>0) then
907 874368 : gnom%x = fvm(ie)%norm_elem_coord_physgrid(1,i,j)
908 874368 : gnom%y = fvm(ie)%norm_elem_coord_physgrid(2,i,j)
909 : !
910 : ! coordinate transform only necessary for points on another panel
911 : !
912 874368 : if (NINT(fvm(ie)%flux_orient_physgrid(1,1,1)).NE.NINT(fvm(ie)%flux_orient_physgrid(1,i,j))) then
913 38016 : tmpcart3d=cubedsphere2cart(gnom,NINT(fvm(ie)%flux_orient_physgrid(1,i,j)))
914 38016 : tmpgnom=cart2cubedsphere(tmpcart3d,NINT(fvm(ie)%flux_orient_physgrid(1,1,1)))
915 : else
916 836352 : tmpgnom%x = fvm(ie)%norm_elem_coord_physgrid(1,i,j)
917 836352 : tmpgnom%y = fvm(ie)%norm_elem_coord_physgrid(2,i,j)
918 : end if
919 : !
920 : ! convert to element normalized coordinates
921 : !
922 1748736 : fvm(ie)%norm_elem_coord_physgrid(1,i,j) =(tmpgnom%x-elem(ie)%corners(1)%x)/&
923 1748736 : (0.5_r8*dble(fv_nphys)*fvm(ie)%dalpha_physgrid)-1.0_r8
924 0 : fvm(ie)%norm_elem_coord_physgrid(2,i,j) =(tmpgnom%y-elem(ie)%corners(1)%y)/&
925 874368 : (0.5_r8*dble(fv_nphys)*fvm(ie)%dalpha_physgrid)-1.0_r8
926 : else
927 432 : fvm(ie)%norm_elem_coord_physgrid(1,i,j) = 1D9
928 432 : fvm(ie)%norm_elem_coord_physgrid(2,i,j) = 1D9
929 : end if
930 : end do
931 : end do
932 : end do
933 : !
934 : ! compute Dinv
935 : !
936 12336 : do ie=nets,nete
937 109536 : do j=1-nhc_phys,fv_nphys+nhc_phys
938 982800 : do i=1-nhc_phys,fv_nphys+nhc_phys
939 874800 : x1 = fvm(ie)%norm_elem_coord_physgrid(1,i,j)
940 874800 : x2 = fvm(ie)%norm_elem_coord_physgrid(2,i,j)
941 6123600 : call Dmap(D(i,j,:,:),x1,x2,elem(ie)%corners3D,cubed_sphere_map,elem(ie)%corners,elem(ie)%u2qmap,elem(ie)%facenum)
942 874800 : detD = D(i,j,1,1)*D(i,j,2,2) - D(i,j,1,2)*D(i,j,2,1)
943 :
944 874800 : fvm(ie)%Dinv_physgrid(i,j,1,1) = D(i,j,2,2)/detD
945 874800 : fvm(ie)%Dinv_physgrid(i,j,1,2) = -D(i,j,1,2)/detD
946 874800 : fvm(ie)%Dinv_physgrid(i,j,2,1) = -D(i,j,2,1)/detD
947 972000 : fvm(ie)%Dinv_physgrid(i,j,2,2) = D(i,j,1,1)/detD
948 : end do
949 : end do
950 : end do
951 : end if
952 :
953 1536 : end subroutine fvm_pg_init
954 :
955 :
956 : end module fvm_mod
|
