Line data Source code
1 : module fvm_consistent_se_cslam
2 : use shr_kind_mod, only: r8=>shr_kind_r8
3 : use dimensions_mod, only: nc, nhe, nlev, ntrac, np, nhr, nhc, ngpc, ns, nht
4 : use dimensions_mod, only: irecons_tracer
5 : use dimensions_mod, only: kmin_jet,kmax_jet
6 : use cam_abortutils, only: endrun
7 : use cam_logfile, only: iulog
8 :
9 : use se_dyn_time_mod, only: timelevel_t
10 : use element_mod, only: element_t
11 : use fvm_control_volume_mod, only: fvm_struct
12 : use hybrid_mod, only: hybrid_t, config_thread_region, get_loop_ranges, threadOwnsVertLevel
13 : use perf_mod, only: t_startf, t_stopf
14 : implicit none
15 : private
16 : save
17 :
18 : real (kind=r8),parameter , private :: eps=1.0e-14_r8
19 : public :: run_consistent_se_cslam
20 : contains
21 : !
22 : !**************************************************************************************
23 : !
24 : ! Consistent CSLAM-SE algorithm documented in
25 : !
26 : ! Lauritzen et al. (2017): CAM-SE-CSLAM: Consistent finite-volume transport with
27 : ! spectral-element dynamics. Mon. Wea. Rev.
28 : !
29 : !
30 : !**************************************************************************************
31 : !
32 738816 : subroutine run_consistent_se_cslam(elem,fvm,hybrid,dt_fvm,tl,nets,nete,hvcoord,&
33 : ghostbufQnhc,ghostBufQ1, ghostBufFlux,kminp,kmaxp)
34 : ! ---------------------------------------------------------------------------------
35 : use fvm_mod , only: fill_halo_fvm
36 : use fvm_reconstruction_mod, only: reconstruction
37 : use fvm_analytic_mod , only: gauss_points
38 : use edge_mod , only: ghostpack, ghostunpack
39 : use edgetype_mod , only: edgebuffer_t
40 : use bndry_mod , only: ghost_exchange
41 : use hybvcoord_mod , only: hvcoord_t
42 : use constituents , only: qmin
43 : use dimensions_mod , only: large_Courant_incr,irecons_tracer_lev
44 : use thread_mod , only: vert_num_threads, omp_set_nested
45 : implicit none
46 : type (element_t) , intent(inout) :: elem(:)
47 : type (fvm_struct), target , intent(inout) :: fvm(:)
48 : type (hybrid_t) , intent(in) :: hybrid ! distributed parallel structure (shared)
49 : type (TimeLevel_t) , intent(in) :: tl ! time level struct
50 : type (hvcoord_t) , intent(in) :: hvcoord
51 : integer , intent(in) :: nets ! starting thread element number (private)
52 : integer , intent(in) :: nete ! ending thread element number (private)
53 : real (kind=r8) , intent(in) :: dt_fvm
54 : type (EdgeBuffer_t) , intent(inout) :: ghostbufQnhc,ghostBufQ1, ghostBufFlux
55 : integer , intent(in) :: kminp,kmaxp
56 :
57 : !high-order air density reconstruction
58 1477632 : real (kind=r8) :: ctracer(irecons_tracer,1-nhe:nc+nhe,1-nhe:nc+nhe,ntrac)
59 : real (kind=r8) :: inv_dp_area(nc,nc)
60 : type (hybrid_t) :: hybridnew
61 :
62 : real (kind=r8), dimension(ngpc) :: gsweights, gspts
63 :
64 1477632 : logical :: llimiter(ntrac)
65 : integer :: i,j,k,ie,itr,kptr,q
66 : integer :: kmin_jet_local,kmax_jet_local
67 : integer :: kmin,kmax
68 : integer :: ir
69 : integer :: kblk ! total number of vertical levels per thread
70 : integer :: klev ! total number of vertical levels in the JET region
71 : integer :: region_num_threads
72 : logical :: inJetCall
73 : logical :: ActiveJetThread
74 :
75 738816 : real(r8), pointer :: fcube(:,:,:,:)
76 738816 : real(r8), pointer :: spherecentroid(:,:,:)
77 :
78 2955264 : llimiter = .true.
79 :
80 738816 : inJetCall = .false.
81 738816 : if(((kminp .ne. 1) .or. (kmaxp .ne. nlev)) .and. vert_num_threads>1) then
82 0 : write(iulog,*)'WARNING: deactivating vertical threading for JET region call'
83 0 : inJetCall = .true.
84 0 : region_num_threads = 1
85 : else
86 738816 : region_num_threads = vert_num_threads
87 : endif
88 :
89 738816 : call omp_set_nested(.true.)
90 : !$OMP PARALLEL NUM_THREADS(region_num_threads), DEFAULT(SHARED), &
91 : !$OMP PRIVATE(hybridnew,kblk,ie,k,kmin,gspts,inv_dp_area,itr), &
92 : !$OMP PRIVATE(kmin_jet_local,kmax,kmax_jet_local,kptr,q,ctracer,ActiveJetThread)
93 738816 : call gauss_points(ngpc,gsweights,gspts) !set gauss points/weights
94 2955264 : gspts = 0.5_r8*(gspts+1.0_r8) !shift location so in [0:1] instead of [-1:1]
95 :
96 738816 : if(inJetCall) then
97 : ! ===============================================================================
98 : ! if this is the reduced Jet region call then do not thread over the vertical....
99 : ! Just use the number of vertical levels that were passed into subroutine
100 : ! ===============================================================================
101 0 : hybridnew = config_thread_region(hybrid,'serial')
102 0 : kmin = kminp
103 0 : kmax = kmaxp
104 : else
105 738816 : hybridnew = config_thread_region(hybrid,'vertical')
106 738816 : call get_loop_ranges(hybridnew,kbeg=kmin,kend=kmax)
107 : endif
108 :
109 738816 : kblk = kmax-kmin+1
110 : !call t_startf('fvm:before_Qnhc')
111 5933616 : do ie=nets,nete
112 140259600 : do k=kmin,kmax
113 7158434400 : elem(ie)%sub_elem_mass_flux(:,:,:,k) = dt_fvm*elem(ie)%sub_elem_mass_flux(:,:,:,k)*fvm(ie)%dp_ref_inverse(k)
114 1761037200 : fvm(ie)%dp_fvm(1:nc,1:nc,k) = fvm(ie)%dp_fvm (1:nc,1:nc,k)*fvm(ie)%dp_ref_inverse(k)
115 : end do
116 5194800 : kptr = kmin-1
117 5194800 : call ghostpack(ghostbufQnhc,fvm(ie)%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax) ,kblk, kptr,ie)
118 21518016 : do q=1,ntrac
119 15584400 : kptr = kptr + nlev
120 20779200 : call ghostpack(ghostbufQnhc,fvm(ie)%c(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax,q),kblk,kptr,ie)
121 : enddo
122 : end do
123 : !call t_stopf('fvm:before_Qnhc')
124 : !call t_startf('fvm:ghost_exchange:Qnhc')
125 738816 : call ghost_exchange(hybridnew,ghostbufQnhc,location='ghostbufQnhc')
126 : !call t_stopf('fvm:ghost_exchange:Qnhc')
127 : !call t_startf('fvm:orthogonal_swept_areas')
128 5933616 : do ie=nets,nete
129 140259600 : do k=kmin,kmax
130 7163629200 : fvm(ie)%se_flux (1:nc,1:nc,:,k) = elem(ie)%sub_elem_mass_flux(:,:,:,k)
131 : end do
132 5194800 : kptr = kmin-1
133 5194800 : call ghostunpack(ghostbufQnhc, fvm(ie)%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax) , kblk ,kptr,ie)
134 20779200 : do q=1,ntrac
135 15584400 : kptr = kptr + nlev
136 20779200 : call ghostunpack(ghostbufQnhc, fvm(ie)%c(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax,q),kblk,kptr,ie)
137 : enddo
138 140259600 : do k=kmin,kmax
139 140259600 : call compute_displacements_for_swept_areas (fvm(ie),fvm(ie)%dp_fvm(:,:,k),k,gsweights,gspts)
140 : end do
141 5194800 : kptr = 4*(kmin-1)
142 5933616 : call ghostpack(ghostBufFlux, fvm(ie)%se_flux(:,:,:,kmin:kmax),4*kblk,kptr,ie)
143 : end do
144 :
145 738816 : call ghost_exchange(hybridnew,ghostBufFlux,location='ghostBufFlux')
146 :
147 5933616 : do ie=nets,nete
148 5194800 : kptr = 4*(kmin-1)
149 5194800 : call ghostunpack(ghostBufFlux, fvm(ie)%se_flux(:,:,:,kmin:kmax),4*kblk,kptr,ie)
150 140998416 : do k=kmin,kmax
151 140259600 : call ghost_flux_unpack(fvm(ie),fvm(ie)%se_flux(:,:,:,k))
152 : end do
153 : enddo
154 :
155 : !call t_stopf('fvm:orthogonal_swept_areas')
156 5933616 : do ie=nets,nete
157 : ! Intel compiler version 2023.0.0 on derecho had significant slowdown on subroutine interface without
158 : ! these pointers.
159 5194800 : fcube => fvm(ie)%c(:,:,:,:)
160 5194800 : spherecentroid => fvm(ie)%spherecentroid(:,1-nhe:nc+nhe,1-nhe:nc+nhe)
161 140998416 : do k=kmin,kmax
162 : !call t_startf('FVM:tracers_reconstruct')
163 : call reconstruction(fcube,nlev,k,&
164 : ctracer(:,:,:,:),irecons_tracer,llimiter,ntrac,&
165 : nc,nhe,nhr,nhc,nht,ns,nhr+(nhe-1),&
166 135064800 : fvm(ie)%jx_min,fvm(ie)%jx_max,fvm(ie)%jy_min,fvm(ie)%jy_max,&
167 : fvm(ie)%cubeboundary,fvm(ie)%halo_interp_weight,fvm(ie)%ibase,&
168 : spherecentroid,&
169 135064800 : fvm(ie)%recons_metrics,fvm(ie)%recons_metrics_integral,&
170 : fvm(ie)%rot_matrix,fvm(ie)%centroid_stretch,&
171 : fvm(ie)%vertex_recons_weights,fvm(ie)%vtx_cart,&
172 405194400 : irecons_tracer_lev(k))
173 : !call t_stopf('FVM:tracers_reconstruct')
174 : !call t_startf('fvm:swept_flux')
175 140259600 : call swept_flux(elem(ie),fvm(ie),k,ctracer,irecons_tracer_lev(k),gsweights,gspts)
176 : !call t_stopf('fvm:swept_flux')
177 : end do
178 : end do
179 : !
180 : !***************************************
181 : !
182 : ! Large Courant number increment
183 : !
184 : !***************************************
185 : !
186 : ! In the jet region the effective Courant number
187 : ! in the cslam trajectory algorithm can be > 1
188 : ! (by up to 20%) in CAM
189 : !
190 : ! We limit the trajectories to < 1 but in this step
191 : ! we do a piecewise constant update for the
192 : ! amount of mass for which the Courant number is >1
193 : !
194 : !
195 738816 : if (large_Courant_incr) then
196 : !call t_startf('fvm:fill_halo_fvm:large_Courant')
197 : !if (kmin_jet<kmin.or.kmax_jet>kmax) then
198 : ! call endrun('ERROR: kmax_jet must be .le. kmax passed to run_consistent_se_cslam')
199 : !end if
200 : ! Determine the extent of the JET that is owned by this thread
201 738816 : ActiveJetThread = threadOwnsVertLevel(hybridnew,kmin_jet) .or. threadOwnsVertLevel(hybridnew,kmax_jet)
202 738816 : kmin_jet_local = max(kmin_jet,kmin)
203 738816 : kmax_jet_local = min(kmax_jet,kmax)
204 738816 : klev = kmax_jet-kmin_jet+1
205 738816 : call fill_halo_fvm(ghostbufQ1,elem,fvm,hybridnew,nets,nete,1,kmin_jet_local,kmax_jet_local,klev,active=ActiveJetThread)
206 : !call t_stopf('fvm:fill_halo_fvm:large_Courant')
207 : !call t_startf('fvm:large_Courant_number_increment')
208 738816 : if(ActiveJetThread) then
209 19948032 : do k=kmin_jet_local,kmax_jet_local !1,nlev
210 155012832 : do ie=nets,nete
211 154274016 : call large_courant_number_increment(fvm(ie),k)
212 : end do
213 : end do
214 : endif
215 : !call t_stopf('fvm:large_Courant_number_increment')
216 : end if
217 :
218 : !call t_startf('fvm:end_of_reconstruct_subroutine')
219 19948032 : do k=kmin,kmax
220 : !
221 : ! convert to mixing ratio
222 : !
223 155012832 : do ie=nets,nete
224 540259200 : do j=1,nc
225 1755842400 : do i=1,nc
226 1620777600 : inv_dp_area(i,j) = 1.0_r8/fvm(ie)%dp_fvm(i,j,k)
227 : end do
228 : end do
229 :
230 540259200 : do itr=1,ntrac
231 1755842400 : do j=1,nc
232 5267527200 : do i=1,nc
233 : ! convert to mixing ratio
234 3646749600 : fvm(ie)%c(i,j,k,itr) = fvm(ie)%c(i,j,k,itr)*inv_dp_area(i,j)
235 : ! remove round-off undershoots
236 4862332800 : fvm(ie)%c(i,j,k,itr) = MAX(fvm(ie)%c(i,j,k,itr),qmin(itr))
237 : end do
238 : end do
239 : end do
240 : !
241 : ! convert to dp and scale back dp
242 : !
243 1755842400 : fvm(ie)%dp_fvm(1:nc,1:nc,k) = fvm(ie)%dp_fvm(1:nc,1:nc,k)*fvm(ie)%dp_ref(k)*fvm(ie)%inv_area_sphere
244 : #ifdef waccm_debug
245 : do j=1,nc
246 : do i=1,nc
247 : fvm(ie)%CSLAM_gamma(i,j,k,1) = MAXVAL(fvm(ie)%CSLAM_gamma(i,j,k,:))
248 : end do
249 : end do
250 : #endif
251 7177643616 : elem(ie)%sub_elem_mass_flux(:,:,:,k)=0
252 : end do
253 : end do
254 : !call t_stopf('fvm:end_of_reconstruct_subroutine')
255 : !$OMP END PARALLEL
256 738816 : call omp_set_nested(.false.)
257 1477632 : end subroutine run_consistent_se_cslam
258 :
259 135064800 : subroutine swept_flux(elem,fvm,ilev,ctracer,irecons_tracer_actual,gsweights,gspts)
260 738816 : use fvm_analytic_mod , only: get_high_order_weights_over_areas
261 : use dimensions_mod, only : kmin_jet,kmax_jet
262 : implicit none
263 : type (element_t) , intent(in) :: elem
264 : type (fvm_struct), intent(inout):: fvm
265 : integer , intent(in) :: ilev, irecons_tracer_actual
266 : real (kind=r8), intent(inout) :: ctracer(irecons_tracer,1-nhe:nc+nhe,1-nhe:nc+nhe,ntrac)
267 : real (kind=r8), dimension(ngpc), intent(in) :: gsweights, gspts
268 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
269 : real (kind=r8) , dimension(0:7 , imin:imax,imin:imax,num_sides) :: displ
270 : integer (kind=r8) , dimension(1:2,11 , imin:imax,imin:imax,num_sides) :: base_vec
271 : real (kind=r8) , dimension(1:2, 6 , imin:imax,imin:imax,num_sides) :: base_vtx
272 : integer , dimension(2,num_area, imin:imax,imin:imax,num_sides) :: idx
273 : real (kind=r8) , dimension(imin:imax,imin:imax,num_sides) :: mass_flux_se
274 : real (kind=r8) , dimension(irecons_tracer,num_area) :: weights
275 : real (kind=r8) :: gamma
276 : integer :: i,j,iside,iarea,iw
277 :
278 : integer, parameter :: num_seg_max=5
279 : REAL(KIND=r8), dimension(2,num_seg_max,num_area) :: x, dx, x_static, dx_static
280 : integer , dimension(num_area) :: num_seg, num_seg_static
281 : REAL(KIND=r8), dimension(2,8) :: x_start, dgam_vec
282 : REAL(KIND=r8) :: gamma_max, displ_first_guess
283 :
284 270129600 : REAL(KIND=r8) :: flux,flux_tracer(ntrac)
285 :
286 : REAL(KIND=r8), dimension(num_area) :: dp_area
287 :
288 : real (kind=r8) :: dp(1-nhc:nc+nhc,1-nhc:nc+nhc)
289 :
290 : logical :: tl1,tl2,tr1,tr2
291 :
292 : integer, dimension(4), parameter :: imin_side = (/1 ,0 ,1 ,1 /)
293 : integer, dimension(4), parameter :: imax_side = (/nc ,nc ,nc ,nc+1/)
294 : integer, dimension(4), parameter :: jmin_side = (/1 ,1 ,0 ,1 /)
295 : integer, dimension(4), parameter :: jmax_side = (/nc+1,nc ,nc ,nc /)
296 :
297 : integer :: iseg, iseg_tmp,flowcase,ii,jj,itr
298 :
299 135064800 : call define_swept_areas(fvm,ilev,displ,base_vec,base_vtx,idx)
300 :
301 7158434400 : mass_flux_se(1:nc,1:nc,1:4) = -elem%sub_elem_mass_flux(1:nc,1:nc,1:4,ilev)
302 540259200 : mass_flux_se(0 ,1:nc,2 ) = elem%sub_elem_mass_flux(1 ,1:nc,4 ,ilev)
303 540259200 : mass_flux_se(nc+1,1:nc,4 ) = elem%sub_elem_mass_flux(nc ,1:nc,2 ,ilev)
304 540259200 : mass_flux_se(1:nc,0 ,3 ) = elem%sub_elem_mass_flux(1:nc,1 ,1 ,ilev)
305 540259200 : mass_flux_se(1:nc,nc+1,1 ) = elem%sub_elem_mass_flux(1:nc,nc ,3 ,ilev)
306 : !
307 : ! prepare for air/tracer update
308 : !
309 : ! dp = fvm%dp_fvm(1-nhe:nc+nhe,1-nhe:nc+nhe,ilev)
310 12290896800 : dp = fvm%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,ilev)
311 1755842400 : fvm%dp_fvm(1:nc,1:nc,ilev) = fvm%dp_fvm(1:nc,1:nc,ilev)*fvm%area_sphere
312 540259200 : do itr=1,ntrac
313 5267527200 : fvm%c(1:nc,1:nc,ilev,itr) = fvm%c(1:nc,1:nc,ilev,itr)*fvm%dp_fvm(1:nc,1:nc,ilev)
314 2971425600 : do iw=1,irecons_tracer_actual
315 2431166400 : ctracer(iw,1-nhe:nc+nhe,1-nhe:nc+nhe,itr)=ctracer(iw,1-nhe:nc+nhe,1-nhe:nc+nhe,itr)*&
316 78202519200 : dp(1-nhe:nc+nhe,1-nhe:nc+nhe)
317 : end do
318 : end do
319 :
320 675324000 : do iside=1,4
321 2566231200 : do j=jmin_side(iside),jmax_side(iside)
322 8914276800 : do i=imin_side(iside),imax_side(iside)
323 : !DO NOT USE MASS_FLUX_SE AS THRESHOLD - THRESHOLD CONDITION MUST BE CONSISTENT WITH
324 : !THE ONE USED IN DEFINE_SWEPT_AREAS
325 : ! if (mass_flux_se(i,j,iside)>eps) then
326 8374017600 : if (fvm%se_flux(i,j,iside,ilev)>eps) then
327 : !
328 : ! || ||
329 : ! tl1 || || tr1
330 : ! || ||
331 : ! =============================
332 : ! || ||
333 : ! tl2 || || tr2
334 : ! || ||
335 : !
336 3241555200 : tl1 = displ(3,i,j,iside)<0.0_r8.and.displ(6,i,j,iside).ge.0.0_r8 !departure point in tl1 quadrant
337 3241555200 : tl2 = displ(6,i,j,iside)<0.0_r8.and.displ(7,i,j,iside) >0.0_r8 !departure point in tl2 quadrant
338 3241555200 : tr1 = displ(2,i,j,iside)<0.0_r8.and.displ(4,i,j,iside).ge.0.0_r8 !departure point in tr1 quadrant
339 3241555200 : tr2 = displ(4,i,j,iside)<0.0_r8.and.displ(5,i,j,iside) >0.0_r8 !departure point in tr2 quadrant
340 :
341 : !
342 : ! pathological cases
343 : !
344 : ! | || || || ||
345 : ! | ||-----------|| ||-----------||
346 : ! | || || || ||
347 : ! ================================ =================================
348 : ! || || | || ||
349 : ! ---------|| || ------|--|| ||
350 : ! || || | || ||
351 : !
352 : ! tl1=tl1.or.tl2
353 : ! tr1=tr1.or.tr2
354 : ! tl1=displ(3,i,j,iside)<0.0_r8.and..not.(tl1.and.tl2)
355 : ! tr1=displ(2,i,j,iside)<0.0_r8.and..not.(tr1.and.tr2)
356 :
357 38898662400 : num_seg=-1; num_seg_static=-1 !initialization
358 3241555200 : if (.not.tl1.and..not.tl2.and..not.tr1.and..not.tr2) then
359 49048163 : flowcase=0
360 : !
361 : ! || || || || || ||
362 : ! || * * || || *----------* |*----------* ||
363 : ! || / \ || || / || || \ ||
364 : ! ||/ \|| ||/ || || \||
365 : ! ============================= ============================= =============================
366 : ! || || || || || ||
367 : !
368 : !
369 : call define_area3_center (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
370 49048163 : num_seg_static,x_start, dgam_vec,fvm%se_flux(i,j,iside,ilev),displ_first_guess)
371 :
372 49048163 : gamma=1.0_r8!fvm%se_flux(i,j,iside,ilev)
373 49048163 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
374 : else
375 3192507037 : if (tl1.and.tr1) then
376 41320433 : flowcase=1
377 : !
378 : !
379 : ! tl1 || || tr1 || || || ||
380 : ! *--||-------------||--* *--||-------------|| ||-------------||--*
381 : ! \ || || / \ || ||\ /|| || /
382 : ! \|| ||/ \|| || \ / || ||/
383 : ! ============================= =========================*=== ==*==========================
384 : ! || || || || || ||
385 : !
386 : call define_area2 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static,&
387 41320433 : num_seg, num_seg_static,x_start, dgam_vec,displ_first_guess)
388 : call define_area3_left_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static,&
389 41320433 : num_seg, num_seg_static,x_start, dgam_vec)
390 : call define_area4 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static,&
391 41320433 : num_seg, num_seg_static,x_start, dgam_vec)
392 41320433 : gamma=1.0_r8
393 41320433 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
394 3151186604 : else if (tl1.and..not.tr1.and..not.tr2) then
395 1517343698 : flowcase=2
396 : !
397 : ! || || || || || ||
398 : ! *--||----------* || /||----------* || *--||-------------*
399 : ! \ || \ || / || \ || \ || ||
400 : ! \|| \|| / || \|| \|| ||
401 : ! ============================= ==*========================== =============================
402 : ! || || || || || ||
403 : !
404 : call define_area2 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, num_seg_static,&
405 1517343698 : x_start, dgam_vec,displ_first_guess)
406 : call define_area3_left(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, num_seg_static,&
407 1517343698 : x_start, dgam_vec)
408 1517343698 : gamma=1.0_r8
409 1517343698 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
410 1633842906 : else if (tr1.and..not.tl1.and..not.tl2) then !displ(3).ge.0.0_r8) then
411 1517586108 : flowcase=3
412 : !
413 : ! || *----------||--* || *----------||\ *-------------||--*
414 : ! || / || / || / || \ || || /
415 : ! ||/ ||/ ||/ || \ || ||/
416 : ! ============================= ==========================*== =============================
417 : ! || || || || || ||
418 : ! || || || || || ||
419 : ! || || || || || ||
420 : !
421 : call define_area3_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, &
422 1517586108 : num_seg_static, x_start, dgam_vec)
423 : call define_area4 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, &
424 1517586108 : num_seg_static, x_start, dgam_vec,displ_first_guess)
425 1517586108 : gamma=1.0_r8
426 1517586108 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
427 116256798 : else if (tl2.and..not.tr1.and..not.tr2) then !displ(2).ge.0.0_r8) then
428 57024381 : flowcase=4
429 : !
430 : ! ||----------* || ||-------------*
431 : ! /|| \ || /|| ||
432 : ! / || \|| / || ||
433 : ! ===/========================= ===/=========================
434 : ! | /|| || | /|| ||
435 : ! |/ || || |/ || ||
436 : ! * || || * || ||
437 : !
438 : call define_area1_area2(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
439 57024381 : num_seg_static,x_start, dgam_vec)
440 : call define_area3_left (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
441 : num_seg_static,&
442 57024381 : x_start, dgam_vec,displ_first_guess)
443 57024381 : gamma = 1.0_r8
444 57024381 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
445 59232417 : else if (tr2.and..not.tl1.and..not.tl2) then !displ(3).ge.0.0_r8) then
446 57014235 : flowcase=5
447 : ! case(5)
448 : !
449 : !
450 : ! || *-----2----||
451 : ! || /1 3||\
452 : ! ||/ 4 || \
453 : ! =============================
454 : ! || ||\ |
455 : ! || || \|
456 : ! || || *
457 : !
458 : call define_area3_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
459 57014235 : num_seg_static,x_start, dgam_vec)
460 : call define_area4_area5(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
461 57014235 : num_seg_static,x_start, dgam_vec,displ_first_guess)
462 57014235 : gamma=1.0_r8
463 57014235 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
464 2218182 : else if (tl2.and.tr1.and..not.tr2) then
465 1078046 : flowcase=6
466 : ! case(6)
467 : !
468 : !
469 : ! ||-------------||--*
470 : ! /|| || /
471 : ! / || ||/
472 : ! ===/=========================
473 : ! | /|| ||
474 : ! |/ || ||
475 : ! * || ||
476 : !
477 : !
478 : call define_area1_area2 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
479 1078046 : num_seg_static,x_start, dgam_vec)
480 : call define_area3_left_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
481 1078046 : num_seg_static,x_start, dgam_vec)
482 : call define_area4 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
483 1078046 : num_seg_static,x_start, dgam_vec,displ_first_guess)
484 :
485 1078046 : gamma=1.0_r8
486 1078046 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
487 1140136 : else if (tr2.and.tl1.and..not.tl2) then
488 1112119 : flowcase=7
489 : ! case(7)
490 : !
491 : !
492 : ! *--||-------------||
493 : ! \ || ||\
494 : ! \|| || \
495 : ! =============================
496 : ! || ||\ |
497 : ! || || \|
498 : ! || || *
499 : !
500 : !
501 : call define_area2 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
502 1112119 : num_seg_static,x_start, dgam_vec,displ_first_guess)
503 : call define_area3_left_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
504 1112119 : num_seg_static,x_start, dgam_vec)
505 : call define_area4_area5 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
506 1112119 : num_seg_static,x_start, dgam_vec)
507 1112119 : gamma = 1.0_r8
508 1112119 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
509 28017 : else if (tl2.and.tr2) then
510 28017 : flowcase=8
511 : ! case(8)
512 : !
513 : !
514 : ! ||-------------||
515 : ! /|| ||\
516 : ! / || || \
517 : ! =============================
518 : ! | /|| ||\ |
519 : ! |/ || || \|
520 : ! * || || *
521 : !
522 : !
523 : !
524 : !
525 : !
526 : call define_area1_area2 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
527 28017 : num_seg_static,x_start, dgam_vec)
528 : call define_area3_left_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
529 28017 : num_seg_static,x_start, dgam_vec)
530 : call define_area4_area5 (i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg,&
531 28017 : num_seg_static,x_start, dgam_vec,displ_first_guess)
532 28017 : gamma = 1.0_r8
533 28017 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
534 : else
535 0 : call endrun('ERROR - unknown flow case')
536 : end if
537 : end if
538 : !
539 : ! iterate to get flux area
540 : !
541 : !call t_startf('fvm:swept_area:get_gamma')
542 19449331200 : do iarea=1,num_area
543 19449331200 : dp_area(iarea) = dp(idx(1,iarea,i,j,iside),idx(2,iarea,i,j,iside))
544 : end do
545 : call get_flux_segments_area_iterate(x,x_static,dx_static,dx,x_start,dgam_vec,num_seg,num_seg_static,&
546 9724665600 : num_seg_max,num_area,dp_area,flowcase,gamma,mass_flux_se(i,j,iside),0.0_r8,gamma_max, &
547 9724665600 : gsweights,gspts,ilev)
548 : !call t_stopf('fvm:swept_area:get_gamma')
549 : !
550 : ! pack segments for high-order weights computation
551 : !
552 19449331200 : do iarea=1,num_area
553 25932441600 : do iseg=1,num_seg_static(iarea)
554 9724665600 : iseg_tmp=num_seg(iarea)+iseg
555 29173996800 : x (:,iseg_tmp,iarea) = x_static (:,iseg,iarea)
556 45381772800 : dx(:,iseg_tmp,iarea) = dx_static(:,iseg,iarea)
557 : end do
558 19449331200 : num_seg(iarea)=num_seg(iarea)+MAX(0,num_seg_static(iarea))
559 : end do
560 : !
561 : ! compute higher-order weights
562 : !
563 : !call t_startf('fvm:swept_area:get_high_order_w')
564 : call get_high_order_weights_over_areas(x,dx,num_seg,num_seg_max,num_area,weights,ngpc,&
565 3241555200 : gsweights, gspts,irecons_tracer)
566 : !call t_stopf('fvm:swept_area:get_high_order_w')
567 : !
568 : !**************************************************
569 : !
570 : ! remap air and tracers
571 : !
572 : !**************************************************
573 : !
574 : !call t_startf('fvm:swept_area:remap')
575 12966220800 : flux=0.0_r8; flux_tracer=0.0_r8
576 19449331200 : do iarea=1,num_area
577 19449331200 : if (num_seg(iarea)>0) then
578 6593885667 : ii=idx(1,iarea,i,j,iside); jj=idx(2,iarea,i,j,iside)
579 6593885667 : flux=flux+weights(1,iarea)*dp(ii,jj)
580 26375542668 : do itr=1,ntrac
581 >14506*10^7 : do iw=1,irecons_tracer_actual
582 >13847*10^7 : flux_tracer(itr) = flux_tracer(itr)+weights(iw,iarea)*ctracer(iw,ii,jj,itr)
583 : end do
584 : end do
585 : end if
586 : end do
587 3241555200 : fvm%se_flux(i,j,iside,ilev) = mass_flux_se(i,j,iside)-flux
588 3241555200 : if (fvm%se_flux(i,j,iside,ilev)>1.0E-13_r8.and.(ilev<kmin_jet.or.ilev>kmax_jet)) then
589 0 : write(iulog,*) "CN excess flux outside of pre-scribed jet region"
590 0 : write(iulog,*) "Increase jet region with kmin_jet and kmax_jet ",&
591 0 : ilev,fvm%se_flux(i,j,iside,ilev),mass_flux_se(i,j,iside),flux,flowcase,&
592 0 : kmin_jet,kmax_jet
593 0 : call endrun('ERROR in CSLAM: local Courant number is > 1; Increase kmin_jet/kmax_jet?')
594 : end if
595 :
596 3241555200 : fvm%dp_fvm(i ,j ,ilev ) = fvm%dp_fvm(i ,j ,ilev )-flux
597 12966220800 : fvm% c(i ,j ,ilev,1:ntrac) = fvm% c(i ,j ,ilev,1:ntrac)-flux_tracer(1:ntrac)
598 : !
599 : ! update flux in nearest neighbor cells
600 : !
601 3241555200 : if (iside==1) then
602 808237483 : fvm%dp_fvm(i,j-1,ilev ) = fvm%dp_fvm(i,j-1,ilev )+flux
603 3232949932 : fvm% c(i,j-1,ilev,1:ntrac) = fvm% c(i,j-1,ilev,1:ntrac)+flux_tracer(1:ntrac)
604 : end if
605 3241555200 : if (iside==2) then
606 1071792916 : fvm%dp_fvm(i+1,j,ilev ) = fvm%dp_fvm(i+1,j,ilev )+flux
607 4287171664 : fvm% c(i+1,j,ilev,1:ntrac) = fvm% c(i+1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
608 : end if
609 3241555200 : if (iside==3) then
610 812540117 : fvm%dp_fvm(i,j+1,ilev ) = fvm%dp_fvm(i,j+1,ilev )+flux
611 3250160468 : fvm% c(i,j+1,ilev,1:ntrac) = fvm% c(i,j+1,ilev,1:ntrac)+flux_tracer(1:ntrac)
612 : end if
613 3241555200 : if (iside==4) then
614 548984684 : fvm%dp_fvm(i-1,j,ilev ) = fvm%dp_fvm(i-1,j,ilev )+flux
615 2195938736 : fvm% c(i-1,j,ilev,1:ntrac) = fvm% c(i-1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
616 : end if
617 : !call t_stopf('fvm:swept_area:remap')
618 : end if
619 : end do
620 : end do
621 : end do
622 135064800 : end subroutine swept_flux
623 :
624 :
625 135064800 : subroutine large_courant_number_increment(fvm,ilev)
626 : implicit none
627 : type (fvm_struct), intent(inout):: fvm
628 : integer , intent(in) :: ilev
629 :
630 : integer, parameter :: num_sides=4, imin= 0, imax=nc+1
631 :
632 : integer, dimension(4), parameter :: imin_side = (/1 ,0 ,1 ,1 /)
633 : integer, dimension(4), parameter :: imax_side = (/nc ,nc ,nc ,nc+1/)
634 : integer, dimension(4), parameter :: jmin_side = (/1 ,1 ,0 ,1 /)
635 : integer, dimension(4), parameter :: jmax_side = (/nc+1,nc ,nc ,nc /)
636 :
637 : integer :: i,j,iside,itr
638 270129600 : real (kind=r8) :: flux,flux_tracer(ntrac)
639 : real (kind=r8), dimension(0:nc+1,0:nc+1) :: inv_dp_area
640 270129600 : real (kind=r8), dimension(0:nc+1,0:nc+1,ntrac):: c_tmp
641 :
642 4187008800 : inv_dp_area=1.0_r8/fvm%dp_fvm(0:nc+1,0:nc+1,ilev)
643 12696091200 : c_tmp = fvm%c(0:nc+1,0:nc+1,ilev,1:ntrac)
644 675324000 : do iside=1,4
645 2566231200 : do j=jmin_side(iside),jmax_side(iside)
646 8914276800 : do i=imin_side(iside),imax_side(iside)
647 8374017600 : if (fvm%se_flux(i,j,iside,ilev)>eps) then
648 : flux = fvm%se_flux(i,j,iside,ilev)
649 : #ifdef waccm_debug
650 : if (i>0.and.j>0.and.i<nc+1.and.j<nc+1) then
651 : fvm%CSLAM_gamma(i,j,ilev,iside) = fvm%CSLAM_gamma(i,j,ilev,iside)+&
652 : fvm%se_flux(i,j,iside,ilev)*inv_dp_area(i,j)
653 : end if
654 : #endif
655 :
656 232724 : do itr=1,ntrac
657 232724 : flux_tracer(itr) = fvm%se_flux(i,j,iside,ilev)*c_tmp(i,j,itr)*inv_dp_area(i,j)
658 : end do
659 58181 : fvm%dp_fvm(i ,j ,ilev ) = fvm%dp_fvm(i ,j ,ilev )-flux
660 232724 : fvm% c(i ,j ,ilev,1:ntrac) = fvm% c(i ,j ,ilev,1:ntrac)-flux_tracer(1:ntrac)
661 : !
662 : ! update flux in nearest neighbor cells
663 : !
664 58181 : if (iside==1) then
665 36298 : fvm%dp_fvm(i,j-1,ilev ) = fvm%dp_fvm(i,j-1,ilev )+flux
666 145192 : fvm% c(i,j-1,ilev,1:ntrac) = fvm% c(i,j-1,ilev,1:ntrac)+flux_tracer(1:ntrac)
667 : end if
668 58181 : if (iside==2) then
669 17602 : fvm%dp_fvm(i+1,j,ilev ) = fvm%dp_fvm(i+1,j,ilev )+flux
670 70408 : fvm% c(i+1,j,ilev,1:ntrac) = fvm% c(i+1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
671 : end if
672 58181 : if (iside==3) then
673 1760 : fvm%dp_fvm(i,j+1,ilev ) = fvm%dp_fvm(i,j+1,ilev )+flux
674 7040 : fvm% c(i,j+1,ilev,1:ntrac) = fvm% c(i,j+1,ilev,1:ntrac)+flux_tracer(1:ntrac)
675 : end if
676 58181 : if (iside==4) then
677 2521 : fvm%dp_fvm(i-1,j,ilev ) = fvm%dp_fvm(i-1,j,ilev )+flux
678 10084 : fvm% c(i-1,j,ilev,1:ntrac) = fvm% c(i-1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
679 : end if
680 : end if
681 : end do
682 : end do
683 : end do
684 135064800 : end subroutine large_courant_number_increment
685 :
686 135064800 : subroutine ghost_flux_unpack(fvm,var)
687 : use control_mod, only : neast, nwest, seast, swest
688 : implicit none
689 : type (fvm_struct), intent(inout) :: fvm
690 : real(kind=r8) :: var(1-nhe:nc+nhe,1-nhe:nc+nhe,4)
691 :
692 : integer :: i,j,ishft
693 : !
694 : ! rotate coordinates if needed
695 : !
696 135064800 : if (fvm%cubeboundary.NE.0) then
697 104450112 : do j=1-nhe,nc+nhe
698 539658912 : do i=1-nhe,nc+nhe
699 435208800 : ishft = NINT(fvm%flux_orient(2,i,j))
700 2263085760 : var(i,j,1:4) = cshift(var(i,j,1:4),shift=ishft)
701 : end do
702 : end do
703 : !
704 : ! non-existent cells in physical space - necessary?
705 : !
706 17408352 : if (fvm%cubeboundary==nwest) then
707 1950936 : var(1-nhe:0,nc+1 :nc+nhe,:) = 0.0_r8
708 17258280 : else if (fvm%cubeboundary==swest) then
709 1950936 : var(1-nhe:0,1-nhe:0 ,:) = 0.0_r8
710 17108208 : else if (fvm%cubeboundary==neast) then
711 1950936 : var(nc+1 :nc+nhe,nc+1 :nc+nhe,:) = 0.0_r8
712 16958136 : else if (fvm%cubeboundary==seast) then
713 1950936 : var(nc+1 :nc+nhe,1-nhe:0,:) = 0.0_r8
714 : end if
715 : end if
716 135064800 : end subroutine ghost_flux_unpack
717 :
718 135064800 : subroutine compute_displacements_for_swept_areas(fvm,cair,k,gsweights,gspts)
719 : use dimensions_mod, only: large_Courant_incr
720 : implicit none
721 : type (fvm_struct), intent(inout) :: fvm
722 : integer, intent(in) :: k
723 : real (kind=r8) :: cair(1-nhc:nc+nhc,1-nhc:nc+nhc) !high-order air density reconstruction
724 : real (kind=r8), dimension(ngpc), intent(in) :: gsweights, gspts
725 : !
726 : ! flux iside 1 flux iside 3 flux iside 2 flux iside 4
727 : !
728 : ! | | | ---1--> | | --2-->| |--1--> |
729 : ! -4----------3- /\ -4----------3- -4----------3- -4----------3- ||
730 : ! | | /||\ |\\\\\\\\\\| || | |\\\\\\| |\\\\\\| |
731 : ! | --2--> | || dv(1) |\\\\\\\\\\| || | |\\\\\\| |\\\\\\| |
732 : ! |----------| || |----------| || dv(3) | |\\\\\\| |\\\\\\| |
733 : ! |\\\\\\\\\\| || | <--2--- | \||/ | |\\\\\\| |\\\\\\| |
734 : ! |\\\\\\\\\\| || | | \/ | |\\\\\\| |\\\\\\| |
735 : ! -1----------2- -1----------2- -1----------2- -1----------2-
736 : ! | <--1-- | | | | <--1--| |<--2--
737 : !
738 : ! / \
739 : ! line-integral <========== =========>
740 : ! from vertex 2 \ dv(2) dv(4)/
741 : ! to 1
742 : !
743 : ! Note vertical
744 : ! lines have
745 : ! zero line-
746 : ! integral!
747 : !
748 : integer :: i,j,iside,ix
749 : integer, parameter :: num_area=1, num_seg_max=2
750 : REAL(KIND=r8), dimension(2,num_seg_max,num_area,4,nc,nc) :: x_static, dx_static
751 : REAL(KIND=r8), dimension(2,num_seg_max,num_area,4,nc,nc) :: x, dx
752 : REAL(KIND=r8), dimension(2,num_seg_max,num_area) :: x_tmp, dx_tmp
753 : integer , dimension( num_area,4 ) :: num_seg, num_seg_static
754 : REAL(KIND=r8), dimension(2,8, 4,nc,nc) :: x_start, dgam_vec
755 : REAL(KIND=r8), dimension(num_area) :: dp_area
756 : integer, dimension(4) :: flowcase
757 : REAL(KIND=r8) :: gamma(4), flux_se
758 :
759 135064800 : num_seg_static(1,1) = 1; num_seg(1,1) = 1; flowcase(1) = -1
760 135064800 : num_seg_static(1,2) = 0; num_seg(1,2) = 2; flowcase(2) = -2
761 135064800 : num_seg_static(1,3) = 1; num_seg(1,3) = 1; flowcase(3) = -1
762 135064800 : num_seg_static(1,4) = 0; num_seg(1,4) = 2; flowcase(4) = -4
763 :
764 540259200 : do j=1,nc
765 1755842400 : do i=1,nc
766 4051944000 : do ix=1,2
767 2431166400 : iside=1;
768 2431166400 : x_static (ix,1,1,iside,i,j) = fvm%vtx_cart(2,ix,i,j)
769 2431166400 : dx_static(ix,1,1,iside,i,j) = fvm%vtx_cart(1,ix,i,j)-fvm%vtx_cart(2,ix,i,j)
770 2431166400 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(1,ix,i,j)
771 2431166400 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(2,ix,i,j)
772 2431166400 : dgam_vec (ix,1, iside,i,j) = fvm%vtx_cart(4,ix,i,j)-fvm%vtx_cart(1,ix,i,j)
773 : !
774 : ! compute first guess
775 : !
776 2431166400 : gamma(iside) = 0.5_r8
777 2431166400 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
778 2431166400 : dx (ix,1,1,iside,i,j) = -dx_static(ix,1,1,iside,i,j)
779 : !
780 : ! side 2
781 : !
782 2431166400 : iside=2;
783 2431166400 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(2,ix,i,j)
784 2431166400 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(3,ix,i,j)
785 2431166400 : dgam_vec (ix,1, iside,i,j) = fvm%vtx_cart(1,ix,i,j)-fvm%vtx_cart(2,ix,i,j)
786 2431166400 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)
787 : !
788 : ! compute first guess - gamma=1
789 : !
790 2431166400 : gamma(iside) = 0.5_r8
791 2431166400 : dx (ix,1,1,iside,i,j) = gamma(iside)*dgam_vec (ix,1, iside,i,j)
792 2431166400 : x (ix,2,1,iside,i,j) = x_start(ix,2,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
793 2431166400 : dx (ix,2,1,iside,i,j) = -gamma(iside)*dgam_vec (ix,1, iside,i,j)
794 : !
795 : ! side 3
796 : !
797 2431166400 : iside=3;
798 2431166400 : x_static (ix,1,1,iside,i,j) = fvm%vtx_cart(4,ix,i,j)
799 2431166400 : dx_static(ix,1,1,iside,i,j) = fvm%vtx_cart(3,ix,i,j)-fvm%vtx_cart(4,ix,i,j)
800 2431166400 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(3,ix,i,j)
801 2431166400 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(4,ix,i,j)
802 2431166400 : dgam_vec (ix,1, iside,i,j) = fvm%vtx_cart(2,ix,i,j)-fvm%vtx_cart(3,ix,i,j)
803 : !
804 : ! compute first guess - gamma(iside)=1
805 : !
806 2431166400 : gamma(iside) = 0.5_r8
807 2431166400 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
808 2431166400 : dx (ix,1,1,iside,i,j) = -dx_static(ix,1,1,iside,i,j)
809 : !
810 : ! side 4
811 : !
812 2431166400 : iside=4;
813 2431166400 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(1,ix,i,j)
814 2431166400 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(4,ix,i,j)
815 2431166400 : dgam_vec (ix,1, iside,i,j) = fvm%vtx_cart(2,ix,i,j)-fvm%vtx_cart(1,ix,i,j)
816 2431166400 : x (ix,2,1,iside,i,j) = x_start(ix,2,iside,i,j)
817 : !
818 : ! compute first guess - gamma(iside)=1
819 : !
820 2431166400 : gamma(iside) = 0.5_r8
821 2431166400 : dx (ix,2,1,iside,i,j) = gamma(iside)*dgam_vec (ix,1, iside,i,j)
822 2431166400 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
823 3646749600 : dx (ix,1,1,iside,i,j) = -gamma(iside)*dgam_vec (ix,1, iside,i,j)
824 : end do
825 : end do
826 : end do
827 :
828 : ! do k=1,nlev
829 540259200 : do j=1,nc
830 1755842400 : do i=1,nc
831 2431166400 : dp_area = cair(i,j)
832 6483110400 : do iside=1,4
833 4862332800 : flux_se = -fvm%se_flux(i,j,iside,k)
834 6077916000 : if (flux_se>eps) then
835 2431166400 : gamma(iside)=0.5_r8
836 : !
837 : ! this copying is necessary since get_flux_segments_area_iterate change x and dx
838 : !
839 15802749276 : x_tmp (:,1:num_seg(1,iside),:)=x (:,1:num_seg(1,iside),:,iside,i,j)
840 15802749276 : dx_tmp(:,1:num_seg(1,iside),:)=dx(:,1:num_seg(1,iside),:,iside,i,j)
841 : call get_flux_segments_area_iterate(&
842 : x_tmp(:,:,:),x_static(:,:,:,iside,i,j),dx_static(:,:,:,iside,i,j),dx_tmp(:,:,:),&
843 : x_start(:,:,iside,i,j),dgam_vec(:,:,iside,i,j),num_seg(:,iside),num_seg_static(:,iside),&
844 : num_seg_max,num_area,dp_area,flowcase(iside),gamma(iside),flux_se,0.0_r8,1.0_r8, &
845 2431166400 : gsweights,gspts,k)
846 7293499200 : fvm%se_flux(i,j,iside,k) = ABS(SUM(gamma(iside)*dgam_vec(:,1,iside,i,j)))
847 : #ifdef waccm_debug
848 : fvm%CSLAM_gamma(i,j,k,iside) = gamma(iside)
849 : #endif
850 2431166400 : if (gamma(iside)>1_r8) then
851 0 : if (.not.large_Courant_incr) then
852 0 : write(iulog,*) 'ERROR in CSLAM: local Courant number is >1: gamma=',gamma(iside),' k=',k
853 0 : call endrun('ERROR in CSLAM: local Courant number is > 1; set namelist se_large_Courant_incr=.true. ')
854 : endif
855 0 : gamma(iside)=1.0_r8-eps
856 : end if
857 : else
858 2431166400 : fvm%se_flux(i,j,iside,k) = 0.0_r8
859 : #ifdef waccm_debug
860 : fvm%CSLAM_gamma(i,j,k,iside) = 0.0_r8
861 : #endif
862 : end if
863 : enddo
864 : end do
865 : end do
866 : ! end do
867 135064800 : end subroutine compute_displacements_for_swept_areas
868 :
869 :
870 :
871 5672721600 : subroutine get_flux_segments_area_iterate(x,x_static,dx_static,dx,x_start,dgam_vec,num_seg,num_seg_static,&
872 5672721600 : num_seg_max,num_area,c,flow_case,gamma,flux,gamma_min,gamma_max,gsweights,gspts,ilev)
873 : implicit none
874 : integer , intent(in) :: num_area, num_seg_max
875 : REAL(KIND=r8), dimension(2,num_seg_max,num_area), intent(in) :: x_static, dx_static
876 : REAL(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx
877 : integer , dimension(num_area ), intent(in) :: num_seg, num_seg_static
878 : REAL(KIND=r8), dimension(2,8) , intent(in) :: x_start, dgam_vec
879 : REAL(KIND=r8) , intent(inout) :: gamma
880 : REAL(KIND=r8) , intent(in) :: flux,gamma_min,gamma_max
881 : integer , intent(in) :: flow_case,ilev
882 :
883 : real (kind=r8), dimension(num_area) , intent(in) :: c
884 : real (kind=r8), dimension(ngpc) , intent(in) :: gsweights, gspts
885 :
886 : real (kind=r8) :: flux_static
887 11345443200 : real (kind=r8) :: weight_area(num_area), xtmp(2), xtmp2(2)
888 : real (kind=r8) :: gamma1, gamma2, gamma3, dgamma, f1, f2
889 :
890 : real (kind=r8), dimension( ngpc ) :: xq,yq
891 : real (kind=r8), dimension( ngpc,1) :: F !linear
892 :
893 5672721600 : real (kind=r8) :: xq2,xq2i, rho, rhoi, yrh, w_static(num_area)
894 :
895 : integer :: iseg,iarea,iter,ipt
896 : integer, parameter :: iter_max=40
897 : logical :: lexit_after_one_more_iteration
898 :
899 5672721600 : lexit_after_one_more_iteration = .false.
900 : !
901 : ! compute static line-integrals (not necessary to recompute them for every iteration)
902 : !
903 5672721600 : flux_static = 0.0_r8
904 24311664000 : w_static = 0.0_r8
905 24311664000 : weight_area = 0.0_r8
906 24311664000 : do iarea=1,num_area
907 29579135308 : do iseg=1,num_seg_static(iarea)
908 :
909 : !rck vector directive needed here
910 : !DIR$ SIMD
911 43760771632 : do ipt=1,ngpc
912 32820578724 : xq(ipt) = x_static(1,iseg,iarea)+dx_static(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
913 32820578724 : yq(ipt) = x_static(2,iseg,iarea)+dx_static(2,iseg,iarea)*gspts(ipt)
914 43760771632 : F(ipt,1) = yq(ipt)/(SQRT(1.0_r8+xq(ipt)*xq(ipt) + yq(ipt)*yq(ipt))*(1.0_r8+xq(ipt)*xq(ipt)))! potential ! potential
915 : enddo
916 62399714032 : weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx_static(1,iseg,iarea) !integral
917 : end do
918 18638942400 : w_static(iarea)= weight_area(iarea)
919 24311664000 : flux_static = flux_static+weight_area(iarea)*c(iarea) !add to swept flux
920 : end do
921 : !
922 : ! initilization
923 : !
924 5672721600 : gamma1=0.0_r8; f1=-flux ! zero flux guess 1
925 : !
926 : ! compute flux integrals of first guess passed to subroutine
927 : !
928 5672721600 : gamma2=gamma
929 5672721600 : f2 = 0.0_r8
930 24311664000 : weight_area=w_static
931 24311664000 : do iarea=1,num_area
932 35633342656 : do iseg=1,num_seg(iarea)
933 : !rck vector directive needed here
934 : !DIR$ SIMD
935 67977601024 : do ipt=1,ngpc
936 50983200768 : xq(ipt) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
937 50983200768 : yq(ipt) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(ipt)
938 50983200768 : xq2 = xq(ipt)*xq(ipt)
939 50983200768 : xq2i = 1.0_r8/(1.0_r8+xq2)
940 50983200768 : rho = SQRT(1.0_r8+xq2+yq(ipt)*yq(ipt))
941 50983200768 : rhoi = 1.0_r8/rho
942 50983200768 : yrh = yq(ipt)*rhoi
943 67977601024 : F(ipt,1) = yrh*xq2i
944 : enddo
945 86616543424 : weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx(1,iseg,iarea)! integral
946 : end do
947 24311664000 : f2 = f2+weight_area(iarea)*c(iarea)
948 : end do
949 5672721600 : f2 = f2-flux !integral error
950 5672721600 : iter=0
951 5672721600 : if (abs(f2-f1)<eps) then
952 : !
953 : ! in case the first guess is converged
954 : !
955 : return
956 : end if
957 :
958 :
959 5672721600 : dgamma=(gamma2-gamma1)*f2/(f2-f1);
960 5672721600 : gamma3 = gamma2-dgamma; ! Newton "guess" for gamma
961 5672721600 : gamma1 = gamma2; f1 = f2; gamma2 = gamma3; ! prepare for iteration
962 16054882493 : do iter=1,iter_max
963 : !
964 : ! update vertex location: flow_case dependent to avoid many zero operations
965 : !
966 17341712049 : select case(flow_case)
967 : case(-4)
968 1286829556 : iarea=1
969 3860488668 : dx (:,2,1) = gamma3*dgam_vec (:,1)
970 3860488668 : x (:,1,1) = x_start(:,1)+gamma3*dgam_vec(:,1)
971 3860488668 : dx (:,1,1) = -gamma3*dgam_vec (:,1)
972 :
973 : case(-2)
974 2523517074 : iarea=1
975 7570551222 : dx (:,1,iarea) = gamma3*dgam_vec (:,1)
976 7570551222 : x (:,2,iarea) = x_start(:,2)+gamma3*dgam_vec(:,1)
977 7570551222 : dx (:,2,iarea) = -gamma3*dgam_vec (:,1)
978 : case(-1)
979 : !
980 : ! to compute first-guess perpendicular displacements for iside=1
981 : !
982 3799076835 : iarea=1
983 11397230505 : x (:,1,iarea) = x_start(:,1)+gamma3*dgam_vec(:,1)
984 11397230505 : dx (:,1,iarea) = -dx_static(:,1,iarea)
985 11397230505 : x (:,2,iarea) = x_start(:,2)+gamma3*dgam_vec(:,1)
986 11397230505 : dx (:,2,iarea) = x_start(:,2)-x(:,2,iarea)
987 : case(0)
988 134566768 : iarea=3
989 403700304 : xtmp = x_start(:,1)+gamma3*dgam_vec(:,1)
990 403700304 : dx (:,1,iarea) = xtmp(: )-x(:,1,iarea) !dynamic - line 2
991 403700304 : x (:,2,iarea) = xtmp(: ) !dynamic - line 3
992 403700304 : dx (:,2,iarea) = x_static(:,2,iarea)-x(:,2,iarea) !dynamic - line 3
993 : case(1)
994 88457026 : iarea=2
995 265371078 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
996 265371078 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
997 265371078 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
998 265371078 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
999 :
1000 88457026 : iarea = 3
1001 265371078 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1002 265371078 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1003 265371078 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1004 265371078 : x (:,2,iarea) = xtmp (:) !dynamic
1005 265371078 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1006 265371078 : x (:,3,iarea) = xtmp2(:) !dynamic
1007 265371078 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1008 :
1009 88457026 : iarea = 4
1010 265371078 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1011 265371078 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1012 265371078 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1013 265371078 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1014 : case(2)
1015 3929822858 : iarea=2
1016 11789468574 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1017 11789468574 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1018 11789468574 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
1019 11789468574 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
1020 :
1021 3929822858 : iarea=3
1022 11789468574 : xtmp(: ) = x_start(:,4)+gamma3*dgam_vec(:,4)!
1023 11789468574 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1024 11789468574 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1025 11789468574 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1026 : case(3)
1027 3920199421 : iarea = 3
1028 11760598263 : xtmp (: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1029 11760598263 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1030 11760598263 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1031 11760598263 : dx (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
1032 :
1033 3920199421 : iarea = 4
1034 11760598263 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1035 11760598263 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1036 11760598263 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1037 11760598263 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1038 : case(4)
1039 183838405 : iarea = 1
1040 551515215 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1041 551515215 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1042 551515215 : x (:,2,iarea) = xtmp(:) !dynamic
1043 551515215 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1044 :
1045 183838405 : iarea = 2
1046 551515215 : xtmp (: ) = x_start(:,2)+gamma3*dgam_vec(:,2)
1047 551515215 : xtmp2 (: ) = x_start(:,3)+gamma3*dgam_vec(:,3)
1048 :
1049 551515215 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1050 :
1051 551515215 : x (:,2,iarea) = xtmp (:) !dynamic
1052 551515215 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1053 :
1054 551515215 : x (:,3,iarea) = xtmp2(:) !dynamic
1055 551515215 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1056 :
1057 183838405 : iarea = 3
1058 551515215 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1059 551515215 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1060 551515215 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1061 551515215 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1062 : case(5)
1063 184157424 : iarea = 3
1064 552472272 : xtmp (: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1065 552472272 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1066 552472272 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1067 552472272 : dx (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
1068 :
1069 184157424 : iarea = 4
1070 552472272 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1071 552472272 : xtmp2 (: ) = x_start(:,7)+gamma3*dgam_vec(:,7)
1072 :
1073 552472272 : dx(:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1074 552472272 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1075 552472272 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic - line 2
1076 552472272 : x (:,3,iarea) = xtmp2(:) !dynamic -line 1
1077 552472272 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic - line 1
1078 :
1079 184157424 : iarea = 5
1080 552472272 : xtmp (: ) = x_start(:,8)+gamma3*dgam_vec(:,8)
1081 :
1082 552472272 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1083 552472272 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1084 552472272 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1085 : case(6)
1086 2149657 : iarea = 1
1087 6448971 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1088 6448971 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1089 6448971 : x (:,2,iarea) = xtmp(:) !dynamic
1090 6448971 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1091 :
1092 2149657 : iarea = 2
1093 6448971 : xtmp (: ) = x_start(:,2)+gamma3*dgam_vec(:,2)
1094 6448971 : xtmp2 (: ) = x_start(:,3)+gamma3*dgam_vec(:,3)
1095 :
1096 6448971 : dx(:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1097 6448971 : x (:,2,iarea) = xtmp (:) !dynamic
1098 6448971 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1099 6448971 : x (:,3,iarea) = xtmp2(:) !dynamic
1100 6448971 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1101 :
1102 2149657 : iarea = 3
1103 6448971 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1104 6448971 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1105 6448971 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1106 6448971 : x (:,2,iarea) = xtmp (:) !dynamic
1107 6448971 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1108 6448971 : x (:,3,iarea) = xtmp2(:) !dynamic
1109 6448971 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1110 :
1111 2149657 : iarea = 4
1112 6448971 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1113 6448971 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1114 6448971 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1115 6448971 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1116 : case(7)
1117 2211760 : iarea=2
1118 6635280 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1119 6635280 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1120 6635280 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
1121 6635280 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
1122 :
1123 2211760 : iarea = 3
1124 6635280 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1125 6635280 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1126 6635280 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1127 6635280 : x (:,2,iarea) = xtmp (:) !dynamic
1128 6635280 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1129 6635280 : x (:,3,iarea) = xtmp2(:) !dynamic
1130 6635280 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1131 :
1132 2211760 : iarea = 4
1133 6635280 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1134 6635280 : xtmp2 (: ) = x_start(:,7)+gamma3*dgam_vec(:,7)
1135 :
1136 6635280 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1137 6635280 : x (:,2,iarea) = xtmp(:) !dynamic
1138 6635280 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1139 6635280 : x (:,3,iarea) = xtmp2(:) !dynamic
1140 6635280 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1141 :
1142 2211760 : iarea = 5
1143 6635280 : xtmp (: ) = x_start(:,8)+gamma3*dgam_vec(:,8)
1144 6635280 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1145 6635280 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1146 6635280 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1147 : case(8)
1148 55709 : iarea = 1
1149 167127 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1150 167127 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1151 167127 : x (:,2,iarea) = xtmp(:) !dynamic
1152 167127 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1153 :
1154 55709 : iarea = 2
1155 167127 : xtmp (: ) = x_start(:,2)+gamma3*dgam_vec(:,2)
1156 167127 : xtmp2 (: ) = x_start(:,3)+gamma3*dgam_vec(:,3)
1157 :
1158 167127 : dx(:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1159 167127 : x (:,2,iarea) = xtmp (:) !dynamic
1160 167127 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1161 167127 : x (:,3,iarea) = xtmp2(:) !dynamic
1162 167127 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1163 :
1164 55709 : iarea = 3
1165 167127 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1166 167127 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1167 167127 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1168 167127 : x (:,2,iarea) = xtmp (:) !dynamic
1169 167127 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1170 167127 : x (:,3,iarea) = xtmp2(:) !dynamic
1171 167127 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1172 :
1173 55709 : iarea = 4
1174 167127 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1175 167127 : xtmp2 (: ) = x_start(:,7)+gamma3*dgam_vec(:,7)
1176 :
1177 167127 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1178 167127 : x (:,2,iarea) = xtmp(:) !dynamic
1179 167127 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1180 167127 : x (:,3,iarea) = xtmp2(:) !dynamic
1181 167127 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1182 :
1183 55709 : iarea = 5
1184 167127 : xtmp (: ) = x_start(:,8)+gamma3*dgam_vec(:,8)
1185 167127 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1186 167127 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1187 167127 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1188 : case default
1189 16054882493 : call endrun('flow case not defined in get_flux_segments_area_iterate')
1190 : end select
1191 : !
1192 : ! compute flux integral
1193 : !
1194 65891601098 : f2 = 0.0_r8
1195 65891601098 : weight_area=w_static
1196 65891601098 : do iarea=1,num_area
1197 96165219724 : do iseg=1,num_seg(iarea)
1198 : !rck vector directive needed here
1199 : !DIR$ SIMD
1200 >18531*10^7 : do ipt=1,ngpc
1201 :
1202 >13898*10^7 : xq(ipt) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
1203 >13898*10^7 : yq(ipt) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(ipt)
1204 :
1205 >13898*10^7 : xq2 = xq(ipt)*xq(ipt)
1206 >13898*10^7 : xq2i = 1.0_r8/(1.0_r8+xq2)
1207 >13898*10^7 : rho = SQRT(1.0_r8+xq2+yq(ipt)*yq(ipt))
1208 >13898*10^7 : rhoi = 1.0_r8/rho
1209 >13898*10^7 : yrh = yq(ipt)*rhoi
1210 >18531*10^7 : F(ipt,1) = yrh*xq2i
1211 : end do
1212 >23515*10^7 : weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx(1,iseg,iarea)! integral
1213 : end do
1214 65891601098 : f2 = f2+weight_area(iarea)*c(iarea)
1215 : end do
1216 16054882493 : f2 = f2-flux !integral error
1217 :
1218 : !
1219 : ! uncommented logic leads to noise in PS in FKESSLER at element boundary
1220 : !
1221 : ! if (ABS(f2)<eps.or.ABS((gamma2-gamma1)*f2)<eps.or.lexit_after_one_more_iteration) then
1222 : !
1223 21727604093 : if (ABS(f2)<eps.or.lexit_after_one_more_iteration) then
1224 5672815930 : gamma=gamma3
1225 5672815930 : if (gamma>gamma_max) then
1226 94330 : lexit_after_one_more_iteration=.true.
1227 94330 : gamma=gamma_max
1228 94330 : gamma3=gamma_max
1229 : else
1230 : exit
1231 : end if
1232 : else
1233 : !
1234 : ! Newton increment
1235 : !
1236 10382066563 : if (abs(f2-f1)<eps) then
1237 : !
1238 : ! if entering here abs(f2)>eps and abs(f1)>eps but abs(f2-f1)<eps
1239 : !
1240 26 : dgamma=-0.5_r8*(gamma2-gamma1)
1241 26 : lexit_after_one_more_iteration=.true.
1242 : else
1243 10382066537 : dgamma=(gamma2-gamma1)*f2/(f2-f1)
1244 : endif
1245 10382066563 : if (ABS(dgamma)>eps) then
1246 10382066563 : gamma3 = gamma2-dgamma;
1247 : else
1248 : !
1249 : ! dgamma set to minimum displacement to avoid f2-f1=0
1250 : !
1251 0 : gamma3=gamma2-SIGN(1.0_r8,dgamma)*eps
1252 : end if
1253 10382066563 : gamma3=MAX(gamma3,gamma_min)
1254 : !
1255 : ! prepare for next iteration
1256 : !
1257 10382066563 : gamma1 = gamma2; f1 = f2; gamma2 = gamma3;
1258 : endif
1259 : end do
1260 5672721600 : if (iter>iter_max) write(iulog,*) "WARNING: iteration not converged",&
1261 0 : ABS(f2),flux,gamma1,gamma2,gamma3,ilev
1262 : end subroutine get_flux_segments_area_iterate
1263 :
1264 135064800 : subroutine define_swept_areas(fvm,ilev,displ,base_vec,base_vtx,idx)
1265 : use control_mod, only : neast, nwest, seast, swest
1266 : implicit none
1267 : type (fvm_struct), intent(inout) :: fvm
1268 : integer , intent(in) :: ilev
1269 :
1270 :
1271 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1272 : real (kind=r8) , dimension(0:7 , imin:imax,imin:imax,num_sides), intent(out) :: displ
1273 : integer (kind=r8) , dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(out) :: base_vec
1274 : real (kind=r8) , dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(out) :: base_vtx
1275 : integer , dimension(2,num_area, imin:imax,imin:imax,num_sides), intent(out) :: idx
1276 :
1277 : real (kind=r8) :: flux_sum (0:nc+1,0:nc+1,2)
1278 : integer :: degenerate (1:nc+1,1:nc+1 )
1279 : integer :: circular_flow(1:nc+1,1:nc+1 )
1280 : integer :: illcond (1:nc+1,1:nc+1)
1281 : integer :: ib,i,j,sgn, iside, iarea
1282 :
1283 : !
1284 : ! set where reconstruction function is as a function of area and side
1285 : !
1286 : integer, dimension(num_area*4), parameter :: idx_shift_tmp = (/-1,-1, 0, 1, 1,& !iside=1
1287 : 1, 0, 0, 0, 1,& !iside=2
1288 : 1, 1, 0,-1,-1,& !iside=3
1289 : -1, 0, 0, 0,-1/) !iside=4
1290 :
1291 : integer, dimension(num_area*4), parameter :: idy_shift_tmp = (/-1, 0, 0, 0,-1,& !iside=1
1292 : -1,-1, 0, 1, 1,& !iside=2
1293 : 1, 0, 0, 0, 1,& !iside=3
1294 : 1, 1, 0,-1,-1/) !iside=4
1295 :
1296 : integer, dimension(num_area,4), parameter :: idx_shift = RESHAPE(idx_shift_tmp,(/num_area,4/))
1297 : integer, dimension(num_area,4), parameter :: idy_shift = RESHAPE(idy_shift_tmp,(/num_area,4/))
1298 :
1299 : integer, dimension(4), parameter :: iside_m1 = (/4,1,2,3/)
1300 : integer, dimension(4), parameter :: iside_p1 = (/2,3,4,1/)
1301 : integer, dimension(4), parameter :: iside_p2 = (/3,4,1,2/)
1302 : integer, dimension(4), parameter :: iside_p3 = (/4,1,2,3/)
1303 :
1304 : integer, dimension(4), parameter :: imin_side = (/1 ,0 ,1 ,1 /)
1305 : integer, dimension(4), parameter :: imax_side = (/nc ,nc ,nc ,nc+1/)
1306 : integer, dimension(4), parameter :: jmin_side = (/1 ,1 ,0 ,1 /)
1307 : integer, dimension(4), parameter :: jmax_side = (/nc+1,nc ,nc ,nc /)
1308 :
1309 :
1310 :
1311 : integer :: iur,jur,ilr,jlr,iul,jul,ill,jll
1312 :
1313 135064800 : ib = fvm%cubeboundary
1314 3376620000 : flux_sum(0:nc+1,1:nc+1,1) = fvm%se_flux(0:nc+1,0:nc ,3,ilev)-fvm%se_flux(0:nc+1,1:nc+1,1,ilev)
1315 3511684800 : flux_sum(1:nc+1,0:nc+1,2) = fvm%se_flux(0:nc ,0:nc+1,2,ilev)-fvm%se_flux(1:nc+1,0:nc+1,4,ilev)
1316 :
1317 : !
1318 : ! Degenerate case ("two departure points")
1319 : !
1320 : ! || | || no change in this situation || no change in this situation
1321 : ! || | || ||
1322 : ! ||-------- ||---------- ||----------
1323 : ! || | || ||
1324 : ! ======================= ======================= =====================
1325 : ! | || | || ||
1326 : ! -----|---|| ------|---|| ---------||
1327 : ! | || | || ||
1328 : ! | || | || ||
1329 : !
1330 : !
1331 2836360800 : where (flux_sum(0:nc,1:nc+1,1)*flux_sum(1:nc+1,1:nc+1,1)<0.0_r8.and.flux_sum(1:nc+1,0:nc,2)*flux_sum(1:nc+1,1:nc+1,2)<0.0_r8)
1332 : degenerate(:,:) = 0
1333 : elsewhere
1334 : degenerate(:,:) = 1
1335 : end where
1336 :
1337 135064800 : if (ib>0) then
1338 17408352 : if (ib==swest) degenerate(1 ,1 ) = 1
1339 17258280 : if (ib==nwest) degenerate(1 ,nc+1) = 1
1340 17258280 : if (ib==neast) degenerate(nc+1,nc+1) = 1
1341 17258280 : if (ib==seast) degenerate(nc+1,1 ) = 1
1342 : end if
1343 :
1344 675324000 : do j=1,nc+1
1345 2836360800 : do i=1,nc+1
1346 4862332800 : do sgn=-1,1,2
1347 : if (&
1348 17288294400 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1349 23771404800 : sgn*flux_sum(i ,j,1)>0.0_r8.and.sgn*flux_sum(i,j ,2)<0.0_r8) then
1350 2606426 : circular_flow(i,j) = 0
1351 : else
1352 4319467174 : circular_flow(i,j) = 1
1353 : end if
1354 : end do
1355 : end do
1356 : end do
1357 : !
1358 : ! wrap around corners
1359 : !
1360 135064800 : if (ib==nwest) then
1361 150072 : flux_sum(0,nc+1,1) = fvm%se_flux(0,nc,3,ilev)-fvm%se_flux(1,nc+1,4,ilev)
1362 150072 : flux_sum(1,nc+1,2) = fvm%se_flux(0,nc,3,ilev)-fvm%se_flux(1,nc+1,4,ilev)
1363 :
1364 150072 : i=1;j=nc+1;
1365 150072 : circular_flow(i,j) = 1
1366 450216 : do sgn=-1,1,2
1367 : if (&
1368 : sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1369 450216 : sgn*flux_sum(i ,j,1)>0.0_r8.and.sgn*flux_sum(i,j ,2)<0.0_r8) then
1370 48 : circular_flow(i,j) = 0
1371 : end if
1372 : end do
1373 134914728 : else if (ib==swest) then
1374 150072 : flux_sum(0,1,1) = fvm%se_flux(1,0,4,ilev)-fvm%se_flux(0,1,1,ilev)
1375 150072 : flux_sum(1,0,2) = fvm%se_flux(0,1,1,ilev)-fvm%se_flux(1,0,4,ilev)
1376 150072 : i=1;j=1;
1377 150072 : circular_flow(i,j) = 1
1378 450216 : do sgn=-1,1,2
1379 : if (&
1380 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.&
1381 450216 : sgn*flux_sum(i ,j,1)>0.0_r8.and.sgn*flux_sum(i,j ,2)<0.0_r8) then
1382 72 : circular_flow(i,j) = 0
1383 : end if
1384 : end do
1385 134764656 : else if (ib==neast) then
1386 150072 : flux_sum(nc+1,nc+1,1) = fvm%se_flux(nc+1,nc,3,ilev)-fvm%se_flux(nc,nc+1,2,ilev)
1387 150072 : flux_sum(nc+1,nc+1,2) = fvm%se_flux(nc,nc+1,2,ilev)-fvm%se_flux(nc+1,nc,3,ilev)
1388 150072 : i=nc+1;j=nc+1;
1389 150072 : circular_flow(i,j) = 1
1390 450216 : do sgn=-1,1,2
1391 : if (&
1392 300144 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1393 150072 : sgn*flux_sum(i,j ,2)<0.0_r8) then
1394 41 : circular_flow(i,j) = 0
1395 : end if
1396 : end do
1397 134614584 : else if (ib==seast) then
1398 150072 : flux_sum(nc+1,1 ,1) = fvm%se_flux(nc,0,2,ilev)-fvm%se_flux(nc+1,1,1,ilev)
1399 150072 : flux_sum(nc+1,0 ,2) = fvm%se_flux(nc,0,2,ilev)-fvm%se_flux(nc+1,1,1,ilev)
1400 150072 : i=nc+1;j=1;
1401 150072 : circular_flow(i,j) = 1
1402 450216 : do sgn=-1,1,2
1403 : if (&
1404 300144 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1405 150072 : sgn*flux_sum(i,j ,2)<0.0_r8) then
1406 64 : circular_flow(i,j) = 0
1407 : end if
1408 : end do
1409 : end if
1410 2836360800 : illcond = circular_flow*degenerate
1411 : !
1412 : !
1413 : !
1414 : !
1415 675324000 : do iside=1,4
1416 2566231200 : do j=jmin_side(iside),jmax_side(iside)
1417 8914276800 : do i=imin_side(iside),imax_side(iside)
1418 8374017600 : if (fvm%se_flux(i,j,iside,ilev)>eps) then
1419 3241555200 : iur = i+idx_shift(4,iside); jur = j+idy_shift(4,iside) !(i,j) index of upper right quadrant
1420 3241555200 : ilr = i+idx_shift(5,iside); jlr = j+idy_shift(5,iside) !(i,j) index of lower left quadrant
1421 3241555200 : iul = i+idx_shift(2,iside); jul = j+idy_shift(2,iside) !(i,j) index of upper right quadrant
1422 3241555200 : ill = i+idx_shift(1,iside); jll = j+idy_shift(1,iside) !(i,j) index of lower left quadrant
1423 :
1424 : !iside=1
1425 3241555200 : if (iside==1) then
1426 808237483 : displ(0,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i,j) !center left
1427 808237483 : displ(1,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i+1,j) !center right
1428 808237483 : displ(2,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j) !c2
1429 808237483 : displ(3,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j) !c3
1430 808237483 : displ(4,i,j,iside) = -flux_sum (i+1,j ,1)*illcond(i+1,j) !r1
1431 808237483 : displ(5,i,j,iside) = -flux_sum (i+1,j-1,2)*illcond(i+1,j) !r2
1432 808237483 : displ(6,i,j,iside) = -flux_sum (i-1,j ,1)*illcond(i ,j) !l1
1433 808237483 : displ(7,i,j,iside) = flux_sum (i ,j-1,2)*illcond(i ,j) !l2
1434 :
1435 : end if
1436 3241555200 : if (iside==2) then
1437 : !iside=2
1438 1071792916 : displ(0,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j ) !center left
1439 1071792916 : displ(1,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j+1) !center right
1440 1071792916 : displ(2,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i+1,j+1) !c2
1441 1071792916 : displ(3,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i+1,j ) !c3
1442 1071792916 : displ(4,i,j,iside) = flux_sum (i+1,j+1,2)*illcond(i+1,j+1) !r1
1443 1071792916 : displ(5,i,j,iside) = -flux_sum (i+1,j+1,1)*illcond(i+1,j+1) !r2
1444 1071792916 : displ(6,i,j,iside) = flux_sum (i+1,j-1,2)*illcond(i+1,j) !l1
1445 1071792916 : displ(7,i,j,iside) = flux_sum (i+1,j ,1)*illcond(i+1,j) !l2
1446 : end if
1447 : !iside=3
1448 3241555200 : if (iside==3) then
1449 812540117 : displ(0,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i+1,j+1) !center left
1450 812540117 : displ(1,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i ,j+1) !center right
1451 812540117 : displ(2,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j+1) !c2
1452 812540117 : displ(3,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j+1) !c3
1453 812540117 : displ(4,i,j,iside) = flux_sum (i-1,j+1,1)*illcond(i ,j+1) !r1
1454 812540117 : displ(5,i,j,iside) = flux_sum (i ,j+1,2)*illcond(i ,j+1) !r2
1455 812540117 : displ(6,i,j,iside) = flux_sum (i+1,j+1,1)*illcond(i+1,j+1) !l1
1456 812540117 : displ(7,i,j,iside) = -flux_sum (i+1,j+1,2)*illcond(i+1,j+1) !l2
1457 : end if
1458 3241555200 : if (iside==4) then
1459 : !iside=4
1460 548984684 : displ(0,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j+1) !center left
1461 548984684 : displ(1,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j ) !center right
1462 548984684 : displ(2,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i ,j ) !c2
1463 548984684 : displ(3,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i ,j+1) !c3
1464 548984684 : displ(4,i,j,iside) = -flux_sum (i ,j-1,2)*illcond(i ,j ) !r1
1465 548984684 : displ(5,i,j,iside) = flux_sum (i-1,j ,1)*illcond(i ,j ) !r2
1466 548984684 : displ(6,i,j,iside) = -flux_sum (i ,j+1,2)*illcond(i ,j+1) !l1
1467 548984684 : displ(7,i,j,iside) = -flux_sum (i-1,j+1,1)*illcond(i ,j+1) !l2
1468 : end if
1469 :
1470 9724665600 : base_vtx(:,1,i,j,iside) = fvm%vtx_cart(iside,:,i ,j ) !vertex center left
1471 9724665600 : base_vtx(:,2,i,j,iside) = fvm%vtx_cart(iside_p1(iside),:,i ,j ) !vertex center right
1472 9724665600 : base_vtx(:,3,i,j,iside) = fvm%vtx_cart(iside,:,iur,jur ) !vertex upper right
1473 9724665600 : base_vtx(:,4,i,j,iside) = fvm%vtx_cart(iside_p3(iside),:,ilr,jlr) !vertex lower right
1474 9724665600 : base_vtx(:,5,i,j,iside) = fvm%vtx_cart(iside_p1(iside),:,iul,jul) !vertex upper left
1475 9724665600 : base_vtx(:,6,i,j,iside) = fvm%vtx_cart(iside_p2(iside),:,ill,jll) !vertex lower left
1476 :
1477 9724665600 : base_vec(:, 1,i,j,iside) = fvm%flux_vec (:,i ,j ,iside ) !vector center
1478 9724665600 : base_vec(:, 2,i,j,iside) = fvm%flux_vec (:,i ,j ,iside_p1(iside)) !vector center right
1479 9724665600 : base_vec(:, 3,i,j,iside) = fvm%flux_vec (:,i ,j ,iside_p3(iside)) !vector center left
1480 9724665600 : base_vec(:, 4,i,j,iside) = fvm%flux_vec (:,iur,jur,iside ) !vector upper right 1
1481 9724665600 : base_vec(:, 5,i,j,iside) = fvm%flux_vec (:,iur,jur,iside_p3(iside)) !vector upper right 2
1482 9724665600 : base_vec(:, 6,i,j,iside) = fvm%flux_vec (:,ilr,jlr,iside_p3(iside)) !vector lower right 1
1483 9724665600 : base_vec(:, 7,i,j,iside) = fvm%flux_vec (:,ilr,jlr,iside_p2(iside)) !vector lower right 2
1484 9724665600 : base_vec(:, 8,i,j,iside) = fvm%flux_vec (:,iul,jul,iside ) !vector upper left 1
1485 9724665600 : base_vec(:, 9,i,j,iside) = fvm%flux_vec (:,iul,jul,iside_p1(iside)) !vector upper left 2
1486 9724665600 : base_vec(:,10,i,j,iside) = fvm%flux_vec (:,ill,jll,iside_p1(iside)) !vector lower left 1
1487 9724665600 : base_vec(:,11,i,j,iside) = fvm%flux_vec (:,ill,jll,iside_p2(iside)) !vector lower left 2
1488 :
1489 19449331200 : do iarea=1,5
1490 16207776000 : idx(1,iarea,i,j,iside) = i+idx_shift(iarea,iside)
1491 19449331200 : idx(2,iarea,i,j,iside) = j+idy_shift(iarea,iside)
1492 : end do
1493 : else
1494 29173996800 : displ(:,i,j,iside) = 9D99!for debugging
1495 : end if
1496 : end do
1497 : end do
1498 : end do
1499 : !
1500 : ! wrap around corners here
1501 : !
1502 :
1503 135064800 : end subroutine define_swept_areas
1504 :
1505 :
1506 : !
1507 : ! Notation conventions used in define_area subroutines
1508 : !
1509 : !
1510 : !
1511 : ! ^ ||---> ^ <---|| ^
1512 : ! /|\ || 3 /|\ 2 || /|\
1513 : ! | 6 || 1 | || | 4
1514 : ! | || | || |
1515 : ! =================================
1516 : ! || ||
1517 : ! || ||
1518 : ! 7 || || 5
1519 : ! <---|| ||--->
1520 : !
1521 :
1522 58130444 : subroutine define_area1_area2(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, num_seg_static,&
1523 : x_start, dgam_vec)
1524 : implicit none
1525 : integer, intent(in) :: i,j,iside
1526 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1527 : real (kind=r8) , dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1528 : integer (kind=r8) , dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1529 : real (kind=r8) , dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1530 : integer, parameter :: num_seg_max=5
1531 : REAL(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1532 : integer , dimension(num_area) , intent(inout) :: num_seg, num_seg_static
1533 : REAL(KIND=r8), dimension(2,8) , intent(inout):: x_start, dgam_vec
1534 :
1535 :
1536 : real (kind=r8) , dimension(2,3) :: xdep !departure points
1537 : real (kind=r8) :: gamma
1538 : integer :: iarea
1539 :
1540 :
1541 : REAL(KIND=r8) :: xtmp(2),xtmp2(2)
1542 : !
1543 : !
1544 : ! ||----- ||
1545 : ! /|| ||
1546 : ! / || ||
1547 : ! ===X=========================
1548 : ! | /|| ||
1549 : ! |/ || ||
1550 : ! * || ||
1551 : !
1552 : !
1553 : ! crossing X
1554 174391332 : if (SUM(ABS(base_vec(:,9,i,j,iside))).NE.0) then
1555 58028860 : gamma = displ(0,i,j,iside)*displ(7,i,j,iside)/(displ(0,i,j,iside)-displ(6,i,j,iside))
1556 : ! gamma = MAX(MIN(gamma,displ(7,i,j,iside),-displ(3,i,j,iside)),0.0_r8)!MWR manuscript
1557 58028860 : gamma = MAX(MIN(gamma,displ(7,i,j,iside),-0.25_r8*displ(3,i,j,iside)),0.0_r8)
1558 : else
1559 : !
1560 : ! corner case
1561 : !
1562 101584 : gamma=displ(0,i,j,iside)
1563 : end if
1564 :
1565 :
1566 174391332 : xdep (:,1) = base_vtx(:, 6,i,j,iside)+displ(7,i,j,iside)*base_vec(:,10,i,j,iside)-displ(6,i,j,iside)*base_vec(:,11,i,j,iside)
1567 174391332 : x_start (:,1) = base_vtx(:, 6,i,j,iside)
1568 174391332 : dgam_vec(:,1) = base_vec(:,10,i,j,iside)*gamma
1569 :
1570 174391332 : xdep(:,2) = base_vtx(:,2,i,j,iside)+displ(1,i,j,iside)*base_vec(:, 1,i,j,iside)+displ(2,i,j,iside)*base_vec(:, 2,i,j,iside)
1571 :
1572 58130444 : iarea = 1
1573 58130444 : num_seg (iarea) = 2
1574 58130444 : num_seg_static(iarea) = 1
1575 :
1576 174391332 : x_static (:,1,iarea) = base_vtx(:,6,i,j,iside) !static
1577 174391332 : dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static
1578 :
1579 174391332 : xtmp(: ) = x_start(:,1)+dgam_vec(:,1)
1580 174391332 : x (:,1,iarea) = xdep(:,1) !static
1581 174391332 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1582 :
1583 174391332 : x (:,2,iarea) = xtmp(:) !dynamic
1584 174391332 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1585 : !
1586 : !
1587 : !
1588 58130444 : iarea = 2
1589 58130444 : num_seg (iarea) = 3
1590 :
1591 174391332 : x_start (:,2) = base_vtx(:,5,i,j,iside)
1592 174391332 : dgam_vec(:,2) = base_vec(:,9,i,j,iside)*gamma
1593 174391332 : xtmp (: ) = x_start(:,2)+dgam_vec(:,2)
1594 :
1595 174391332 : x_start (:,3) = base_vtx(:,5,i,j,iside)
1596 174391332 : dgam_vec(:,3) = base_vec(:,8,i,j,iside)*displ(0,i,j,iside)
1597 174391332 : xtmp2 (: ) = x_start(:,3)+dgam_vec(:,3)
1598 :
1599 174391332 : x (:,1,iarea) = base_vtx(:,5,i,j,iside) !static
1600 174391332 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1601 :
1602 174391332 : x (:,2,iarea) = xtmp (:) !dynamic
1603 174391332 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1604 :
1605 174391332 : x (:,3,iarea) = xtmp2(:) !dynamic
1606 174391332 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1607 58130444 : end subroutine define_area1_area2
1608 :
1609 :
1610 1559776250 : subroutine define_area2(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, num_seg_static,x_start, dgam_vec,&
1611 : displ_first_guess)
1612 : implicit none
1613 : integer, intent(in) :: i,j,iside
1614 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1615 : real (kind=r8) , dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1616 : integer (kind=r8) , dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1617 : real (kind=r8) , dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1618 : integer, parameter :: num_seg_max=5
1619 : REAL(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1620 : integer , dimension(num_area) , intent(inout) :: num_seg, num_seg_static
1621 : REAL(KIND=r8), dimension(2,8) , intent(inout):: x_start, dgam_vec
1622 :
1623 :
1624 : real (kind=r8) , dimension(2,3) :: xdep !departure points
1625 : real (kind=r8), optional, intent(out) :: displ_first_guess
1626 : real (kind=r8) :: gamma
1627 : integer :: iarea
1628 :
1629 :
1630 : REAL(KIND=r8) :: xtmp(2)
1631 : ! *: xdep(:,1)
1632 : ! x: xtmp
1633 : !
1634 : ! 2 || ||
1635 : ! *--x ||
1636 : ! 1\3||1 ||
1637 : ! \|| ||
1638 : ! =============================
1639 : ! || ||
1640 : !
1641 : !
1642 : ! compute departure points (xdep(1) is left; xdep(3) is right and xdep(2) is midway
1643 : !
1644 : xdep(:,1) = base_vtx(:,5,i,j,iside)+&
1645 4679328750 : MAX(0.0_r8,displ(6,i,j,iside))*base_vec(:,8,i,j,iside)-displ(3,i,j,iside)*base_vec(:,9,i,j,iside)
1646 4679328750 : x_start (:,1) = base_vtx(:,5,i,j,iside)
1647 1559776250 : gamma = displ(0,i,j,iside)
1648 4679328750 : dgam_vec(:,1) = base_vec(:,8,i,j,iside)*gamma
1649 1559776250 : if (present(displ_first_guess)) displ_first_guess = gamma
1650 :
1651 1559776250 : iarea = 2
1652 1559776250 : num_seg (iarea) = 2
1653 1559776250 : num_seg_static(iarea) = 1
1654 :
1655 4679328750 : x_static (:,1,iarea) = base_vtx(:,5,i,j,iside) !static - line 1
1656 4679328750 : dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static - line 1
1657 :
1658 4679328750 : xtmp (: ) = x_start(:,1)+dgam_vec(:,1)
1659 4679328750 : x (:,1,iarea) = xdep(:,1) !static - line 2
1660 4679328750 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1661 :
1662 4679328750 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
1663 4679328750 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
1664 1559776250 : end subroutine define_area2
1665 :
1666 :
1667 1574368079 : subroutine define_area3_left(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, &
1668 : num_seg, num_seg_static,x_start, dgam_vec,displ_first_guess)
1669 : implicit none
1670 : integer, intent(in) :: i,j,iside
1671 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1672 : real (kind=r8) , dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1673 : integer (kind=r8) , dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1674 : real (kind=r8) , dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1675 : integer, parameter :: num_seg_max=5
1676 : REAL(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1677 : integer , dimension(num_area) , intent(inout) :: num_seg, num_seg_static
1678 : REAL(KIND=r8), dimension(2,8) , intent(inout):: x_start, dgam_vec
1679 : real (kind=r8), optional, intent(out) :: displ_first_guess
1680 :
1681 : real (kind=r8) , dimension(2,3) :: xdep !departure points
1682 : real (kind=r8) :: gamma
1683 : integer :: iarea
1684 :
1685 :
1686 : REAL(KIND=r8) :: xtmp(2)
1687 :
1688 : ! iarea = 3
1689 : !-------------------------------------------------------------------------------------------
1690 : !
1691 : ! xtmp xdep(2)
1692 : ! |x-----2------* ||
1693 : ! || \ ||
1694 : ! |1 3 ||
1695 : ! || \||
1696 : ! ===========4==============
1697 : !
1698 : !
1699 1574368079 : xdep(:,2) = base_vtx(:,2,i,j,iside)+displ(1,i,j,iside)*base_vec(:,1,i,j,iside)&
1700 6297472316 : +MAX(0.0_r8,displ(2,i,j,iside))*base_vec(:,2,i,j,iside)
1701 4723104237 : x_start (:,4) = base_vtx(:,1,i,j,iside)
1702 1574368079 : gamma = displ(0,i,j,iside)
1703 4723104237 : dgam_vec(:,4) = base_vec(:,1,i,j,iside)*gamma
1704 4723104237 : xtmp (: ) = x_start(:,4)+dgam_vec(:,4)
1705 :
1706 1574368079 : if (present(displ_first_guess)) displ_first_guess = gamma
1707 :
1708 1574368079 : iarea = 3
1709 1574368079 : num_seg (iarea) = 2
1710 1574368079 : num_seg_static(iarea) = 2
1711 :
1712 4723104237 : x_static (:,1,iarea) = xdep(:,2) !static - line 3
1713 4723104237 : dx_static(:,1,iarea) = base_vtx(:,2,i,j,iside)-xdep(:,2) !static - line 3
1714 :
1715 4723104237 : x_static (:,2,iarea) = base_vtx(:,2,i,j,iside) !static - line 4
1716 4723104237 : dx_static(:,2,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 4
1717 :
1718 4723104237 : x (:,1,iarea) = base_vtx(:,1,i,j,iside) !static - line 1
1719 4723104237 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1720 :
1721 4723104237 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1722 4723104237 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1723 1574368079 : end subroutine define_area3_left
1724 :
1725 1574600343 : subroutine define_area3_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, &
1726 : num_seg_static,x_start, dgam_vec)
1727 : implicit none
1728 : integer, intent(in) :: i,j,iside
1729 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1730 : real (kind=r8) , dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1731 : integer (kind=r8) , dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1732 : real (kind=r8) , dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1733 : integer, parameter :: num_seg_max=5
1734 : REAL(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1735 : integer , dimension(num_area) , intent(inout) :: num_seg, num_seg_static
1736 : REAL(KIND=r8), dimension(2,8) , intent(inout):: x_start, dgam_vec
1737 :
1738 :
1739 : real (kind=r8) , dimension(2,3) :: xdep !departure points
1740 : real (kind=r8) :: gamma
1741 : integer :: iarea
1742 :
1743 : REAL(KIND=r8) :: xtmp(2)
1744 : !
1745 : !
1746 : ! || *-----2----||\
1747 : ! || /1 3|| \
1748 : ! ||/ 4 ||
1749 : ! =============================
1750 : ! || ||
1751 : ! || ||
1752 : ! || ||
1753 : !
1754 1574600343 : xdep(:,1) = base_vtx(:,1,i,j,iside)+displ(0,i,j,iside)*base_vec(:,1,i,j,iside)&
1755 6298401372 : +MAX(0.0_r8,displ(3,i,j,iside))*base_vec(:,3,i,j,iside)
1756 4723801029 : x_start (:,5) = base_vtx(:,2,i,j,iside)
1757 1574600343 : gamma = displ(1,i,j,iside)
1758 4723801029 : dgam_vec(:,5) = base_vec(:,1,i,j,iside)*gamma
1759 4723801029 : xtmp (: ) = x_start(:,5)+dgam_vec(:,5)
1760 :
1761 1574600343 : iarea = 3
1762 1574600343 : num_seg (iarea) = 2
1763 1574600343 : num_seg_static(iarea) = 2
1764 :
1765 4723801029 : x_static (:,1,iarea) = base_vtx(:,1,i,j,iside) !static - line 1
1766 4723801029 : dx_static(:,1,iarea) = xdep(:,1)-base_vtx(:,1,i,j,iside) !static - line 1
1767 :
1768 4723801029 : x_static (:,2,iarea) = base_vtx(:,2,i,j,iside) !static - line 4
1769 4723801029 : dx_static(:,2,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 4
1770 :
1771 4723801029 : x (:,1,iarea) = xdep(:,1) !static - line 2
1772 4723801029 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1773 :
1774 4723801029 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1775 4723801029 : dx (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
1776 1574600343 : end subroutine define_area3_right
1777 :
1778 :
1779 43538615 : subroutine define_area3_left_right(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, &
1780 : num_seg_static,x_start, dgam_vec)
1781 : implicit none
1782 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1783 : integer, parameter :: num_seg_max=5
1784 : integer, intent(in) :: i,j,iside
1785 : real (kind=r8), dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout):: displ
1786 : integer (kind=r8), dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout):: base_vec
1787 : real (kind=r8), dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout):: base_vtx
1788 : real(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout):: x, dx, x_static, dx_static
1789 : integer, dimension(num_area), intent(inout):: num_seg, num_seg_static
1790 : real(KIND=r8), dimension(2,8), intent(inout):: x_start, dgam_vec
1791 :
1792 : real (kind=r8) :: gamma
1793 : integer :: iarea
1794 : real(KIND=r8) :: xtmp(2),xtmp2(2)
1795 : !
1796 : ! ||-------------||
1797 : ! /|| ||\
1798 : ! || ||
1799 : ! =============================
1800 : ! || ||
1801 : ! || ||
1802 : ! || ||
1803 : !
1804 130615845 : x_start (:,4) = base_vtx(:,1,i,j,iside)
1805 130615845 : x_start (:,5) = base_vtx(:,2,i,j,iside)
1806 43538615 : gamma = displ(0,i,j,iside)
1807 130615845 : dgam_vec(:,4) = base_vec(:,1,i,j,iside)*gamma
1808 130615845 : dgam_vec(:,5) = base_vec(:,1,i,j,iside)*gamma
1809 130615845 : xtmp (: ) = x_start(:,4)+dgam_vec(:,4)
1810 130615845 : xtmp2 (: ) = x_start(:,5)+dgam_vec(:,5)
1811 :
1812 43538615 : iarea = 3
1813 43538615 : num_seg (iarea) = 3
1814 43538615 : num_seg_static(iarea) = 1
1815 :
1816 130615845 : x_static (:,1,iarea) = base_vtx(:,2,i,j,iside) !static
1817 130615845 : dx_static(:,1,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static
1818 :
1819 130615845 : x (:,1,iarea) = base_vtx(:,1,i,j,iside) !static
1820 130615845 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1821 :
1822 130615845 : x (:,2,iarea) = xtmp (:) !dynamic
1823 130615845 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1824 :
1825 130615845 : x (:,3,iarea) = xtmp2(:) !dynamic
1826 130615845 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1827 43538615 : end subroutine define_area3_left_right
1828 :
1829 58154371 : subroutine define_area4_area5(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, &
1830 : num_seg_static,x_start, dgam_vec,displ_first_guess)
1831 : implicit none
1832 : integer, intent(in) :: i,j,iside
1833 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1834 : integer, parameter :: num_seg_max=5
1835 : real (kind=r8), dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1836 : integer (kind=r8), dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1837 : real (kind=r8), dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1838 : real(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1839 : integer, dimension(num_area), intent(inout) :: num_seg, num_seg_static
1840 : real(KIND=r8), dimension(2,8), intent(inout) :: x_start, dgam_vec
1841 : real(KIND=r8), optional, intent(out) :: displ_first_guess
1842 :
1843 :
1844 : real (kind=r8) , dimension(2,3) :: xdep !departure points
1845 : real (kind=r8) :: gamma
1846 : integer :: iarea
1847 :
1848 : real(KIND=r8) :: xtmp(2),xtmp2(2)
1849 : !
1850 : ! || --------||
1851 : ! || ||\
1852 : ! || || \
1853 : ! =============================
1854 : ! || ||\ |
1855 : ! || || \|
1856 : ! || || *
1857 : !
1858 : !
1859 : ! iarea = 4
1860 : !
1861 58154371 : iarea = 4
1862 58154371 : num_seg (iarea) = 3
1863 :
1864 174463113 : if (SUM(ABS(base_vec(:,5,i,j,iside))).NE.0) then
1865 58055934 : gamma = displ(1,i,j,iside)*displ(5,i,j,iside)/(displ(1,i,j,iside)-displ(4,i,j,iside))
1866 : ! gamma = MAX(MIN(gamma,displ(5,i,j,iside),-displ(2,i,j,iside)),0.0_r8)!MWR manuscript
1867 58055934 : gamma = MAX(MIN(gamma,displ(5,i,j,iside),-0.25_r8*displ(2,i,j,iside)),0.0_r8)
1868 : else
1869 : !
1870 : ! corner case
1871 : !
1872 98437 : gamma = displ(1,i,j,iside)
1873 : end if
1874 :
1875 58154371 : if (present(displ_first_guess)) displ_first_guess = displ(1,i,j,iside)
1876 :
1877 174463113 : x_start (:,6) = base_vtx(:,3,i,j,iside)
1878 174463113 : dgam_vec(:,6) = base_vec(:,4,i,j,iside)*displ(1,i,j,iside)
1879 174463113 : xtmp (: ) = x_start(:,6)+dgam_vec(:,6)
1880 174463113 : x_start (:,7) = base_vtx(:,3,i,j,iside)
1881 174463113 : dgam_vec(:,7) = base_vec(:,5,i,j,iside)*gamma
1882 174463113 : xtmp2 (: ) = x_start(:,7)+dgam_vec(:,7)
1883 :
1884 174463113 : x (:,1,iarea) = base_vtx(:,3,i,j,iside)!static -line 1
1885 174463113 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1886 :
1887 174463113 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1888 174463113 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic - line 2
1889 :
1890 174463113 : x (:,3,iarea) = xtmp2(:) !static -line 1
1891 174463113 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic - line 1
1892 : !
1893 : !iarea = 5
1894 : !
1895 : xdep(:,1) = base_vtx(:,4,i,j,iside)+displ(5,i,j,iside)*base_vec(:,6,i,j,iside)&
1896 174463113 : -displ(4,i,j,iside)*base_vec(:,7,i,j,iside)
1897 174463113 : x_start (:,8) = base_vtx(:,4,i,j,iside)
1898 174463113 : dgam_vec(:,8) = base_vec(:,6,i,j,iside)*gamma
1899 174463113 : xtmp (: ) = x_start(:,8)+dgam_vec(:,8)
1900 :
1901 58154371 : iarea = 5
1902 58154371 : num_seg (iarea) = 2
1903 58154371 : num_seg_static(iarea) = 1
1904 :
1905 174463113 : x (:,1,iarea) = base_vtx(:,4,i,j,iside)!static -line 1
1906 174463113 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1907 :
1908 174463113 : x_static (:,1,iarea) = xdep(:,1) !static - line 1
1909 174463113 : dx_static(:,1,iarea) = x(:,1,iarea)-x_static(:,1,iarea) !static - line 1
1910 :
1911 174463113 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1912 174463113 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1913 58154371 : end subroutine define_area4_area5
1914 :
1915 :
1916 1559984587 : subroutine define_area4(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, &
1917 : num_seg_static,x_start, dgam_vec,displ_first_guess)
1918 : implicit none
1919 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1920 : integer, parameter :: num_seg_max=5
1921 : integer, intent(in) :: i,j,iside
1922 : real (kind=r8), dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1923 : integer (kind=r8), dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1924 : real (kind=r8), dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1925 :
1926 : real(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1927 : integer, dimension(num_area) , intent(inout) :: num_seg, num_seg_static
1928 : real(KIND=r8), dimension(2,8) , intent(inout) :: x_start, dgam_vec
1929 : real(KIND=r8), optional, intent(out) :: displ_first_guess
1930 :
1931 :
1932 :
1933 : real (kind=r8), dimension(2,3) :: xdep !departure points
1934 : real (kind=r8) :: gamma
1935 : integer :: iarea
1936 : real(KIND=r8) :: xtmp(2)
1937 :
1938 1559984587 : iarea = 4
1939 1559984587 : num_seg (iarea) = 2
1940 1559984587 : num_seg_static(iarea) = 1
1941 :
1942 1559984587 : xdep(:,1) = base_vtx(:,3,i,j,iside)+MAX(0.0_r8,displ(4,i,j,iside))*base_vec(:,4,i,j,iside)&
1943 6239938348 : -displ(2,i,j,iside)*base_vec(:,5,i,j,iside)
1944 4679953761 : x_start (:,6) = base_vtx(:,3,i,j,iside)
1945 1559984587 : gamma = displ(1,i,j,iside)
1946 4679953761 : dgam_vec(:,6) = base_vec(:,4,i,j,iside)*gamma
1947 4679953761 : xtmp (: ) = x_start(:,6)+dgam_vec(:,6)
1948 :
1949 1559984587 : if (present(displ_first_guess)) displ_first_guess = gamma
1950 :
1951 4679953761 : x_static (:,1,iarea) = xdep(:,1) !static
1952 4679953761 : dx_static(:,1,iarea) = base_vtx(:,3,i,j,iside)-xdep(:,1) !static
1953 :
1954 4679953761 : x (:,1,iarea) = base_vtx(:,3,i,j,iside) !static - line 2
1955 4679953761 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1956 :
1957 4679953761 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1958 4679953761 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1959 1559984587 : end subroutine define_area4
1960 :
1961 49048163 : subroutine define_area3_center(i,j,iside,displ,base_vec,base_vtx,x, dx, x_static, dx_static, num_seg, num_seg_static,&
1962 : x_start, dgam_vec,se_flux_center,displ_first_guess)
1963 : implicit none
1964 : integer, intent(in) :: i,j,iside
1965 : integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
1966 : integer, parameter :: num_seg_max=5
1967 : real (kind=r8), dimension(0:7 , imin:imax,imin:imax,num_sides), intent(inout) :: displ
1968 : integer (kind=r8), dimension(1:2,11 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vec
1969 : real (kind=r8), dimension(1:2, 6 , imin:imax,imin:imax,num_sides), intent(inout) :: base_vtx
1970 :
1971 : real(KIND=r8), dimension(2,num_seg_max,num_area), intent(inout) :: x, dx, x_static, dx_static
1972 : integer, dimension(num_area), intent(inout) :: num_seg, num_seg_static
1973 : real(KIND=r8), dimension(2,8), intent(inout) :: x_start, dgam_vec
1974 : real(KIND=r8) , intent(in ) :: se_flux_center
1975 : real(KIND=r8), optional, intent(out) :: displ_first_guess
1976 :
1977 : real (kind=r8) , dimension(2,3) :: xdep !departure points
1978 : real (kind=r8) :: gamma
1979 : integer :: iarea
1980 : !
1981 : ! xdep(2)
1982 : ! ______X______
1983 : ! || / \ ||
1984 : ! || *--/ \--* ||
1985 : ! || /xdep(1) xdep(3)\ ||
1986 : ! ||/ \||
1987 : ! ========================================
1988 : ! || ||
1989 : !
1990 : !
1991 : ! compute departure points (xdep(1) is left; xdep(3) is right and xdep(2) is midway
1992 : !
1993 :
1994 : xdep(:,1) = base_vtx(:,1,i,j,iside)+&
1995 147144489 : displ(0,i,j,iside)*base_vec(:,1,i,j,iside)+displ(3,i,j,iside)*base_vec(:,3,i,j,iside)
1996 : xdep(:,3) = base_vtx(:,2,i,j,iside)+&
1997 147144489 : displ(1,i,j,iside)*base_vec(:,1,i,j,iside)+displ(2,i,j,iside)*base_vec(:,2,i,j,iside)
1998 147144489 : xdep(:,2) = 0.5_r8*(xdep(:,1)+xdep(:,3))
1999 :
2000 49048163 : gamma= se_flux_center
2001 : x_start(:,1) = ABS(base_vec(:,3,i,j,iside))*((xdep(:,2)-base_vtx(:,1,i,j,iside)))+&
2002 147144489 : base_vtx(:,1,i,j,iside) !xdep(2) - midway between departure points projected to side 1
2003 :
2004 147144489 : dgam_vec(:,1) = gamma*base_vec(:,1,i,j,iside)
2005 :
2006 49048163 : if (present(displ_first_guess)) displ_first_guess = gamma
2007 :
2008 147144489 : xdep(:,2) = x_start(:,1)+dgam_vec(:,1)
2009 49048163 : iarea = 3
2010 49048163 : num_seg (iarea) = 2
2011 49048163 : num_seg_static(iarea) = 3
2012 :
2013 : ! ______X______
2014 : ! || 2 / \ 3 ||
2015 : ! || *--/ \--* ||
2016 : ! || / \ ||
2017 : ! ||/ 1 5 4\||
2018 : ! ========================================
2019 : ! || ||
2020 : !
2021 147144489 : x_static (:,1,iarea) = base_vtx(:,1,i,j,iside) !static - line 1
2022 147144489 : dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static - line 1
2023 :
2024 147144489 : x (:,1,iarea) = xdep(:,1) !static - line 2
2025 147144489 : dx (:,1,iarea) = xdep(:,2)-x(:,1,iarea) !dynamic - line 2
2026 :
2027 147144489 : x (:,2,iarea) = xdep(:,2) !dynamic - line 3
2028 147144489 : dx (:,2,iarea) = xdep(:,3)-x(:,2,iarea) !dynamic - line 3
2029 :
2030 147144489 : x_static (:,2,iarea) = xdep(:,3) !static - line 4
2031 147144489 : dx_static(:,2,iarea) = base_vtx(:,2,i,j,iside)-x_static(:,2,iarea)!static - line 4
2032 :
2033 147144489 : x_static (:,3,iarea) = base_vtx(:,2,i,j,iside) !static - line 5
2034 147144489 : dx_static(:,3,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 5
2035 :
2036 49048163 : end subroutine define_area3_center
2037 : end module fvm_consistent_se_cslam
|
