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 8865792 : 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 488311200 : do k=kmin,kmax
113 25605169200 : elem(ie)%sub_elem_mass_flux(:,:,:,k) = dt_fvm*elem(ie)%sub_elem_mass_flux(:,:,:,k)*fvm(ie)%dp_ref_inverse(k)
114 6285708000 : 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 63076416 : do q=1,ntrac
119 57142800 : kptr = kptr + nlev
120 62337600 : 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 488311200 : do k=kmin,kmax
130 25610364000 : 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 62337600 : do q=1,ntrac
135 57142800 : kptr = kptr + nlev
136 62337600 : call ghostunpack(ghostbufQnhc, fvm(ie)%c(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax,q),kblk,kptr,ie)
137 : enddo
138 488311200 : do k=kmin,kmax
139 488311200 : 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 489050016 : do k=kmin,kmax
151 488311200 : 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 489050016 : 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 483116400 : 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 483116400 : 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 1449349200 : irecons_tracer_lev(k))
173 : !call t_stopf('FVM:tracers_reconstruct')
174 : !call t_startf('fvm:swept_flux')
175 488311200 : 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 69448704 : do k=kmin_jet_local,kmax_jet_local !1,nlev
210 552565104 : do ie=nets,nete
211 551826288 : 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 69448704 : do k=kmin,kmax
220 : !
221 : ! convert to mixing ratio
222 : !
223 552565104 : do ie=nets,nete
224 1932465600 : do j=1,nc
225 6280513200 : do i=1,nc
226 5797396800 : inv_dp_area(i,j) = 1.0_r8/fvm(ie)%dp_fvm(i,j,k)
227 : end do
228 : end do
229 :
230 5797396800 : do itr=1,ntrac
231 21740238000 : do j=1,nc
232 69085645200 : do i=1,nc
233 : ! convert to mixing ratio
234 47828523600 : 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 63771364800 : 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 6280513200 : 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 25673879088 : 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 483116400 : 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 966232800 : 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 483116400 : call define_swept_areas(fvm,ilev,displ,base_vec,base_vtx,idx)
300 :
301 25605169200 : mass_flux_se(1:nc,1:nc,1:4) = -elem%sub_elem_mass_flux(1:nc,1:nc,1:4,ilev)
302 1932465600 : mass_flux_se(0 ,1:nc,2 ) = elem%sub_elem_mass_flux(1 ,1:nc,4 ,ilev)
303 1932465600 : mass_flux_se(nc+1,1:nc,4 ) = elem%sub_elem_mass_flux(nc ,1:nc,2 ,ilev)
304 1932465600 : mass_flux_se(1:nc,0 ,3 ) = elem%sub_elem_mass_flux(1:nc,1 ,1 ,ilev)
305 1932465600 : 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 43963592400 : dp = fvm%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,ilev)
311 6280513200 : fvm%dp_fvm(1:nc,1:nc,ilev) = fvm%dp_fvm(1:nc,1:nc,ilev)*fvm%area_sphere
312 5797396800 : do itr=1,ntrac
313 69085645200 : 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 37683079200 : do iw=1,irecons_tracer_actual
315 31885682400 : ctracer(iw,1-nhe:nc+nhe,1-nhe:nc+nhe,itr)=ctracer(iw,1-nhe:nc+nhe,1-nhe:nc+nhe,itr)*&
316 >10256*10^8 : dp(1-nhe:nc+nhe,1-nhe:nc+nhe)
317 : end do
318 : end do
319 :
320 2415582000 : do iside=1,4
321 9179211600 : do j=jmin_side(iside),jmax_side(iside)
322 31885682400 : 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 29953216800 : if (fvm%se_flux(i,j,iside,ilev)>eps) then
327 : !
328 : ! || ||
329 : ! tl1 || || tr1
330 : ! || ||
331 : ! =============================
332 : ! || ||
333 : ! tl2 || || tr2
334 : ! || ||
335 : !
336 11594793592 : 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 11594793592 : tl2 = displ(6,i,j,iside)<0.0_r8.and.displ(7,i,j,iside) >0.0_r8 !departure point in tl2 quadrant
338 11594793592 : 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 11594793592 : 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 >13913*10^7 : num_seg=-1; num_seg_static=-1 !initialization
358 11594793592 : if (.not.tl1.and..not.tl2.and..not.tr1.and..not.tr2) then
359 212177348 : 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 212177348 : num_seg_static,x_start, dgam_vec,fvm%se_flux(i,j,iside,ilev),displ_first_guess)
371 :
372 212177348 : gamma=1.0_r8!fvm%se_flux(i,j,iside,ilev)
373 212177348 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
374 : else
375 11382616244 : if (tl1.and.tr1) then
376 173890477 : 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 173890477 : 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 173890477 : 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 173890477 : num_seg, num_seg_static,x_start, dgam_vec)
392 173890477 : gamma=1.0_r8
393 173890477 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
394 11208725767 : else if (tl1.and..not.tr1.and..not.tr2) then
395 5381012119 : 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 5381012119 : 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 5381012119 : x_start, dgam_vec)
408 5381012119 : gamma=1.0_r8
409 5381012119 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
410 5827713648 : else if (tr1.and..not.tl1.and..not.tl2) then !displ(3).ge.0.0_r8) then
411 5378587637 : 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 5378587637 : 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 5378587637 : num_seg_static, x_start, dgam_vec,displ_first_guess)
425 5378587637 : gamma=1.0_r8
426 5378587637 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
427 449126011 : else if (tl2.and..not.tr1.and..not.tr2) then !displ(2).ge.0.0_r8) then
428 217727352 : 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 217727352 : 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 217727352 : x_start, dgam_vec,displ_first_guess)
443 217727352 : gamma = 1.0_r8
444 217727352 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
445 231398659 : else if (tr2.and..not.tl1.and..not.tl2) then !displ(3).ge.0.0_r8) then
446 219036528 : 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 219036528 : 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 219036528 : num_seg_static,x_start, dgam_vec,displ_first_guess)
462 219036528 : gamma=1.0_r8
463 219036528 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
464 12362131 : else if (tl2.and.tr1.and..not.tr2) then
465 6146877 : 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 6146877 : 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 6146877 : 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 6146877 : num_seg_static,x_start, dgam_vec,displ_first_guess)
484 :
485 6146877 : gamma=1.0_r8
486 6146877 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
487 6215254 : else if (tr2.and.tl1.and..not.tl2) then
488 6055204 : 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 6055204 : 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 6055204 : 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 6055204 : num_seg_static,x_start, dgam_vec)
507 6055204 : gamma = 1.0_r8
508 6055204 : gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
509 160050 : else if (tl2.and.tr2) then
510 160050 : 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 160050 : 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 160050 : 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 160050 : num_seg_static,x_start, dgam_vec,displ_first_guess)
532 160050 : gamma = 1.0_r8
533 160050 : 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 69568761552 : do iarea=1,num_area
543 69568761552 : 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 34784380776 : num_seg_max,num_area,dp_area,flowcase,gamma,mass_flux_se(i,j,iside),0.0_r8,gamma_max, &
547 34784380776 : gsweights,gspts,ilev)
548 : !call t_stopf('fvm:swept_area:get_gamma')
549 : !
550 : ! pack segments for high-order weights computation
551 : !
552 69568761552 : do iarea=1,num_area
553 92758348736 : do iseg=1,num_seg_static(iarea)
554 34784380776 : iseg_tmp=num_seg(iarea)+iseg
555 >10435*10^7 : x (:,iseg_tmp,iarea) = x_static (:,iseg,iarea)
556 >16232*10^7 : dx(:,iseg_tmp,iarea) = dx_static(:,iseg,iarea)
557 : end do
558 69568761552 : 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 11594793592 : 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 >13913*10^7 : flux=0.0_r8; flux_tracer=0.0_r8
576 69568761552 : do iarea=1,num_area
577 69568761552 : if (num_seg(iarea)>0) then
578 23612948505 : ii=idx(1,iarea,i,j,iside); jj=idx(2,iarea,i,j,iside)
579 23612948505 : flux=flux+weights(1,iarea)*dp(ii,jj)
580 >28335*10^7 : do itr=1,ntrac
581 >18418*10^8 : do iw=1,irecons_tracer_actual
582 >18181*10^8 : 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 11594793592 : fvm%se_flux(i,j,iside,ilev) = mass_flux_se(i,j,iside)-flux
588 11594793592 : 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 11594793592 : fvm%dp_fvm(i ,j ,ilev ) = fvm%dp_fvm(i ,j ,ilev )-flux
597 >13913*10^7 : 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 11594793592 : if (iside==1) then
602 2941165316 : fvm%dp_fvm(i,j-1,ilev ) = fvm%dp_fvm(i,j-1,ilev )+flux
603 35293983792 : fvm% c(i,j-1,ilev,1:ntrac) = fvm% c(i,j-1,ilev,1:ntrac)+flux_tracer(1:ntrac)
604 : end if
605 11594793592 : if (iside==2) then
606 2851815306 : fvm%dp_fvm(i+1,j,ilev ) = fvm%dp_fvm(i+1,j,ilev )+flux
607 34221783672 : fvm% c(i+1,j,ilev,1:ntrac) = fvm% c(i+1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
608 : end if
609 11594793592 : if (iside==3) then
610 2856231478 : fvm%dp_fvm(i,j+1,ilev ) = fvm%dp_fvm(i,j+1,ilev )+flux
611 34274777736 : fvm% c(i,j+1,ilev,1:ntrac) = fvm% c(i,j+1,ilev,1:ntrac)+flux_tracer(1:ntrac)
612 : end if
613 11594793592 : if (iside==4) then
614 2945581492 : fvm%dp_fvm(i-1,j,ilev ) = fvm%dp_fvm(i-1,j,ilev )+flux
615 35346977904 : 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 483116400 : end subroutine swept_flux
623 :
624 :
625 483116400 : 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 966232800 : real (kind=r8) :: flux,flux_tracer(ntrac)
639 : real (kind=r8), dimension(0:nc+1,0:nc+1) :: inv_dp_area
640 966232800 : real (kind=r8), dimension(0:nc+1,0:nc+1,ntrac):: c_tmp
641 :
642 14976608400 : inv_dp_area=1.0_r8/fvm%dp_fvm(0:nc+1,0:nc+1,ilev)
643 >16522*10^7 : c_tmp = fvm%c(0:nc+1,0:nc+1,ilev,1:ntrac)
644 2415582000 : do iside=1,4
645 9179211600 : do j=jmin_side(iside),jmax_side(iside)
646 31885682400 : do i=imin_side(iside),imax_side(iside)
647 29953216800 : 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 91042920 : do itr=1,ntrac
657 91042920 : flux_tracer(itr) = fvm%se_flux(i,j,iside,ilev)*c_tmp(i,j,itr)*inv_dp_area(i,j)
658 : end do
659 7586910 : fvm%dp_fvm(i ,j ,ilev ) = fvm%dp_fvm(i ,j ,ilev )-flux
660 91042920 : 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 7586910 : if (iside==1) then
665 351172 : fvm%dp_fvm(i,j-1,ilev ) = fvm%dp_fvm(i,j-1,ilev )+flux
666 4214064 : fvm% c(i,j-1,ilev,1:ntrac) = fvm% c(i,j-1,ilev,1:ntrac)+flux_tracer(1:ntrac)
667 : end if
668 7586910 : if (iside==2) then
669 4617836 : fvm%dp_fvm(i+1,j,ilev ) = fvm%dp_fvm(i+1,j,ilev )+flux
670 55414032 : fvm% c(i+1,j,ilev,1:ntrac) = fvm% c(i+1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
671 : end if
672 7586910 : if (iside==3) then
673 591688 : fvm%dp_fvm(i,j+1,ilev ) = fvm%dp_fvm(i,j+1,ilev )+flux
674 7100256 : fvm% c(i,j+1,ilev,1:ntrac) = fvm% c(i,j+1,ilev,1:ntrac)+flux_tracer(1:ntrac)
675 : end if
676 7586910 : if (iside==4) then
677 2026214 : fvm%dp_fvm(i-1,j,ilev ) = fvm%dp_fvm(i-1,j,ilev )+flux
678 24314568 : 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 483116400 : end subroutine large_courant_number_increment
685 :
686 483116400 : 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 483116400 : if (fvm%cubeboundary.NE.0) then
697 373610016 : do j=1-nhe,nc+nhe
698 1930318416 : do i=1-nhe,nc+nhe
699 1556708400 : ishft = NINT(fvm%flux_orient(2,i,j))
700 8094883680 : 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 62268336 : if (fvm%cubeboundary==nwest) then
707 6978348 : var(1-nhe:0,nc+1 :nc+nhe,:) = 0.0_r8
708 61731540 : else if (fvm%cubeboundary==swest) then
709 6978348 : var(1-nhe:0,1-nhe:0 ,:) = 0.0_r8
710 61194744 : else if (fvm%cubeboundary==neast) then
711 6978348 : var(nc+1 :nc+nhe,nc+1 :nc+nhe,:) = 0.0_r8
712 60657948 : else if (fvm%cubeboundary==seast) then
713 6978348 : var(nc+1 :nc+nhe,1-nhe:0,:) = 0.0_r8
714 : end if
715 : end if
716 483116400 : end subroutine ghost_flux_unpack
717 :
718 483116400 : 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 483116400 : num_seg_static(1,1) = 1; num_seg(1,1) = 1; flowcase(1) = -1
760 483116400 : num_seg_static(1,2) = 0; num_seg(1,2) = 2; flowcase(2) = -2
761 483116400 : num_seg_static(1,3) = 1; num_seg(1,3) = 1; flowcase(3) = -1
762 483116400 : num_seg_static(1,4) = 0; num_seg(1,4) = 2; flowcase(4) = -4
763 :
764 1932465600 : do j=1,nc
765 6280513200 : do i=1,nc
766 14493492000 : do ix=1,2
767 8696095200 : iside=1;
768 8696095200 : x_static (ix,1,1,iside,i,j) = fvm%vtx_cart(2,ix,i,j)
769 8696095200 : dx_static(ix,1,1,iside,i,j) = fvm%vtx_cart(1,ix,i,j)-fvm%vtx_cart(2,ix,i,j)
770 8696095200 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(1,ix,i,j)
771 8696095200 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(2,ix,i,j)
772 8696095200 : 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 8696095200 : gamma(iside) = 0.5_r8
777 8696095200 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
778 8696095200 : dx (ix,1,1,iside,i,j) = -dx_static(ix,1,1,iside,i,j)
779 : !
780 : ! side 2
781 : !
782 8696095200 : iside=2;
783 8696095200 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(2,ix,i,j)
784 8696095200 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(3,ix,i,j)
785 8696095200 : dgam_vec (ix,1, iside,i,j) = fvm%vtx_cart(1,ix,i,j)-fvm%vtx_cart(2,ix,i,j)
786 8696095200 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)
787 : !
788 : ! compute first guess - gamma=1
789 : !
790 8696095200 : gamma(iside) = 0.5_r8
791 8696095200 : dx (ix,1,1,iside,i,j) = gamma(iside)*dgam_vec (ix,1, iside,i,j)
792 8696095200 : x (ix,2,1,iside,i,j) = x_start(ix,2,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
793 8696095200 : dx (ix,2,1,iside,i,j) = -gamma(iside)*dgam_vec (ix,1, iside,i,j)
794 : !
795 : ! side 3
796 : !
797 8696095200 : iside=3;
798 8696095200 : x_static (ix,1,1,iside,i,j) = fvm%vtx_cart(4,ix,i,j)
799 8696095200 : dx_static(ix,1,1,iside,i,j) = fvm%vtx_cart(3,ix,i,j)-fvm%vtx_cart(4,ix,i,j)
800 8696095200 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(3,ix,i,j)
801 8696095200 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(4,ix,i,j)
802 8696095200 : 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 8696095200 : gamma(iside) = 0.5_r8
807 8696095200 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
808 8696095200 : dx (ix,1,1,iside,i,j) = -dx_static(ix,1,1,iside,i,j)
809 : !
810 : ! side 4
811 : !
812 8696095200 : iside=4;
813 8696095200 : x_start (ix,1, iside,i,j) = fvm%vtx_cart(1,ix,i,j)
814 8696095200 : x_start (ix,2, iside,i,j) = fvm%vtx_cart(4,ix,i,j)
815 8696095200 : dgam_vec (ix,1, iside,i,j) = fvm%vtx_cart(2,ix,i,j)-fvm%vtx_cart(1,ix,i,j)
816 8696095200 : x (ix,2,1,iside,i,j) = x_start(ix,2,iside,i,j)
817 : !
818 : ! compute first guess - gamma(iside)=1
819 : !
820 8696095200 : gamma(iside) = 0.5_r8
821 8696095200 : dx (ix,2,1,iside,i,j) = gamma(iside)*dgam_vec (ix,1, iside,i,j)
822 8696095200 : x (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
823 13044142800 : 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 1932465600 : do j=1,nc
830 6280513200 : do i=1,nc
831 8696095200 : dp_area = cair(i,j)
832 23189587200 : do iside=1,4
833 17392190400 : flux_se = -fvm%se_flux(i,j,iside,k)
834 21740238000 : if (flux_se>eps) then
835 8696095193 : gamma(iside)=0.5_r8
836 : !
837 : ! this copying is necessary since get_flux_segments_area_iterate change x and dx
838 : !
839 56527519252 : x_tmp (:,1:num_seg(1,iside),:)=x (:,1:num_seg(1,iside),:,iside,i,j)
840 56527519252 : 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 8696095193 : gsweights,gspts,k)
846 26088285579 : 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 8696095193 : 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 8696095207 : 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 483116400 : end subroutine compute_displacements_for_swept_areas
868 :
869 :
870 :
871 20290888785 : subroutine get_flux_segments_area_iterate(x,x_static,dx_static,dx,x_start,dgam_vec,num_seg,num_seg_static,&
872 20290888785 : 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 40581777570 : 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 20290888785 : 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 20290888785 : lexit_after_one_more_iteration = .false.
900 : !
901 : ! compute static line-integrals (not necessary to recompute them for every iteration)
902 : !
903 20290888785 : flux_static = 0.0_r8
904 86960951938 : w_static = 0.0_r8
905 86960951938 : weight_area = 0.0_r8
906 86960951938 : do iarea=1,num_area
907 >10580*10^7 : do iseg=1,num_seg_static(iarea)
908 :
909 : !rck vector directive needed here
910 : !DIR$ SIMD
911 >15652*10^7 : do ipt=1,ngpc
912 >11739*10^7 : xq(ipt) = x_static(1,iseg,iarea)+dx_static(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
913 >11739*10^7 : yq(ipt) = x_static(2,iseg,iarea)+dx_static(2,iseg,iarea)*gspts(ipt)
914 >15652*10^7 : 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 >22319*10^7 : weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx_static(1,iseg,iarea) !integral
917 : end do
918 66670063153 : w_static(iarea)= weight_area(iarea)
919 86960951938 : flux_static = flux_static+weight_area(iarea)*c(iarea) !add to swept flux
920 : end do
921 : !
922 : ! initilization
923 : !
924 20290888785 : gamma1=0.0_r8; f1=-flux ! zero flux guess 1
925 : !
926 : ! compute flux integrals of first guess passed to subroutine
927 : !
928 20290888785 : gamma2=gamma
929 20290888785 : f2 = 0.0_r8
930 86960951938 : weight_area=w_static
931 86960951938 : do iarea=1,num_area
932 >12757*10^7 : do iseg=1,num_seg(iarea)
933 : !rck vector directive needed here
934 : !DIR$ SIMD
935 >24362*10^7 : do ipt=1,ngpc
936 >18271*10^7 : xq(ipt) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
937 >18271*10^7 : yq(ipt) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(ipt)
938 >18271*10^7 : xq2 = xq(ipt)*xq(ipt)
939 >18271*10^7 : xq2i = 1.0_r8/(1.0_r8+xq2)
940 >18271*10^7 : rho = SQRT(1.0_r8+xq2+yq(ipt)*yq(ipt))
941 >18271*10^7 : rhoi = 1.0_r8/rho
942 >18271*10^7 : yrh = yq(ipt)*rhoi
943 >24362*10^7 : F(ipt,1) = yrh*xq2i
944 : enddo
945 >31029*10^7 : weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx(1,iseg,iarea)! integral
946 : end do
947 86960951938 : f2 = f2+weight_area(iarea)*c(iarea)
948 : end do
949 20290888785 : f2 = f2-flux !integral error
950 20290888785 : iter=0
951 20290888785 : if (abs(f2-f1)<eps) then
952 : !
953 : ! in case the first guess is converged
954 : !
955 : return
956 : end if
957 :
958 :
959 20290888785 : dgamma=(gamma2-gamma1)*f2/(f2-f1);
960 20290888785 : gamma3 = gamma2-dgamma; ! Newton "guess" for gamma
961 20290888785 : gamma1 = gamma2; f1 = f2; gamma2 = gamma3; ! prepare for iteration
962 57409129081 : do iter=1,iter_max
963 : !
964 : ! update vertex location: flow_case dependent to avoid many zero operations
965 : !
966 64391737382 : select case(flow_case)
967 : case(-4)
968 6982608301 : iarea=1
969 20947824903 : dx (:,2,1) = gamma3*dgam_vec (:,1)
970 20947824903 : x (:,1,1) = x_start(:,1)+gamma3*dgam_vec(:,1)
971 20947824903 : dx (:,1,1) = -gamma3*dgam_vec (:,1)
972 :
973 : case(-2)
974 6660415100 : iarea=1
975 19981245300 : dx (:,1,iarea) = gamma3*dgam_vec (:,1)
976 19981245300 : x (:,2,iarea) = x_start(:,2)+gamma3*dgam_vec(:,1)
977 19981245300 : dx (:,2,iarea) = -gamma3*dgam_vec (:,1)
978 : case(-1)
979 : !
980 : ! to compute first-guess perpendicular displacements for iside=1
981 : !
982 13445886496 : iarea=1
983 40337659488 : x (:,1,iarea) = x_start(:,1)+gamma3*dgam_vec(:,1)
984 40337659488 : dx (:,1,iarea) = -dx_static(:,1,iarea)
985 40337659488 : x (:,2,iarea) = x_start(:,2)+gamma3*dgam_vec(:,1)
986 40337659488 : dx (:,2,iarea) = x_start(:,2)-x(:,2,iarea)
987 : case(0)
988 589505452 : iarea=3
989 1768516356 : xtmp = x_start(:,1)+gamma3*dgam_vec(:,1)
990 1768516356 : dx (:,1,iarea) = xtmp(: )-x(:,1,iarea) !dynamic - line 2
991 1768516356 : x (:,2,iarea) = xtmp(: ) !dynamic - line 3
992 1768516356 : dx (:,2,iarea) = x_static(:,2,iarea)-x(:,2,iarea) !dynamic - line 3
993 : case(1)
994 371327807 : iarea=2
995 1113983421 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
996 1113983421 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
997 1113983421 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
998 1113983421 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
999 :
1000 371327807 : iarea = 3
1001 1113983421 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1002 1113983421 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1003 1113983421 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1004 1113983421 : x (:,2,iarea) = xtmp (:) !dynamic
1005 1113983421 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1006 1113983421 : x (:,3,iarea) = xtmp2(:) !dynamic
1007 1113983421 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1008 :
1009 371327807 : iarea = 4
1010 1113983421 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1011 1113983421 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1012 1113983421 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1013 1113983421 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1014 : case(2)
1015 13977182686 : iarea=2
1016 41931548058 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1017 41931548058 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1018 41931548058 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
1019 41931548058 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
1020 :
1021 13977182686 : iarea=3
1022 41931548058 : xtmp(: ) = x_start(:,4)+gamma3*dgam_vec(:,4)!
1023 41931548058 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1024 41931548058 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1025 41931548058 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1026 : case(3)
1027 13952940944 : iarea = 3
1028 41858822832 : xtmp (: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1029 41858822832 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1030 41858822832 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1031 41858822832 : dx (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
1032 :
1033 13952940944 : iarea = 4
1034 41858822832 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1035 41858822832 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1036 41858822832 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1037 41858822832 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1038 : case(4)
1039 699563641 : iarea = 1
1040 2098690923 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1041 2098690923 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1042 2098690923 : x (:,2,iarea) = xtmp(:) !dynamic
1043 2098690923 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1044 :
1045 699563641 : iarea = 2
1046 2098690923 : xtmp (: ) = x_start(:,2)+gamma3*dgam_vec(:,2)
1047 2098690923 : xtmp2 (: ) = x_start(:,3)+gamma3*dgam_vec(:,3)
1048 :
1049 2098690923 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1050 :
1051 2098690923 : x (:,2,iarea) = xtmp (:) !dynamic
1052 2098690923 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1053 :
1054 2098690923 : x (:,3,iarea) = xtmp2(:) !dynamic
1055 2098690923 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1056 :
1057 699563641 : iarea = 3
1058 2098690923 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1059 2098690923 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1060 2098690923 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1061 2098690923 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1062 : case(5)
1063 704443260 : iarea = 3
1064 2113329780 : xtmp (: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1065 2113329780 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1066 2113329780 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1067 2113329780 : dx (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
1068 :
1069 704443260 : iarea = 4
1070 2113329780 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1071 2113329780 : xtmp2 (: ) = x_start(:,7)+gamma3*dgam_vec(:,7)
1072 :
1073 2113329780 : dx(:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1074 2113329780 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1075 2113329780 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic - line 2
1076 2113329780 : x (:,3,iarea) = xtmp2(:) !dynamic -line 1
1077 2113329780 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic - line 1
1078 :
1079 704443260 : iarea = 5
1080 2113329780 : xtmp (: ) = x_start(:,8)+gamma3*dgam_vec(:,8)
1081 :
1082 2113329780 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1083 2113329780 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1084 2113329780 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1085 : case(6)
1086 12571833 : iarea = 1
1087 37715499 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1088 37715499 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1089 37715499 : x (:,2,iarea) = xtmp(:) !dynamic
1090 37715499 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1091 :
1092 12571833 : iarea = 2
1093 37715499 : xtmp (: ) = x_start(:,2)+gamma3*dgam_vec(:,2)
1094 37715499 : xtmp2 (: ) = x_start(:,3)+gamma3*dgam_vec(:,3)
1095 :
1096 37715499 : dx(:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1097 37715499 : x (:,2,iarea) = xtmp (:) !dynamic
1098 37715499 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1099 37715499 : x (:,3,iarea) = xtmp2(:) !dynamic
1100 37715499 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1101 :
1102 12571833 : iarea = 3
1103 37715499 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1104 37715499 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1105 37715499 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1106 37715499 : x (:,2,iarea) = xtmp (:) !dynamic
1107 37715499 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1108 37715499 : x (:,3,iarea) = xtmp2(:) !dynamic
1109 37715499 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1110 :
1111 12571833 : iarea = 4
1112 37715499 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1113 37715499 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1114 37715499 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1115 37715499 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1116 : case(7)
1117 12357202 : iarea=2
1118 37071606 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1119 37071606 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1120 37071606 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
1121 37071606 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
1122 :
1123 12357202 : iarea = 3
1124 37071606 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1125 37071606 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1126 37071606 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1127 37071606 : x (:,2,iarea) = xtmp (:) !dynamic
1128 37071606 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1129 37071606 : x (:,3,iarea) = xtmp2(:) !dynamic
1130 37071606 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1131 :
1132 12357202 : iarea = 4
1133 37071606 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1134 37071606 : xtmp2 (: ) = x_start(:,7)+gamma3*dgam_vec(:,7)
1135 :
1136 37071606 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1137 37071606 : x (:,2,iarea) = xtmp(:) !dynamic
1138 37071606 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1139 37071606 : x (:,3,iarea) = xtmp2(:) !dynamic
1140 37071606 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1141 :
1142 12357202 : iarea = 5
1143 37071606 : xtmp (: ) = x_start(:,8)+gamma3*dgam_vec(:,8)
1144 37071606 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1145 37071606 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1146 37071606 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1147 : case(8)
1148 326359 : iarea = 1
1149 979077 : xtmp(: ) = x_start(:,1)+gamma3*dgam_vec(:,1)
1150 979077 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1151 979077 : x (:,2,iarea) = xtmp(:) !dynamic
1152 979077 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1153 :
1154 326359 : iarea = 2
1155 979077 : xtmp (: ) = x_start(:,2)+gamma3*dgam_vec(:,2)
1156 979077 : xtmp2 (: ) = x_start(:,3)+gamma3*dgam_vec(:,3)
1157 :
1158 979077 : dx(:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1159 979077 : x (:,2,iarea) = xtmp (:) !dynamic
1160 979077 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1161 979077 : x (:,3,iarea) = xtmp2(:) !dynamic
1162 979077 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1163 :
1164 326359 : iarea = 3
1165 979077 : xtmp (: ) = x_start(:,4)+gamma3*dgam_vec(:,4)
1166 979077 : xtmp2(: ) = x_start(:,5)+gamma3*dgam_vec(:,5)
1167 979077 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1168 979077 : x (:,2,iarea) = xtmp (:) !dynamic
1169 979077 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1170 979077 : x (:,3,iarea) = xtmp2(:) !dynamic
1171 979077 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1172 :
1173 326359 : iarea = 4
1174 979077 : xtmp (: ) = x_start(:,6)+gamma3*dgam_vec(:,6)
1175 979077 : xtmp2 (: ) = x_start(:,7)+gamma3*dgam_vec(:,7)
1176 :
1177 979077 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1178 979077 : x (:,2,iarea) = xtmp(:) !dynamic
1179 979077 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1180 979077 : x (:,3,iarea) = xtmp2(:) !dynamic
1181 979077 : dx (:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1182 :
1183 326359 : iarea = 5
1184 979077 : xtmp (: ) = x_start(:,8)+gamma3*dgam_vec(:,8)
1185 979077 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1186 979077 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1187 979077 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1188 : case default
1189 57409129081 : call endrun('flow case not defined in get_flux_segments_area_iterate')
1190 : end select
1191 : !
1192 : ! compute flux integral
1193 : !
1194 >23609*10^7 : f2 = 0.0_r8
1195 >23609*10^7 : weight_area=w_static
1196 >23609*10^7 : do iarea=1,num_area
1197 >34500*10^7 : do iseg=1,num_seg(iarea)
1198 : !rck vector directive needed here
1199 : !DIR$ SIMD
1200 >66524*10^7 : do ipt=1,ngpc
1201 :
1202 >49893*10^7 : xq(ipt) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
1203 >49893*10^7 : yq(ipt) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(ipt)
1204 :
1205 >49893*10^7 : xq2 = xq(ipt)*xq(ipt)
1206 >49893*10^7 : xq2i = 1.0_r8/(1.0_r8+xq2)
1207 >49893*10^7 : rho = SQRT(1.0_r8+xq2+yq(ipt)*yq(ipt))
1208 >49893*10^7 : rhoi = 1.0_r8/rho
1209 >49893*10^7 : yrh = yq(ipt)*rhoi
1210 >66524*10^7 : F(ipt,1) = yrh*xq2i
1211 : end do
1212 >84393*10^7 : weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx(1,iseg,iarea)! integral
1213 : end do
1214 >23609*10^7 : f2 = f2+weight_area(iarea)*c(iarea)
1215 : end do
1216 57409129081 : 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 77700017866 : if (ABS(f2)<eps.or.lexit_after_one_more_iteration) then
1224 20303958125 : gamma=gamma3
1225 20303958125 : if (gamma>gamma_max) then
1226 13069340 : lexit_after_one_more_iteration=.true.
1227 13069340 : gamma=gamma_max
1228 13069340 : gamma3=gamma_max
1229 : else
1230 : exit
1231 : end if
1232 : else
1233 : !
1234 : ! Newton increment
1235 : !
1236 37105170956 : 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 79 : dgamma=-0.5_r8*(gamma2-gamma1)
1241 79 : lexit_after_one_more_iteration=.true.
1242 : else
1243 37105170877 : dgamma=(gamma2-gamma1)*f2/(f2-f1)
1244 : endif
1245 37105170956 : if (ABS(dgamma)>eps) then
1246 37105170956 : 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 37105170956 : gamma3=MAX(gamma3,gamma_min)
1254 : !
1255 : ! prepare for next iteration
1256 : !
1257 37105170956 : gamma1 = gamma2; f1 = f2; gamma2 = gamma3;
1258 : endif
1259 : end do
1260 20290888785 : 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 483116400 : 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 483116400 : ib = fvm%cubeboundary
1314 12077910000 : 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 12561026400 : 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 10145444400 : 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 483116400 : if (ib>0) then
1338 62268336 : if (ib==swest) degenerate(1 ,1 ) = 1
1339 61731540 : if (ib==nwest) degenerate(1 ,nc+1) = 1
1340 61731540 : if (ib==neast) degenerate(nc+1,nc+1) = 1
1341 61731540 : if (ib==seast) degenerate(nc+1,1 ) = 1
1342 : end if
1343 :
1344 2415582000 : do j=1,nc+1
1345 10145444400 : do i=1,nc+1
1346 17392190400 : do sgn=-1,1,2
1347 : if (&
1348 61838899200 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1349 85028486400 : sgn*flux_sum(i ,j,1)>0.0_r8.and.sgn*flux_sum(i,j ,2)<0.0_r8) then
1350 12082326 : circular_flow(i,j) = 0
1351 : else
1352 15447642474 : circular_flow(i,j) = 1
1353 : end if
1354 : end do
1355 : end do
1356 : end do
1357 : !
1358 : ! wrap around corners
1359 : !
1360 483116400 : if (ib==nwest) then
1361 536796 : flux_sum(0,nc+1,1) = fvm%se_flux(0,nc,3,ilev)-fvm%se_flux(1,nc+1,4,ilev)
1362 536796 : flux_sum(1,nc+1,2) = fvm%se_flux(0,nc,3,ilev)-fvm%se_flux(1,nc+1,4,ilev)
1363 :
1364 536796 : i=1;j=nc+1;
1365 536796 : circular_flow(i,j) = 1
1366 1610388 : do sgn=-1,1,2
1367 : if (&
1368 : sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1369 1610388 : sgn*flux_sum(i ,j,1)>0.0_r8.and.sgn*flux_sum(i,j ,2)<0.0_r8) then
1370 2097 : circular_flow(i,j) = 0
1371 : end if
1372 : end do
1373 482579604 : else if (ib==swest) then
1374 536796 : flux_sum(0,1,1) = fvm%se_flux(1,0,4,ilev)-fvm%se_flux(0,1,1,ilev)
1375 536796 : flux_sum(1,0,2) = fvm%se_flux(0,1,1,ilev)-fvm%se_flux(1,0,4,ilev)
1376 536796 : i=1;j=1;
1377 536796 : circular_flow(i,j) = 1
1378 1610388 : do sgn=-1,1,2
1379 : if (&
1380 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.&
1381 1610388 : sgn*flux_sum(i ,j,1)>0.0_r8.and.sgn*flux_sum(i,j ,2)<0.0_r8) then
1382 2227 : circular_flow(i,j) = 0
1383 : end if
1384 : end do
1385 482042808 : else if (ib==neast) then
1386 536796 : flux_sum(nc+1,nc+1,1) = fvm%se_flux(nc+1,nc,3,ilev)-fvm%se_flux(nc,nc+1,2,ilev)
1387 536796 : flux_sum(nc+1,nc+1,2) = fvm%se_flux(nc,nc+1,2,ilev)-fvm%se_flux(nc+1,nc,3,ilev)
1388 536796 : i=nc+1;j=nc+1;
1389 536796 : circular_flow(i,j) = 1
1390 1610388 : do sgn=-1,1,2
1391 : if (&
1392 1073592 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1393 536796 : sgn*flux_sum(i,j ,2)<0.0_r8) then
1394 828 : circular_flow(i,j) = 0
1395 : end if
1396 : end do
1397 481506012 : else if (ib==seast) then
1398 536796 : flux_sum(nc+1,1 ,1) = fvm%se_flux(nc,0,2,ilev)-fvm%se_flux(nc+1,1,1,ilev)
1399 536796 : flux_sum(nc+1,0 ,2) = fvm%se_flux(nc,0,2,ilev)-fvm%se_flux(nc+1,1,1,ilev)
1400 536796 : i=nc+1;j=1;
1401 536796 : circular_flow(i,j) = 1
1402 1610388 : do sgn=-1,1,2
1403 : if (&
1404 1073592 : sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
1405 536796 : sgn*flux_sum(i,j ,2)<0.0_r8) then
1406 1742 : circular_flow(i,j) = 0
1407 : end if
1408 : end do
1409 : end if
1410 10145444400 : illcond = circular_flow*degenerate
1411 : !
1412 : !
1413 : !
1414 : !
1415 2415582000 : do iside=1,4
1416 9179211600 : do j=jmin_side(iside),jmax_side(iside)
1417 31885682400 : do i=imin_side(iside),imax_side(iside)
1418 29953216800 : if (fvm%se_flux(i,j,iside,ilev)>eps) then
1419 11594793592 : iur = i+idx_shift(4,iside); jur = j+idy_shift(4,iside) !(i,j) index of upper right quadrant
1420 11594793592 : ilr = i+idx_shift(5,iside); jlr = j+idy_shift(5,iside) !(i,j) index of lower left quadrant
1421 11594793592 : iul = i+idx_shift(2,iside); jul = j+idy_shift(2,iside) !(i,j) index of upper right quadrant
1422 11594793592 : ill = i+idx_shift(1,iside); jll = j+idy_shift(1,iside) !(i,j) index of lower left quadrant
1423 :
1424 : !iside=1
1425 11594793592 : if (iside==1) then
1426 2941165316 : displ(0,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i,j) !center left
1427 2941165316 : displ(1,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i+1,j) !center right
1428 2941165316 : displ(2,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j) !c2
1429 2941165316 : displ(3,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j) !c3
1430 2941165316 : displ(4,i,j,iside) = -flux_sum (i+1,j ,1)*illcond(i+1,j) !r1
1431 2941165316 : displ(5,i,j,iside) = -flux_sum (i+1,j-1,2)*illcond(i+1,j) !r2
1432 2941165316 : displ(6,i,j,iside) = -flux_sum (i-1,j ,1)*illcond(i ,j) !l1
1433 2941165316 : displ(7,i,j,iside) = flux_sum (i ,j-1,2)*illcond(i ,j) !l2
1434 :
1435 : end if
1436 11594793592 : if (iside==2) then
1437 : !iside=2
1438 2851815306 : displ(0,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j ) !center left
1439 2851815306 : displ(1,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j+1) !center right
1440 2851815306 : displ(2,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i+1,j+1) !c2
1441 2851815306 : displ(3,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i+1,j ) !c3
1442 2851815306 : displ(4,i,j,iside) = flux_sum (i+1,j+1,2)*illcond(i+1,j+1) !r1
1443 2851815306 : displ(5,i,j,iside) = -flux_sum (i+1,j+1,1)*illcond(i+1,j+1) !r2
1444 2851815306 : displ(6,i,j,iside) = flux_sum (i+1,j-1,2)*illcond(i+1,j) !l1
1445 2851815306 : displ(7,i,j,iside) = flux_sum (i+1,j ,1)*illcond(i+1,j) !l2
1446 : end if
1447 : !iside=3
1448 11594793592 : if (iside==3) then
1449 2856231478 : displ(0,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i+1,j+1) !center left
1450 2856231478 : displ(1,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i ,j+1) !center right
1451 2856231478 : displ(2,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j+1) !c2
1452 2856231478 : displ(3,i,j,iside) = flux_sum (i+1,j ,2)*illcond(i+1,j+1) !c3
1453 2856231478 : displ(4,i,j,iside) = flux_sum (i-1,j+1,1)*illcond(i ,j+1) !r1
1454 2856231478 : displ(5,i,j,iside) = flux_sum (i ,j+1,2)*illcond(i ,j+1) !r2
1455 2856231478 : displ(6,i,j,iside) = flux_sum (i+1,j+1,1)*illcond(i+1,j+1) !l1
1456 2856231478 : displ(7,i,j,iside) = -flux_sum (i+1,j+1,2)*illcond(i+1,j+1) !l2
1457 : end if
1458 11594793592 : if (iside==4) then
1459 : !iside=4
1460 2945581492 : displ(0,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j+1) !center left
1461 2945581492 : displ(1,i,j,iside) = -flux_sum (i ,j ,2)*illcond(i ,j ) !center right
1462 2945581492 : displ(2,i,j,iside) = -flux_sum (i ,j ,1)*illcond(i ,j ) !c2
1463 2945581492 : displ(3,i,j,iside) = flux_sum (i ,j+1,1)*illcond(i ,j+1) !c3
1464 2945581492 : displ(4,i,j,iside) = -flux_sum (i ,j-1,2)*illcond(i ,j ) !r1
1465 2945581492 : displ(5,i,j,iside) = flux_sum (i-1,j ,1)*illcond(i ,j ) !r2
1466 2945581492 : displ(6,i,j,iside) = -flux_sum (i ,j+1,2)*illcond(i ,j+1) !l1
1467 2945581492 : displ(7,i,j,iside) = -flux_sum (i-1,j+1,1)*illcond(i ,j+1) !l2
1468 : end if
1469 :
1470 34784380776 : base_vtx(:,1,i,j,iside) = fvm%vtx_cart(iside,:,i ,j ) !vertex center left
1471 34784380776 : base_vtx(:,2,i,j,iside) = fvm%vtx_cart(iside_p1(iside),:,i ,j ) !vertex center right
1472 34784380776 : base_vtx(:,3,i,j,iside) = fvm%vtx_cart(iside,:,iur,jur ) !vertex upper right
1473 34784380776 : base_vtx(:,4,i,j,iside) = fvm%vtx_cart(iside_p3(iside),:,ilr,jlr) !vertex lower right
1474 34784380776 : base_vtx(:,5,i,j,iside) = fvm%vtx_cart(iside_p1(iside),:,iul,jul) !vertex upper left
1475 34784380776 : base_vtx(:,6,i,j,iside) = fvm%vtx_cart(iside_p2(iside),:,ill,jll) !vertex lower left
1476 :
1477 34784380776 : base_vec(:, 1,i,j,iside) = fvm%flux_vec (:,i ,j ,iside ) !vector center
1478 34784380776 : base_vec(:, 2,i,j,iside) = fvm%flux_vec (:,i ,j ,iside_p1(iside)) !vector center right
1479 34784380776 : base_vec(:, 3,i,j,iside) = fvm%flux_vec (:,i ,j ,iside_p3(iside)) !vector center left
1480 34784380776 : base_vec(:, 4,i,j,iside) = fvm%flux_vec (:,iur,jur,iside ) !vector upper right 1
1481 34784380776 : base_vec(:, 5,i,j,iside) = fvm%flux_vec (:,iur,jur,iside_p3(iside)) !vector upper right 2
1482 34784380776 : base_vec(:, 6,i,j,iside) = fvm%flux_vec (:,ilr,jlr,iside_p3(iside)) !vector lower right 1
1483 34784380776 : base_vec(:, 7,i,j,iside) = fvm%flux_vec (:,ilr,jlr,iside_p2(iside)) !vector lower right 2
1484 34784380776 : base_vec(:, 8,i,j,iside) = fvm%flux_vec (:,iul,jul,iside ) !vector upper left 1
1485 34784380776 : base_vec(:, 9,i,j,iside) = fvm%flux_vec (:,iul,jul,iside_p1(iside)) !vector upper left 2
1486 34784380776 : base_vec(:,10,i,j,iside) = fvm%flux_vec (:,ill,jll,iside_p1(iside)) !vector lower left 1
1487 34784380776 : base_vec(:,11,i,j,iside) = fvm%flux_vec (:,ill,jll,iside_p2(iside)) !vector lower left 2
1488 :
1489 69568761552 : do iarea=1,5
1490 57973967960 : idx(1,iarea,i,j,iside) = i+idx_shift(iarea,iside)
1491 69568761552 : idx(2,iarea,i,j,iside) = j+idy_shift(iarea,iside)
1492 : end do
1493 : else
1494 >10435*10^7 : 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 483116400 : 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 224034279 : 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 672102837 : if (SUM(ABS(base_vec(:,9,i,j,iside))).NE.0) then
1555 223574651 : 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 223574651 : 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 459628 : gamma=displ(0,i,j,iside)
1563 : end if
1564 :
1565 :
1566 672102837 : 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 672102837 : x_start (:,1) = base_vtx(:, 6,i,j,iside)
1568 672102837 : dgam_vec(:,1) = base_vec(:,10,i,j,iside)*gamma
1569 :
1570 672102837 : 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 224034279 : iarea = 1
1573 224034279 : num_seg (iarea) = 2
1574 224034279 : num_seg_static(iarea) = 1
1575 :
1576 672102837 : x_static (:,1,iarea) = base_vtx(:,6,i,j,iside) !static
1577 672102837 : dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static
1578 :
1579 672102837 : xtmp(: ) = x_start(:,1)+dgam_vec(:,1)
1580 672102837 : x (:,1,iarea) = xdep(:,1) !static
1581 672102837 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1582 :
1583 672102837 : x (:,2,iarea) = xtmp(:) !dynamic
1584 672102837 : dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic
1585 : !
1586 : !
1587 : !
1588 224034279 : iarea = 2
1589 224034279 : num_seg (iarea) = 3
1590 :
1591 672102837 : x_start (:,2) = base_vtx(:,5,i,j,iside)
1592 672102837 : dgam_vec(:,2) = base_vec(:,9,i,j,iside)*gamma
1593 672102837 : xtmp (: ) = x_start(:,2)+dgam_vec(:,2)
1594 :
1595 672102837 : x_start (:,3) = base_vtx(:,5,i,j,iside)
1596 672102837 : dgam_vec(:,3) = base_vec(:,8,i,j,iside)*displ(0,i,j,iside)
1597 672102837 : xtmp2 (: ) = x_start(:,3)+dgam_vec(:,3)
1598 :
1599 672102837 : x (:,1,iarea) = base_vtx(:,5,i,j,iside) !static
1600 672102837 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1601 :
1602 672102837 : x (:,2,iarea) = xtmp (:) !dynamic
1603 672102837 : dx(:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1604 :
1605 672102837 : x (:,3,iarea) = xtmp2(:) !dynamic
1606 672102837 : dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:) !dynamic
1607 224034279 : end subroutine define_area1_area2
1608 :
1609 :
1610 5560957800 : 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 16682873400 : 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 16682873400 : x_start (:,1) = base_vtx(:,5,i,j,iside)
1647 5560957800 : gamma = displ(0,i,j,iside)
1648 16682873400 : dgam_vec(:,1) = base_vec(:,8,i,j,iside)*gamma
1649 5560957800 : if (present(displ_first_guess)) displ_first_guess = gamma
1650 :
1651 5560957800 : iarea = 2
1652 5560957800 : num_seg (iarea) = 2
1653 5560957800 : num_seg_static(iarea) = 1
1654 :
1655 16682873400 : x_static (:,1,iarea) = base_vtx(:,5,i,j,iside) !static - line 1
1656 16682873400 : dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static - line 1
1657 :
1658 16682873400 : xtmp (: ) = x_start(:,1)+dgam_vec(:,1)
1659 16682873400 : x (:,1,iarea) = xdep(:,1) !static - line 2
1660 16682873400 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1661 :
1662 16682873400 : x (:,2,iarea) = xtmp(:) !dynamic - line 3
1663 16682873400 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
1664 5560957800 : end subroutine define_area2
1665 :
1666 :
1667 5598739471 : 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 5598739471 : xdep(:,2) = base_vtx(:,2,i,j,iside)+displ(1,i,j,iside)*base_vec(:,1,i,j,iside)&
1700 22394957884 : +MAX(0.0_r8,displ(2,i,j,iside))*base_vec(:,2,i,j,iside)
1701 16796218413 : x_start (:,4) = base_vtx(:,1,i,j,iside)
1702 5598739471 : gamma = displ(0,i,j,iside)
1703 16796218413 : dgam_vec(:,4) = base_vec(:,1,i,j,iside)*gamma
1704 16796218413 : xtmp (: ) = x_start(:,4)+dgam_vec(:,4)
1705 :
1706 5598739471 : if (present(displ_first_guess)) displ_first_guess = gamma
1707 :
1708 5598739471 : iarea = 3
1709 5598739471 : num_seg (iarea) = 2
1710 5598739471 : num_seg_static(iarea) = 2
1711 :
1712 16796218413 : x_static (:,1,iarea) = xdep(:,2) !static - line 3
1713 16796218413 : dx_static(:,1,iarea) = base_vtx(:,2,i,j,iside)-xdep(:,2) !static - line 3
1714 :
1715 16796218413 : x_static (:,2,iarea) = base_vtx(:,2,i,j,iside) !static - line 4
1716 16796218413 : dx_static(:,2,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 4
1717 :
1718 16796218413 : x (:,1,iarea) = base_vtx(:,1,i,j,iside) !static - line 1
1719 16796218413 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1720 :
1721 16796218413 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1722 16796218413 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1723 5598739471 : end subroutine define_area3_left
1724 :
1725 5597624165 : 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 5597624165 : xdep(:,1) = base_vtx(:,1,i,j,iside)+displ(0,i,j,iside)*base_vec(:,1,i,j,iside)&
1755 22390496660 : +MAX(0.0_r8,displ(3,i,j,iside))*base_vec(:,3,i,j,iside)
1756 16792872495 : x_start (:,5) = base_vtx(:,2,i,j,iside)
1757 5597624165 : gamma = displ(1,i,j,iside)
1758 16792872495 : dgam_vec(:,5) = base_vec(:,1,i,j,iside)*gamma
1759 16792872495 : xtmp (: ) = x_start(:,5)+dgam_vec(:,5)
1760 :
1761 5597624165 : iarea = 3
1762 5597624165 : num_seg (iarea) = 2
1763 5597624165 : num_seg_static(iarea) = 2
1764 :
1765 16792872495 : x_static (:,1,iarea) = base_vtx(:,1,i,j,iside) !static - line 1
1766 16792872495 : dx_static(:,1,iarea) = xdep(:,1)-base_vtx(:,1,i,j,iside) !static - line 1
1767 :
1768 16792872495 : x_static (:,2,iarea) = base_vtx(:,2,i,j,iside) !static - line 4
1769 16792872495 : dx_static(:,2,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 4
1770 :
1771 16792872495 : x (:,1,iarea) = xdep(:,1) !static - line 2
1772 16792872495 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1773 :
1774 16792872495 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1775 16792872495 : dx (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
1776 5597624165 : end subroutine define_area3_right
1777 :
1778 :
1779 186252608 : 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 558757824 : x_start (:,4) = base_vtx(:,1,i,j,iside)
1805 558757824 : x_start (:,5) = base_vtx(:,2,i,j,iside)
1806 186252608 : gamma = displ(0,i,j,iside)
1807 558757824 : dgam_vec(:,4) = base_vec(:,1,i,j,iside)*gamma
1808 558757824 : dgam_vec(:,5) = base_vec(:,1,i,j,iside)*gamma
1809 558757824 : xtmp (: ) = x_start(:,4)+dgam_vec(:,4)
1810 558757824 : xtmp2 (: ) = x_start(:,5)+dgam_vec(:,5)
1811 :
1812 186252608 : iarea = 3
1813 186252608 : num_seg (iarea) = 3
1814 186252608 : num_seg_static(iarea) = 1
1815 :
1816 558757824 : x_static (:,1,iarea) = base_vtx(:,2,i,j,iside) !static
1817 558757824 : dx_static(:,1,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static
1818 :
1819 558757824 : x (:,1,iarea) = base_vtx(:,1,i,j,iside) !static
1820 558757824 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
1821 :
1822 558757824 : x (:,2,iarea) = xtmp (:) !dynamic
1823 558757824 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
1824 :
1825 558757824 : x (:,3,iarea) = xtmp2(:) !dynamic
1826 558757824 : dx (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
1827 186252608 : end subroutine define_area3_left_right
1828 :
1829 225251782 : 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 225251782 : iarea = 4
1862 225251782 : num_seg (iarea) = 3
1863 :
1864 675755346 : if (SUM(ABS(base_vec(:,5,i,j,iside))).NE.0) then
1865 224997980 : 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 224997980 : 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 253802 : gamma = displ(1,i,j,iside)
1873 : end if
1874 :
1875 225251782 : if (present(displ_first_guess)) displ_first_guess = displ(1,i,j,iside)
1876 :
1877 675755346 : x_start (:,6) = base_vtx(:,3,i,j,iside)
1878 675755346 : dgam_vec(:,6) = base_vec(:,4,i,j,iside)*displ(1,i,j,iside)
1879 675755346 : xtmp (: ) = x_start(:,6)+dgam_vec(:,6)
1880 675755346 : x_start (:,7) = base_vtx(:,3,i,j,iside)
1881 675755346 : dgam_vec(:,7) = base_vec(:,5,i,j,iside)*gamma
1882 675755346 : xtmp2 (: ) = x_start(:,7)+dgam_vec(:,7)
1883 :
1884 675755346 : x (:,1,iarea) = base_vtx(:,3,i,j,iside)!static -line 1
1885 675755346 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1886 :
1887 675755346 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1888 675755346 : dx (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic - line 2
1889 :
1890 675755346 : x (:,3,iarea) = xtmp2(:) !static -line 1
1891 675755346 : 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 675755346 : -displ(4,i,j,iside)*base_vec(:,7,i,j,iside)
1897 675755346 : x_start (:,8) = base_vtx(:,4,i,j,iside)
1898 675755346 : dgam_vec(:,8) = base_vec(:,6,i,j,iside)*gamma
1899 675755346 : xtmp (: ) = x_start(:,8)+dgam_vec(:,8)
1900 :
1901 225251782 : iarea = 5
1902 225251782 : num_seg (iarea) = 2
1903 225251782 : num_seg_static(iarea) = 1
1904 :
1905 675755346 : x (:,1,iarea) = base_vtx(:,4,i,j,iside)!static -line 1
1906 675755346 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 1
1907 :
1908 675755346 : x_static (:,1,iarea) = xdep(:,1) !static - line 1
1909 675755346 : dx_static(:,1,iarea) = x(:,1,iarea)-x_static(:,1,iarea) !static - line 1
1910 :
1911 675755346 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1912 675755346 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1913 225251782 : end subroutine define_area4_area5
1914 :
1915 :
1916 5558624991 : 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 5558624991 : iarea = 4
1939 5558624991 : num_seg (iarea) = 2
1940 5558624991 : num_seg_static(iarea) = 1
1941 :
1942 5558624991 : xdep(:,1) = base_vtx(:,3,i,j,iside)+MAX(0.0_r8,displ(4,i,j,iside))*base_vec(:,4,i,j,iside)&
1943 22234499964 : -displ(2,i,j,iside)*base_vec(:,5,i,j,iside)
1944 16675874973 : x_start (:,6) = base_vtx(:,3,i,j,iside)
1945 5558624991 : gamma = displ(1,i,j,iside)
1946 16675874973 : dgam_vec(:,6) = base_vec(:,4,i,j,iside)*gamma
1947 16675874973 : xtmp (: ) = x_start(:,6)+dgam_vec(:,6)
1948 :
1949 5558624991 : if (present(displ_first_guess)) displ_first_guess = gamma
1950 :
1951 16675874973 : x_static (:,1,iarea) = xdep(:,1) !static
1952 16675874973 : dx_static(:,1,iarea) = base_vtx(:,3,i,j,iside)-xdep(:,1) !static
1953 :
1954 16675874973 : x (:,1,iarea) = base_vtx(:,3,i,j,iside) !static - line 2
1955 16675874973 : dx (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
1956 :
1957 16675874973 : x (:,2,iarea) = xtmp(:) !dynamic -line 2
1958 16675874973 : dx (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
1959 5558624991 : end subroutine define_area4
1960 :
1961 212177348 : 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 636532044 : 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 636532044 : displ(1,i,j,iside)*base_vec(:,1,i,j,iside)+displ(2,i,j,iside)*base_vec(:,2,i,j,iside)
1998 636532044 : xdep(:,2) = 0.5_r8*(xdep(:,1)+xdep(:,3))
1999 :
2000 212177348 : gamma= se_flux_center
2001 : x_start(:,1) = ABS(base_vec(:,3,i,j,iside))*((xdep(:,2)-base_vtx(:,1,i,j,iside)))+&
2002 636532044 : base_vtx(:,1,i,j,iside) !xdep(2) - midway between departure points projected to side 1
2003 :
2004 636532044 : dgam_vec(:,1) = gamma*base_vec(:,1,i,j,iside)
2005 :
2006 212177348 : if (present(displ_first_guess)) displ_first_guess = gamma
2007 :
2008 636532044 : xdep(:,2) = x_start(:,1)+dgam_vec(:,1)
2009 212177348 : iarea = 3
2010 212177348 : num_seg (iarea) = 2
2011 212177348 : num_seg_static(iarea) = 3
2012 :
2013 : ! ______X______
2014 : ! || 2 / \ 3 ||
2015 : ! || *--/ \--* ||
2016 : ! || / \ ||
2017 : ! ||/ 1 5 4\||
2018 : ! ========================================
2019 : ! || ||
2020 : !
2021 636532044 : x_static (:,1,iarea) = base_vtx(:,1,i,j,iside) !static - line 1
2022 636532044 : dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static - line 1
2023 :
2024 636532044 : x (:,1,iarea) = xdep(:,1) !static - line 2
2025 636532044 : dx (:,1,iarea) = xdep(:,2)-x(:,1,iarea) !dynamic - line 2
2026 :
2027 636532044 : x (:,2,iarea) = xdep(:,2) !dynamic - line 3
2028 636532044 : dx (:,2,iarea) = xdep(:,3)-x(:,2,iarea) !dynamic - line 3
2029 :
2030 636532044 : x_static (:,2,iarea) = xdep(:,3) !static - line 4
2031 636532044 : dx_static(:,2,iarea) = base_vtx(:,2,i,j,iside)-x_static(:,2,iarea)!static - line 4
2032 :
2033 636532044 : x_static (:,3,iarea) = base_vtx(:,2,i,j,iside) !static - line 5
2034 636532044 : dx_static(:,3,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 5
2035 :
2036 212177348 : end subroutine define_area3_center
2037 : end module fvm_consistent_se_cslam
|