LCOV - code coverage report
Current view: | top level - dynamics/se/dycore - fvm_consistent_se_cslam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 848 | 866 | 97.9 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 15 | 15 | 100.0 %

              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

Generated by: LCOV version 1.14