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: | 2024-12-17 22:39:59 |  | 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    31030272 :     llimiter = .true.
          79             :
          80      738816 :     inJetCall = .false.
          81      738816 :     if(((kminp .ne. 1) .or. (kmaxp .ne. nlev)) .and. vert_num_threads>1) then
          82           0 :        write(iulog,*)'WARNING: deactivating vertical threading for JET region call'
          83           0 :        inJetCall = .true.
          84           0 :        region_num_threads = 1
          85             :     else
          86      738816 :        region_num_threads = vert_num_threads
          87             :     endif
          88             :
          89      738816 :     call omp_set_nested(.true.)
          90             :     !$OMP PARALLEL NUM_THREADS(region_num_threads), DEFAULT(SHARED), &
          91             :     !$OMP PRIVATE(hybridnew,kblk,ie,k,kmin,gspts,inv_dp_area,itr), &
          92             :     !$OMP PRIVATE(kmin_jet_local,kmax,kmax_jet_local,kptr,q,ctracer,ActiveJetThread)
          93      738816 :     call gauss_points(ngpc,gsweights,gspts) !set gauss points/weights
          94     2955264 :     gspts = 0.5_r8*(gspts+1.0_r8) !shift location so in [0:1] instead of [-1:1]
          95             :
          96      738816 :     if(inJetCall) then
          97             :       ! ===============================================================================
          98             :       ! if this is the reduced Jet region call then do not thread over the vertical....
          99             :       ! Just use the number of vertical levels that were passed into subroutine
         100             :       ! ===============================================================================
         101           0 :       hybridnew = config_thread_region(hybrid,'serial')
         102           0 :       kmin = kminp
         103           0 :       kmax = kmaxp
         104             :     else
         105      738816 :       hybridnew = config_thread_region(hybrid,'vertical')
         106      738816 :       call get_loop_ranges(hybridnew,kbeg=kmin,kend=kmax)
         107             :     endif
         108             :
         109      738816 :     kblk = kmax-kmin+1
         110             :     !call t_startf('fvm:before_Qnhc')
         111     5933616 :     do ie=nets,nete
         112   488311200 :        do k=kmin,kmax
         113 25605169200 :           elem(ie)%sub_elem_mass_flux(:,:,:,k) = dt_fvm*elem(ie)%sub_elem_mass_flux(:,:,:,k)*fvm(ie)%dp_ref_inverse(k)
         114  6285708000 :           fvm(ie)%dp_fvm(1:nc,1:nc,k)          =         fvm(ie)%dp_fvm (1:nc,1:nc,k)*fvm(ie)%dp_ref_inverse(k)
         115             :        end do
         116     5194800 :        kptr = kmin-1
         117     5194800 :        call ghostpack(ghostbufQnhc,fvm(ie)%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax)   ,kblk,      kptr,ie)
         118   218920416 :        do q=1,ntrac
         119   212986800 :           kptr = kptr + nlev
         120   218181600 :           call ghostpack(ghostbufQnhc,fvm(ie)%c(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax,q),kblk,kptr,ie)
         121             :        enddo
         122             :     end do
         123             :     !call t_stopf('fvm:before_Qnhc')
         124             :     !call t_startf('fvm:ghost_exchange:Qnhc')
         125      738816 :     call ghost_exchange(hybridnew,ghostbufQnhc,location='ghostbufQnhc')
         126             :     !call t_stopf('fvm:ghost_exchange:Qnhc')
         127             :     !call t_startf('fvm:orthogonal_swept_areas')
         128     5933616 :     do ie=nets,nete
         129   488311200 :       do k=kmin,kmax
         130 25610364000 :         fvm(ie)%se_flux    (1:nc,1:nc,:,k) = elem(ie)%sub_elem_mass_flux(:,:,:,k)
         131             :       end do
         132     5194800 :       kptr = kmin-1
         133     5194800 :       call ghostunpack(ghostbufQnhc, fvm(ie)%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax)   , kblk      ,kptr,ie)
         134   218181600 :       do q=1,ntrac
         135   212986800 :          kptr = kptr + nlev
         136   218181600 :          call ghostunpack(ghostbufQnhc, fvm(ie)%c(1-nhc:nc+nhc,1-nhc:nc+nhc,kmin:kmax,q),kblk,kptr,ie)
         137             :       enddo
         138   488311200 :       do k=kmin,kmax
         139   488311200 :         call compute_displacements_for_swept_areas (fvm(ie),fvm(ie)%dp_fvm(:,:,k),k,gsweights,gspts)
         140             :       end do
         141     5194800 :       kptr = 4*(kmin-1)
         142     5933616 :       call ghostpack(ghostBufFlux, fvm(ie)%se_flux(:,:,:,kmin:kmax),4*kblk,kptr,ie)
         143             :     end do
         144             :
         145      738816 :     call ghost_exchange(hybridnew,ghostBufFlux,location='ghostBufFlux')
         146             :
         147     5933616 :     do ie=nets,nete
         148     5194800 :       kptr = 4*(kmin-1)
         149     5194800 :       call ghostunpack(ghostBufFlux, fvm(ie)%se_flux(:,:,:,kmin:kmax),4*kblk,kptr,ie)
         150   489050016 :       do k=kmin,kmax
         151   488311200 :          call ghost_flux_unpack(fvm(ie),fvm(ie)%se_flux(:,:,:,k))
         152             :       end do
         153             :     enddo
         154             :
         155             :     !call t_stopf('fvm:orthogonal_swept_areas')
         156     5933616 :     do ie=nets,nete
         157             :        ! Intel compiler version 2023.0.0 on derecho had significant slowdown on subroutine interface without
         158             :        ! these pointers.
         159     5194800 :       fcube => fvm(ie)%c(:,:,:,:)
         160     5194800 :       spherecentroid => fvm(ie)%spherecentroid(:,1-nhe:nc+nhe,1-nhe:nc+nhe)
         161   489050016 :       do k=kmin,kmax
         162             :          !call t_startf('FVM:tracers_reconstruct')
         163             :          call reconstruction(fcube,nlev,k,&
         164             :              ctracer(:,:,:,:),irecons_tracer,llimiter,ntrac,&
         165             :              nc,nhe,nhr,nhc,nht,ns,nhr+(nhe-1),&
         166   483116400 :              fvm(ie)%jx_min,fvm(ie)%jx_max,fvm(ie)%jy_min,fvm(ie)%jy_max,&
         167             :              fvm(ie)%cubeboundary,fvm(ie)%halo_interp_weight,fvm(ie)%ibase,&
         168             :              spherecentroid,&
         169   483116400 :              fvm(ie)%recons_metrics,fvm(ie)%recons_metrics_integral,&
         170             :              fvm(ie)%rot_matrix,fvm(ie)%centroid_stretch,&
         171             :              fvm(ie)%vertex_recons_weights,fvm(ie)%vtx_cart,&
         172  1449349200 :              irecons_tracer_lev(k))
         173             :          !call t_stopf('FVM:tracers_reconstruct')
         174             :          !call t_startf('fvm:swept_flux')
         175   488311200 :          call swept_flux(elem(ie),fvm(ie),k,ctracer,irecons_tracer_lev(k),gsweights,gspts)
         176             :          !call t_stopf('fvm:swept_flux')
         177             :       end do
         178             :     end do
         179             :     !
         180             :     !***************************************
         181             :     !
         182             :     ! Large Courant number increment
         183             :     !
         184             :     !***************************************
         185             :     !
         186             :     ! In the jet region the effective Courant number
         187             :     ! in the cslam trajectory algorithm can be > 1
         188             :     ! (by up to 20%) in CAM
         189             :     !
         190             :     ! We limit the trajectories to < 1 but in this step
         191             :     ! we do a piecewise constant update for the
         192             :     ! amount of mass for which the Courant number is >1
         193             :     !
         194             :     !
         195      738816 :     if (large_Courant_incr) then
         196             :       !call t_startf('fvm:fill_halo_fvm:large_Courant')
         197             :       !if (kmin_jet<kmin.or.kmax_jet>kmax) then
         198             :       !  call endrun('ERROR: kmax_jet must be .le. kmax passed to run_consistent_se_cslam')
         199             :       !end if
         200             :       ! Determine the extent of the JET that is owned by this thread
         201      738816 :       ActiveJetThread = threadOwnsVertLevel(hybridnew,kmin_jet) .or. threadOwnsVertLevel(hybridnew,kmax_jet)
         202      738816 :       kmin_jet_local = max(kmin_jet,kmin)
         203      738816 :       kmax_jet_local = min(kmax_jet,kmax)
         204      738816 :       klev = kmax_jet-kmin_jet+1
         205      738816 :       call fill_halo_fvm(ghostbufQ1,elem,fvm,hybridnew,nets,nete,1,kmin_jet_local,kmax_jet_local,klev,active=ActiveJetThread)
         206             :       !call t_stopf('fvm:fill_halo_fvm:large_Courant')
         207             :       !call t_startf('fvm:large_Courant_number_increment')
         208      738816 :       if(ActiveJetThread) then
         209    69448704 :         do k=kmin_jet_local,kmax_jet_local !1,nlev
         210   552565104 :           do ie=nets,nete
         211   551826288 :             call large_courant_number_increment(fvm(ie),k)
         212             :           end do
         213             :         end do
         214             :       endif
         215             :       !call t_stopf('fvm:large_Courant_number_increment')
         216             :     end if
         217             :
         218             :     !call t_startf('fvm:end_of_reconstruct_subroutine')
         219    69448704 :     do k=kmin,kmax
         220             :       !
         221             :       ! convert to mixing ratio
         222             :       !
         223   552565104 :       do ie=nets,nete
         224  1932465600 :         do j=1,nc
         225  6280513200 :           do i=1,nc
         226  5797396800 :             inv_dp_area(i,j) = 1.0_r8/fvm(ie)%dp_fvm(i,j,k)
         227             :           end do
         228             :         end do
         229             :
         230 20290888800 :         do itr=1,ntrac
         231 79714206000 :           do j=1,nc
         232 >25750*10^7 :             do i=1,nc
         233             :               ! convert to mixing ratio
         234 >17826*10^7 :               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 >23769*10^7 :               fvm(ie)%c(i,j,k,itr) = MAX(fvm(ie)%c(i,j,k,itr),qmin(itr))
         237             :             end do
         238             :           end do
         239             :         end do
         240             :         !
         241             :         ! convert to dp and scale back dp
         242             :         !
         243  6280513200 :         fvm(ie)%dp_fvm(1:nc,1:nc,k) = fvm(ie)%dp_fvm(1:nc,1:nc,k)*fvm(ie)%dp_ref(k)*fvm(ie)%inv_area_sphere
         244             : #ifdef waccm_debug
         245             :         do j=1,nc
         246             :           do i=1,nc
         247             :             fvm(ie)%CSLAM_gamma(i,j,k,1) = MAXVAL(fvm(ie)%CSLAM_gamma(i,j,k,:))
         248             :           end do
         249             :         end do
         250             : #endif
         251 25673879088 :         elem(ie)%sub_elem_mass_flux(:,:,:,k)=0
         252             :       end do
         253             :     end do
         254             :     !call t_stopf('fvm:end_of_reconstruct_subroutine')
         255             :     !$OMP END PARALLEL
         256      738816 :     call omp_set_nested(.false.)
         257     1477632 :   end subroutine run_consistent_se_cslam
         258             :
         259   483116400 :   subroutine swept_flux(elem,fvm,ilev,ctracer,irecons_tracer_actual,gsweights,gspts)
         260      738816 :     use fvm_analytic_mod      , only: get_high_order_weights_over_areas
         261             :     use dimensions_mod, only : kmin_jet,kmax_jet
         262             :     implicit none
         263             :     type (element_t) , intent(in)   :: elem
         264             :     type (fvm_struct), intent(inout):: fvm
         265             :     integer          , intent(in) :: ilev, irecons_tracer_actual
         266             :     real (kind=r8), intent(inout) :: ctracer(irecons_tracer,1-nhe:nc+nhe,1-nhe:nc+nhe,ntrac)
         267             :     real (kind=r8), dimension(ngpc), intent(in) :: gsweights, gspts
         268             :     integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
         269             :     real (kind=r8)    , dimension(0:7       , imin:imax,imin:imax,num_sides) :: displ
         270             :     integer (kind=r8) , dimension(1:2,11    , imin:imax,imin:imax,num_sides) :: base_vec
         271             :     real (kind=r8)    , dimension(1:2, 6    , imin:imax,imin:imax,num_sides) :: base_vtx
         272             :     integer                  , dimension(2,num_area, imin:imax,imin:imax,num_sides) :: idx
         273             :     real (kind=r8)    , dimension(imin:imax,imin:imax,num_sides)             :: mass_flux_se
         274             :     real (kind=r8)    , dimension(irecons_tracer,num_area) :: weights
         275             :     real (kind=r8)                     :: gamma
         276             :     integer :: i,j,iside,iarea,iw
         277             :
         278             :     integer, parameter :: num_seg_max=5
         279             :     REAL(KIND=r8), dimension(2,num_seg_max,num_area) :: x, dx, x_static, dx_static
         280             :     integer             , dimension(num_area)               :: num_seg, num_seg_static
         281             :     REAL(KIND=r8), dimension(2,8) :: x_start, dgam_vec
         282             :     REAL(KIND=r8) :: gamma_max, displ_first_guess
         283             :
         284   966232800 :     REAL(KIND=r8) :: flux,flux_tracer(ntrac)
         285             :
         286             :     REAL(KIND=r8), dimension(num_area) :: dp_area
         287             :
         288             :     real (kind=r8) :: dp(1-nhc:nc+nhc,1-nhc:nc+nhc)
         289             :
         290             :     logical :: tl1,tl2,tr1,tr2
         291             :
         292             :     integer, dimension(4), parameter :: imin_side = (/1   ,0   ,1   ,1   /)
         293             :     integer, dimension(4), parameter :: imax_side = (/nc  ,nc  ,nc  ,nc+1/)
         294             :     integer, dimension(4), parameter :: jmin_side = (/1   ,1   ,0   ,1   /)
         295             :     integer, dimension(4), parameter :: jmax_side = (/nc+1,nc  ,nc  ,nc  /)
         296             :
         297             :     integer :: iseg, iseg_tmp,flowcase,ii,jj,itr
         298             :
         299   483116400 :     call define_swept_areas(fvm,ilev,displ,base_vec,base_vtx,idx)
         300             :
         301 25605169200 :     mass_flux_se(1:nc,1:nc,1:4)  = -elem%sub_elem_mass_flux(1:nc,1:nc,1:4,ilev)
         302  1932465600 :     mass_flux_se(0   ,1:nc,2  )  =  elem%sub_elem_mass_flux(1   ,1:nc,4  ,ilev)
         303  1932465600 :     mass_flux_se(nc+1,1:nc,4  )  =  elem%sub_elem_mass_flux(nc  ,1:nc,2  ,ilev)
         304  1932465600 :     mass_flux_se(1:nc,0   ,3  )  =  elem%sub_elem_mass_flux(1:nc,1   ,1  ,ilev)
         305  1932465600 :     mass_flux_se(1:nc,nc+1,1  )  =  elem%sub_elem_mass_flux(1:nc,nc  ,3  ,ilev)
         306             :     !
         307             :     ! prepare for air/tracer update
         308             :     !
         309             : !    dp = fvm%dp_fvm(1-nhe:nc+nhe,1-nhe:nc+nhe,ilev)
         310 43963592400 :     dp = fvm%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,ilev)
         311  6280513200 :     fvm%dp_fvm(1:nc,1:nc,ilev) = fvm%dp_fvm(1:nc,1:nc,ilev)*fvm%area_sphere
         312 20290888800 :     do itr=1,ntrac
         313 >25750*10^7 :       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 >13913*10^7 :       do iw=1,irecons_tracer_actual
         315 >11884*10^7 :         ctracer(iw,1-nhe:nc+nhe,1-nhe:nc+nhe,itr)=ctracer(iw,1-nhe:nc+nhe,1-nhe:nc+nhe,itr)*&
         316 >38229*10^8 :              dp(1-nhe:nc+nhe,1-nhe:nc+nhe)
         317             :       end do
         318             :     end do
         319             :
         320  2415582000 :     do iside=1,4
         321  9179211600 :       do j=jmin_side(iside),jmax_side(iside)
         322 31885682400 :         do i=imin_side(iside),imax_side(iside)
         323             :            !DO NOT USE MASS_FLUX_SE AS THRESHOLD - THRESHOLD CONDITION MUST BE CONSISTENT WITH
         324             :            !THE ONE USED IN DEFINE_SWEPT_AREAS
         325             : !          if (mass_flux_se(i,j,iside)>eps) then
         326 29953216800 :           if (fvm%se_flux(i,j,iside,ilev)>eps) then
         327             :             !
         328             :             !        ||             ||
         329             :             !  tl1   ||             || tr1
         330             :             !        ||             ||
         331             :             !  =============================
         332             :             !        ||             ||
         333             :             !  tl2   ||             || tr2
         334             :             !        ||             ||
         335             :             !
         336 11594793589 :             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 11594793589 :             tl2 = displ(6,i,j,iside)<0.0_r8.and.displ(7,i,j,iside)   >0.0_r8 !departure point in tl2 quadrant
         338 11594793589 :             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 11594793589 :             tr2 = displ(4,i,j,iside)<0.0_r8.and.displ(5,i,j,iside)   >0.0_r8 !departure point in tr2 quadrant
         340             :
         341             :             !
         342             :             ! pathological cases
         343             :             !
         344             :             !        |  ||           ||                      ||           ||
         345             :             !        |  ||-----------||                      ||-----------||
         346             :             !        |  ||           ||                      ||           ||
         347             :             !  ================================     =================================
         348             :             !           ||           ||                   |  ||           ||
         349             :             !  ---------||           ||             ------|--||           ||
         350             :             !           ||           ||                   |  ||           ||
         351             :             !
         352             :             !                tl1=tl1.or.tl2
         353             :             !                tr1=tr1.or.tr2
         354             :             !                tl1=displ(3,i,j,iside)<0.0_r8.and..not.(tl1.and.tl2)
         355             :             !                tr1=displ(2,i,j,iside)<0.0_r8.and..not.(tr1.and.tr2)
         356             :
         357 >13913*10^7 :             num_seg=-1; num_seg_static=-1 !initialization
         358 11594793589 :             if (.not.tl1.and..not.tl2.and..not.tr1.and..not.tr2) then
         359   212983856 :               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   212983856 :                    num_seg_static,x_start, dgam_vec,fvm%se_flux(i,j,iside,ilev),displ_first_guess)
         371             :
         372   212983856 :               gamma=1.0_r8!fvm%se_flux(i,j,iside,ilev)
         373   212983856 :               gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         374             :             else
         375 11381809733 :               if (tl1.and.tr1) then
         376   174458550 :                 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   174458550 :                      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   174458550 :                      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   174458550 :                      num_seg, num_seg_static,x_start, dgam_vec)
         392   174458550 :                 gamma=1.0_r8
         393   174458550 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         394 11207351183 :               else if (tl1.and..not.tr1.and..not.tr2) then
         395  5380007344 :                 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  5380007344 :                      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  5380007344 :                      x_start, dgam_vec)
         408  5380007344 :                 gamma=1.0_r8
         409  5380007344 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         410  5827343839 :               else if (tr1.and..not.tl1.and..not.tl2) then !displ(3).ge.0.0_r8) then
         411  5377464893 :                 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  5377464893 :                      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  5377464893 :                      num_seg_static, x_start, dgam_vec,displ_first_guess)
         425  5377464893 :                 gamma=1.0_r8
         426  5377464893 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         427   449878946 :               else if (tl2.and..not.tr1.and..not.tr2) then !displ(2).ge.0.0_r8) then
         428   218162473 :                 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   218162473 :                      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   218162473 :                      x_start, dgam_vec,displ_first_guess)
         443   218162473 :                 gamma = 1.0_r8
         444   218162473 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         445   231716473 :               else if (tr2.and..not.tl1.and..not.tl2) then !displ(3).ge.0.0_r8) then
         446   219327193 :                 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   219327193 :                      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   219327193 :                      num_seg_static,x_start, dgam_vec,displ_first_guess)
         462   219327193 :                 gamma=1.0_r8
         463   219327193 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         464    12389280 :               else if (tl2.and.tr1.and..not.tr2) then
         465     6119836 :                 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     6119836 :                      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     6119836 :                      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     6119836 :                      num_seg_static,x_start, dgam_vec,displ_first_guess)
         484             :
         485     6119836 :                 gamma=1.0_r8
         486     6119836 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         487     6269444 :               else if (tr2.and.tl1.and..not.tl2) then
         488     6108861 :                 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     6108861 :                      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     6108861 :                      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     6108861 :                      num_seg_static,x_start, dgam_vec)
         507     6108861 :                 gamma =  1.0_r8
         508     6108861 :                 gamma_max = fvm%displ_max(i,j,iside)/displ_first_guess
         509      160583 :               else if (tl2.and.tr2) then
         510      160583 :                 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      160583 :                      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      160583 :                      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      160583 :                      num_seg_static,x_start, dgam_vec,displ_first_guess)
         532      160583 :                 gamma =  1.0_r8
         533      160583 :                 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 69568761534 :             do iarea=1,num_area
         543 69568761534 :               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 34784380767 :                  num_seg_max,num_area,dp_area,flowcase,gamma,mass_flux_se(i,j,iside),0.0_r8,gamma_max,          &
         547 34784380767 :                  gsweights,gspts,ilev)
         548             :             !call t_stopf('fvm:swept_area:get_gamma')
         549             :             !
         550             :             ! pack segments for high-order weights computation
         551             :             !
         552 69568761534 :             do iarea=1,num_area
         553 92758348712 :               do iseg=1,num_seg_static(iarea)
         554 34784380767 :                 iseg_tmp=num_seg(iarea)+iseg
         555 >10435*10^7 :                 x (:,iseg_tmp,iarea)  = x_static (:,iseg,iarea)
         556 >16232*10^7 :                 dx(:,iseg_tmp,iarea)  = dx_static(:,iseg,iarea)
         557             :               end do
         558 69568761534 :               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 11594793589 :                  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 >48698*10^7 :             flux=0.0_r8; flux_tracer=0.0_r8
         576 69568761534 :             do iarea=1,num_area
         577 69568761534 :               if (num_seg(iarea)>0) then
         578 23613490681 :                 ii=idx(1,iarea,i,j,iside); jj=idx(2,iarea,i,j,iside)
         579 23613490681 :                 flux=flux+weights(1,iarea)*dp(ii,jj)
         580 >99176*10^7 :                 do itr=1,ntrac
         581 >68006*10^8 :                   do iw=1,irecons_tracer_actual
         582 >67770*10^8 :                     flux_tracer(itr) = flux_tracer(itr)+weights(iw,iarea)*ctracer(iw,ii,jj,itr)
         583             :                   end do
         584             :                 end do
         585             :               end if
         586             :             end do
         587 11594793589 :             fvm%se_flux(i,j,iside,ilev) = mass_flux_se(i,j,iside)-flux
         588 11594793589 :             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 11594793589 :             fvm%dp_fvm(i  ,j  ,ilev        ) = fvm%dp_fvm(i  ,j  ,ilev        )-flux
         597 >48698*10^7 :             fvm%     c(i  ,j  ,ilev,1:ntrac) = fvm%     c(i  ,j  ,ilev,1:ntrac)-flux_tracer(1:ntrac)
         598             :             !
         599             :             ! update flux in nearest neighbor cells
         600             :             !
         601 11594793589 :             if (iside==1) then
         602  2936257703 :               fvm%dp_fvm(i,j-1,ilev        ) = fvm%dp_fvm(i,j-1,ilev        )+flux
         603 >12332*10^7 :               fvm%     c(i,j-1,ilev,1:ntrac) = fvm%     c(i,j-1,ilev,1:ntrac)+flux_tracer(1:ntrac)
         604             :             end if
         605 11594793589 :             if (iside==2) then
         606  2852956151 :               fvm%dp_fvm(i+1,j,ilev        ) = fvm%dp_fvm(i+1,j,ilev        )+flux
         607 >11982*10^7 :               fvm%     c(i+1,j,ilev,1:ntrac) = fvm%     c(i+1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
         608             :             end if
         609 11594793589 :             if (iside==3) then
         610  2861139088 :               fvm%dp_fvm(i,j+1,ilev        ) = fvm%dp_fvm(i,j+1,ilev        )+flux
         611 >12016*10^7 :               fvm%     c(i,j+1,ilev,1:ntrac) = fvm%     c(i,j+1,ilev,1:ntrac)+flux_tracer(1:ntrac)
         612             :             end if
         613 11594793589 :             if (iside==4) then
         614  2944440647 :               fvm%dp_fvm(i-1,j,ilev        ) = fvm%dp_fvm(i-1,j,ilev        )+flux
         615 >12366*10^7 :               fvm%     c(i-1,j,ilev,1:ntrac) = fvm%     c(i-1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
         616             :             end if
         617             :             !call t_stopf('fvm:swept_area:remap')
         618             :           end if
         619             :         end do
         620             :       end do
         621             :     end do
         622   483116400 :   end subroutine swept_flux
         623             :
         624             :
         625   483116400 :   subroutine large_courant_number_increment(fvm,ilev)
         626             :     implicit none
         627             :     type (fvm_struct), intent(inout):: fvm
         628             :     integer          , intent(in) :: ilev
         629             :
         630             :     integer, parameter :: num_sides=4, imin= 0, imax=nc+1
         631             :
         632             :     integer, dimension(4), parameter :: imin_side = (/1   ,0   ,1   ,1   /)
         633             :     integer, dimension(4), parameter :: imax_side = (/nc  ,nc  ,nc  ,nc+1/)
         634             :     integer, dimension(4), parameter :: jmin_side = (/1   ,1   ,0   ,1   /)
         635             :     integer, dimension(4), parameter :: jmax_side = (/nc+1,nc  ,nc  ,nc  /)
         636             :
         637             :     integer :: i,j,iside,itr
         638   966232800 :     real (kind=r8)    :: flux,flux_tracer(ntrac)
         639             :     real (kind=r8), dimension(0:nc+1,0:nc+1)      :: inv_dp_area
         640   966232800 :     real (kind=r8), dimension(0:nc+1,0:nc+1,ntrac):: c_tmp
         641             :
         642 14976608400 :     inv_dp_area=1.0_r8/fvm%dp_fvm(0:nc+1,0:nc+1,ilev)
         643 >61452*10^7 :     c_tmp      = fvm%c(0:nc+1,0:nc+1,ilev,1:ntrac)
         644  2415582000 :     do iside=1,4
         645  9179211600 :       do j=jmin_side(iside),jmax_side(iside)
         646 31885682400 :         do i=imin_side(iside),imax_side(iside)
         647 29953216800 :           if (fvm%se_flux(i,j,iside,ilev)>eps) then
         648             :             flux = fvm%se_flux(i,j,iside,ilev)
         649             : #ifdef waccm_debug
         650             :             if (i>0.and.j>0.and.i<nc+1.and.j<nc+1) then
         651             :                fvm%CSLAM_gamma(i,j,ilev,iside) = fvm%CSLAM_gamma(i,j,ilev,iside)+&
         652             :                     fvm%se_flux(i,j,iside,ilev)*inv_dp_area(i,j)
         653             :             end if
         654             : #endif
         655             :
         656   297859926 :             do itr=1,ntrac
         657   297859926 :               flux_tracer(itr) = fvm%se_flux(i,j,iside,ilev)*c_tmp(i,j,itr)*inv_dp_area(i,j)
         658             :             end do
         659     7091903 :             fvm%dp_fvm(i  ,j  ,ilev        ) = fvm%dp_fvm(i  ,j  ,ilev        )-flux
         660   297859926 :             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     7091903 :             if (iside==1) then
         665      278999 :               fvm%dp_fvm(i,j-1,ilev        ) = fvm%dp_fvm(i,j-1,ilev        )+flux
         666    11717958 :               fvm%     c(i,j-1,ilev,1:ntrac) = fvm%     c(i,j-1,ilev,1:ntrac)+flux_tracer(1:ntrac)
         667             :             end if
         668     7091903 :             if (iside==2) then
         669     4336167 :               fvm%dp_fvm(i+1,j,ilev        ) = fvm%dp_fvm(i+1,j,ilev        )+flux
         670   182119014 :               fvm%     c(i+1,j,ilev,1:ntrac) = fvm%     c(i+1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
         671             :             end if
         672     7091903 :             if (iside==3) then
         673      632290 :               fvm%dp_fvm(i,j+1,ilev        ) = fvm%dp_fvm(i,j+1,ilev        )+flux
         674    26556180 :               fvm%     c(i,j+1,ilev,1:ntrac) = fvm%     c(i,j+1,ilev,1:ntrac)+flux_tracer(1:ntrac)
         675             :             end if
         676     7091903 :             if (iside==4) then
         677     1844447 :               fvm%dp_fvm(i-1,j,ilev        ) = fvm%dp_fvm(i-1,j,ilev        )+flux
         678    77466774 :               fvm%     c(i-1,j,ilev,1:ntrac) = fvm%     c(i-1,j,ilev,1:ntrac)+flux_tracer(1:ntrac)
         679             :             end if
         680             :           end if
         681             :         end do
         682             :       end do
         683             :     end do
         684   483116400 :   end subroutine large_courant_number_increment
         685             :
         686   483116400 :   subroutine ghost_flux_unpack(fvm,var)
         687             :     use control_mod, only : neast, nwest, seast, swest
         688             :     implicit none
         689             :     type (fvm_struct), intent(inout) :: fvm
         690             :     real(kind=r8)                    :: var(1-nhe:nc+nhe,1-nhe:nc+nhe,4)
         691             :
         692             :     integer :: i,j,ishft
         693             :     !
         694             :     ! rotate coordinates if needed
         695             :     !
         696   483116400 :     if (fvm%cubeboundary.NE.0) then
         697   373610016 :       do j=1-nhe,nc+nhe
         698  1930318416 :         do i=1-nhe,nc+nhe
         699  1556708400 :           ishft = NINT(fvm%flux_orient(2,i,j))
         700  8094883680 :           var(i,j,1:4) = cshift(var(i,j,1:4),shift=ishft)
         701             :         end do
         702             :       end do
         703             :       !
         704             :       ! non-existent cells in physical space - necessary?
         705             :       !
         706    62268336 :       if (fvm%cubeboundary==nwest) then
         707     6978348 :         var(1-nhe:0,nc+1 :nc+nhe,:) = 0.0_r8
         708    61731540 :       else if (fvm%cubeboundary==swest) then
         709     6978348 :         var(1-nhe:0,1-nhe:0     ,:) = 0.0_r8
         710    61194744 :       else if (fvm%cubeboundary==neast) then
         711     6978348 :         var(nc+1 :nc+nhe,nc+1 :nc+nhe,:) = 0.0_r8
         712    60657948 :       else if (fvm%cubeboundary==seast) then
         713     6978348 :         var(nc+1 :nc+nhe,1-nhe:0,:) = 0.0_r8
         714             :       end if
         715             :     end if
         716   483116400 :   end subroutine ghost_flux_unpack
         717             :
         718   483116400 :   subroutine compute_displacements_for_swept_areas(fvm,cair,k,gsweights,gspts)
         719             :     use dimensions_mod, only: large_Courant_incr
         720             :     implicit none
         721             :     type (fvm_struct), intent(inout)     :: fvm
         722             :     integer, intent(in) :: k
         723             :     real (kind=r8)      :: cair(1-nhc:nc+nhc,1-nhc:nc+nhc) !high-order air density reconstruction
         724             :     real (kind=r8), dimension(ngpc), intent(in) :: gsweights, gspts
         725             :     !
         726             :     !   flux iside 1                     flux iside 3                    flux iside 2       flux iside 4
         727             :     !
         728             :     !   |          |                     |  ---1--> |                    |    --2-->|       |--1-->    |
         729             :     !  -4----------3-   /\              -4----------3-                  -4----------3-     -4----------3-   ||
         730             :     !   |          |   /||\              |\\\\\\\\\\|    ||              |   |\\\\\\|       |\\\\\\|   |
         731             :     !   |  --2-->  |    || dv(1)         |\\\\\\\\\\|    ||              |   |\\\\\\|       |\\\\\\|   |
         732             :     !   |----------|    ||               |----------|    || dv(3)        |   |\\\\\\|       |\\\\\\|   |
         733             :     !   |\\\\\\\\\\|    ||               | <--2---  |   \||/             |   |\\\\\\|       |\\\\\\|   |
         734             :     !   |\\\\\\\\\\|    ||               |          |    \/              |   |\\\\\\|       |\\\\\\|   |
         735             :     !  -1----------2-                   -1----------2-                  -1----------2-     -1----------2-
         736             :     !   |  <--1--  |                     |          |                    |    <--1--|       |<--2--
         737             :     !
         738             :     !                                                                     /                          \
         739             :     !   line-integral                                                    <==========         =========>
         740             :     !   from vertex 2                                                     \  dv(2)              dv(4)/
         741             :     !   to 1
         742             :     !
         743             :     !   Note vertical
         744             :     !   lines have
         745             :     !   zero line-
         746             :     !   integral!
         747             :     !
         748             :     integer               :: i,j,iside,ix
         749             :     integer, parameter :: num_area=1, num_seg_max=2
         750             :     REAL(KIND=r8), dimension(2,num_seg_max,num_area,4,nc,nc) :: x_static, dx_static
         751             :     REAL(KIND=r8), dimension(2,num_seg_max,num_area,4,nc,nc) :: x, dx
         752             :     REAL(KIND=r8), dimension(2,num_seg_max,num_area)         :: x_tmp, dx_tmp
         753             :     integer             , dimension(              num_area,4      ) :: num_seg, num_seg_static
         754             :     REAL(KIND=r8), dimension(2,8,                   4,nc,nc) :: x_start, dgam_vec
         755             :     REAL(KIND=r8), dimension(num_area) :: dp_area
         756             :     integer, dimension(4) :: flowcase
         757             :     REAL(KIND=r8)  :: gamma(4), flux_se
         758             :
         759   483116400 :     num_seg_static(1,1) =  1; num_seg(1,1) = 1; flowcase(1) = -1
         760   483116400 :     num_seg_static(1,2) =  0; num_seg(1,2) = 2; flowcase(2) = -2
         761   483116400 :     num_seg_static(1,3) =  1; num_seg(1,3) = 1; flowcase(3) = -1
         762   483116400 :     num_seg_static(1,4) =  0; num_seg(1,4) = 2; flowcase(4) = -4
         763             :
         764  1932465600 :     do j=1,nc
         765  6280513200 :        do i=1,nc
         766 14493492000 :           do ix=1,2
         767  8696095200 :              iside=1;
         768  8696095200 :              x_static (ix,1,1,iside,i,j) = fvm%vtx_cart(2,ix,i,j)
         769  8696095200 :              dx_static(ix,1,1,iside,i,j) = fvm%vtx_cart(1,ix,i,j)-fvm%vtx_cart(2,ix,i,j)
         770  8696095200 :              x_start  (ix,1,  iside,i,j) = fvm%vtx_cart(1,ix,i,j)
         771  8696095200 :              x_start  (ix,2,  iside,i,j) = fvm%vtx_cart(2,ix,i,j)
         772  8696095200 :              dgam_vec (ix,1,  iside,i,j) = fvm%vtx_cart(4,ix,i,j)-fvm%vtx_cart(1,ix,i,j)
         773             :              !
         774             :              ! compute first guess
         775             :              !
         776  8696095200 :              gamma(iside)                       = 0.5_r8
         777  8696095200 :              x        (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
         778  8696095200 :              dx       (ix,1,1,iside,i,j) = -dx_static(ix,1,1,iside,i,j)
         779             :              !
         780             :              ! side 2
         781             :              !
         782  8696095200 :              iside=2;
         783  8696095200 :              x_start  (ix,1,  iside,i,j) = fvm%vtx_cart(2,ix,i,j)
         784  8696095200 :              x_start  (ix,2,  iside,i,j) = fvm%vtx_cart(3,ix,i,j)
         785  8696095200 :              dgam_vec (ix,1,  iside,i,j) = fvm%vtx_cart(1,ix,i,j)-fvm%vtx_cart(2,ix,i,j)
         786  8696095200 :              x        (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)
         787             :              !
         788             :              ! compute first guess - gamma=1
         789             :              !
         790  8696095200 :              gamma(iside)                       = 0.5_r8
         791  8696095200 :              dx       (ix,1,1,iside,i,j) =  gamma(iside)*dgam_vec (ix,1,  iside,i,j)
         792  8696095200 :              x        (ix,2,1,iside,i,j) =  x_start(ix,2,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
         793  8696095200 :              dx       (ix,2,1,iside,i,j) = -gamma(iside)*dgam_vec (ix,1,  iside,i,j)
         794             :              !
         795             :              ! side 3
         796             :              !
         797  8696095200 :              iside=3;
         798  8696095200 :              x_static (ix,1,1,iside,i,j) = fvm%vtx_cart(4,ix,i,j)
         799  8696095200 :              dx_static(ix,1,1,iside,i,j) = fvm%vtx_cart(3,ix,i,j)-fvm%vtx_cart(4,ix,i,j)
         800  8696095200 :              x_start  (ix,1,  iside,i,j) = fvm%vtx_cart(3,ix,i,j)
         801  8696095200 :              x_start  (ix,2,  iside,i,j) = fvm%vtx_cart(4,ix,i,j)
         802  8696095200 :              dgam_vec (ix,1,  iside,i,j) = fvm%vtx_cart(2,ix,i,j)-fvm%vtx_cart(3,ix,i,j)
         803             :              !
         804             :              ! compute first guess - gamma(iside)=1
         805             :              !
         806  8696095200 :              gamma(iside)                       = 0.5_r8
         807  8696095200 :              x        (ix,1,1,iside,i,j) = x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
         808  8696095200 :              dx       (ix,1,1,iside,i,j) = -dx_static(ix,1,1,iside,i,j)
         809             :              !
         810             :              ! side 4
         811             :              !
         812  8696095200 :              iside=4;
         813  8696095200 :              x_start  (ix,1,  iside,i,j) = fvm%vtx_cart(1,ix,i,j)
         814  8696095200 :              x_start  (ix,2,  iside,i,j) = fvm%vtx_cart(4,ix,i,j)
         815  8696095200 :              dgam_vec (ix,1,  iside,i,j) = fvm%vtx_cart(2,ix,i,j)-fvm%vtx_cart(1,ix,i,j)
         816  8696095200 :              x        (ix,2,1,iside,i,j) = x_start(ix,2,iside,i,j)
         817             :              !
         818             :              ! compute first guess - gamma(iside)=1
         819             :              !
         820  8696095200 :              gamma(iside)                       = 0.5_r8
         821  8696095200 :              dx       (ix,2,1,iside,i,j) =  gamma(iside)*dgam_vec (ix,1,  iside,i,j)
         822  8696095200 :              x        (ix,1,1,iside,i,j) =  x_start(ix,1,iside,i,j)+gamma(iside)*dgam_vec(ix,1,iside,i,j)
         823 13044142800 :              dx       (ix,1,1,iside,i,j) = -gamma(iside)*dgam_vec (ix,1,  iside,i,j)
         824             :           end do
         825             :        end do
         826             :     end do
         827             :
         828             : !    do k=1,nlev
         829  1932465600 :       do j=1,nc
         830  6280513200 :         do i=1,nc
         831  8696095200 :           dp_area = cair(i,j)
         832 23189587200 :           do iside=1,4
         833 17392190400 :             flux_se = -fvm%se_flux(i,j,iside,k)
         834 21740238000 :             if (flux_se>eps) then
         835  8696095193 :               gamma(iside)=0.5_r8
         836             :               !
         837             :               ! this copying is necessary since get_flux_segments_area_iterate change x and dx
         838             :               !
         839 56527099828 :               x_tmp (:,1:num_seg(1,iside),:)=x (:,1:num_seg(1,iside),:,iside,i,j)
         840 56527099828 :               dx_tmp(:,1:num_seg(1,iside),:)=dx(:,1:num_seg(1,iside),:,iside,i,j)
         841             :               call get_flux_segments_area_iterate(&
         842             :                    x_tmp(:,:,:),x_static(:,:,:,iside,i,j),dx_static(:,:,:,iside,i,j),dx_tmp(:,:,:),&
         843             :                    x_start(:,:,iside,i,j),dgam_vec(:,:,iside,i,j),num_seg(:,iside),num_seg_static(:,iside),&
         844             :                    num_seg_max,num_area,dp_area,flowcase(iside),gamma(iside),flux_se,0.0_r8,1.0_r8,        &
         845  8696095193 :                    gsweights,gspts,k)
         846 26088285579 :               fvm%se_flux(i,j,iside,k) = ABS(SUM(gamma(iside)*dgam_vec(:,1,iside,i,j)))
         847             : #ifdef waccm_debug
         848             :               fvm%CSLAM_gamma(i,j,k,iside) = gamma(iside)
         849             : #endif
         850  8696095193 :               if (gamma(iside)>1_r8) then
         851           0 :                  if (.not.large_Courant_incr) then
         852           0 :                     write(iulog,*) 'ERROR in CSLAM: local Courant number is >1: gamma=',gamma(iside),' k=',k
         853           0 :                     call endrun('ERROR in CSLAM: local Courant number is > 1; set namelist se_large_Courant_incr=.true. ')
         854             :                  endif
         855           0 :                 gamma(iside)=1.0_r8-eps
         856             :               end if
         857             :             else
         858  8696095207 :               fvm%se_flux(i,j,iside,k) = 0.0_r8
         859             : #ifdef waccm_debug
         860             :               fvm%CSLAM_gamma(i,j,k,iside) = 0.0_r8
         861             : #endif
         862             :             end if
         863             :           enddo
         864             :         end do
         865             :       end do
         866             : !    end do
         867   483116400 :   end subroutine compute_displacements_for_swept_areas
         868             :
         869             :
         870             :
         871 20290888782 :   subroutine get_flux_segments_area_iterate(x,x_static,dx_static,dx,x_start,dgam_vec,num_seg,num_seg_static,&
         872 20290888782 :        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 40581777564 :     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 20290888782 :     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 20290888782 :     lexit_after_one_more_iteration = .false.
         900             :     !
         901             :     ! compute static line-integrals (not necessary to recompute them for every iteration)
         902             :     !
         903 20290888782 :     flux_static = 0.0_r8
         904 86960951920 :     w_static    = 0.0_r8
         905 86960951920 :     weight_area = 0.0_r8
         906 86960951920 :     do iarea=1,num_area
         907 >10580*10^7 :        do iseg=1,num_seg_static(iarea)
         908             :
         909             : !rck vector directive needed here
         910             : !DIR$ SIMD
         911 >15652*10^7 :           do ipt=1,ngpc
         912 >11739*10^7 :              xq(ipt) = x_static(1,iseg,iarea)+dx_static(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
         913 >11739*10^7 :              yq(ipt) = x_static(2,iseg,iarea)+dx_static(2,iseg,iarea)*gspts(ipt)
         914 >15652*10^7 :              F(ipt,1) = yq(ipt)/(SQRT(1.0_r8+xq(ipt)*xq(ipt) + yq(ipt)*yq(ipt))*(1.0_r8+xq(ipt)*xq(ipt)))! potential ! potential
         915             :           enddo
         916 >22319*10^7 :           weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx_static(1,iseg,iarea) !integral
         917             :        end do
         918 66670063138 :        w_static(iarea)= weight_area(iarea)
         919 86960951920 :        flux_static = flux_static+weight_area(iarea)*c(iarea)      !add to swept flux
         920             :     end do
         921             :     !
         922             :     ! initilization
         923             :     !
         924 20290888782 :     gamma1=0.0_r8; f1=-flux   ! zero flux guess 1
         925             :     !
         926             :     ! compute flux integrals of first guess passed to subroutine
         927             :     !
         928 20290888782 :     gamma2=gamma
         929 20290888782 :     f2 = 0.0_r8
         930 86960951920 :     weight_area=w_static
         931 86960951920 :     do iarea=1,num_area
         932 >12757*10^7 :        do iseg=1,num_seg(iarea)
         933             : !rck vector directive needed here
         934             : !DIR$ SIMD
         935 >24363*10^7 :           do ipt=1,ngpc
         936 >18272*10^7 :              xq(ipt)  = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
         937 >18272*10^7 :              yq(ipt)  = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(ipt)
         938 >18272*10^7 :              xq2      =  xq(ipt)*xq(ipt)
         939 >18272*10^7 :              xq2i     =  1.0_r8/(1.0_r8+xq2)
         940 >18272*10^7 :              rho      =  SQRT(1.0_r8+xq2+yq(ipt)*yq(ipt))
         941 >18272*10^7 :              rhoi     =  1.0_r8/rho
         942 >18272*10^7 :              yrh      =  yq(ipt)*rhoi
         943 >24363*10^7 :              F(ipt,1) =  yrh*xq2i
         944             :           enddo
         945 >31030*10^7 :           weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx(1,iseg,iarea)! integral
         946             :        end do
         947 86960951920 :        f2 = f2+weight_area(iarea)*c(iarea)
         948             :     end do
         949 20290888782 :     f2 = f2-flux !integral error
         950 20290888782 :     iter=0
         951 20290888782 :     if (abs(f2-f1)<eps) then
         952             :       !
         953             :       ! in case the first guess is converged
         954             :       !
         955             :       return
         956             :     end if
         957             :
         958             :
         959 20290888782 :     dgamma=(gamma2-gamma1)*f2/(f2-f1);
         960 20290888782 :     gamma3 = gamma2-dgamma;                    ! Newton "guess" for gamma
         961 20290888782 :     gamma1 = gamma2; f1 = f2; gamma2 = gamma3; ! prepare for iteration
         962 57430365546 :     do iter=1,iter_max
         963             :        !
         964             :        ! update vertex location: flow_case dependent to avoid many zero operations
         965             :        !
         966 64410694657 :        select case(flow_case)
         967             :        case(-4)
         968  6980329111 :           iarea=1
         969 20940987333 :           dx       (:,2,1) =  gamma3*dgam_vec (:,1)
         970 20940987333 :           x        (:,1,1) =  x_start(:,1)+gamma3*dgam_vec(:,1)
         971 20940987333 :           dx       (:,1,1) = -gamma3*dgam_vec (:,1)
         972             :
         973             :        case(-2)
         974  6665394593 :           iarea=1
         975 19996183779 :           dx       (:,1,iarea) =  gamma3*dgam_vec (:,1)
         976 19996183779 :           x        (:,2,iarea) =  x_start(:,2)+gamma3*dgam_vec(:,1)
         977 19996183779 :           dx       (:,2,iarea) = -gamma3*dgam_vec (:,1)
         978             :        case(-1)
         979             :           !
         980             :           ! to compute first-guess perpendicular displacements for iside=1
         981             :           !
         982 13447426355 :           iarea=1
         983 40342279065 :           x        (:,1,iarea) = x_start(:,1)+gamma3*dgam_vec(:,1)
         984 40342279065 :           dx       (:,1,iarea) = -dx_static(:,1,iarea)
         985 40342279065 :           x        (:,2,iarea) = x_start(:,2)+gamma3*dgam_vec(:,1)
         986 40342279065 :           dx       (:,2,iarea) = x_start(:,2)-x(:,2,iarea)
         987             :        case(0)
         988   592002766 :           iarea=3
         989  1776008298 :           xtmp = x_start(:,1)+gamma3*dgam_vec(:,1)
         990  1776008298 :           dx       (:,1,iarea) = xtmp(:  )-x(:,1,iarea)           !dynamic - line 2
         991  1776008298 :           x        (:,2,iarea) = xtmp(:  )                        !dynamic - line 3
         992  1776008298 :           dx       (:,2,iarea) = x_static(:,2,iarea)-x(:,2,iarea) !dynamic - line 3
         993             :        case(1)
         994   372714989 :           iarea=2
         995  1118144967 :           xtmp(:        ) = x_start(:,1)+gamma3*dgam_vec(:,1)
         996  1118144967 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)        !dynamic - line 2
         997  1118144967 :           x   (:,2,iarea) = xtmp(:)                     !dynamic  - line 3
         998  1118144967 :           dx  (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
         999             :
        1000   372714989 :           iarea            = 3
        1001  1118144967 :           xtmp (:  )       = x_start(:,4)+gamma3*dgam_vec(:,4)
        1002  1118144967 :           xtmp2(:  )       = x_start(:,5)+gamma3*dgam_vec(:,5)
        1003  1118144967 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1004  1118144967 :           x    (:,2,iarea) = xtmp (:)         !dynamic
        1005  1118144967 :           dx   (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
        1006  1118144967 :           x    (:,3,iarea) = xtmp2(:)              !dynamic
        1007  1118144967 :           dx   (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
        1008             :
        1009   372714989 :           iarea         = 4
        1010  1118144967 :           xtmp    (:  ) = x_start(:,6)+gamma3*dgam_vec(:,6)
        1011  1118144967 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic - line 2
        1012  1118144967 :           x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1013  1118144967 :           dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1014             :        case(2)
        1015 13981278031 :           iarea=2
        1016 41943834093 :           xtmp(:        ) = x_start(:,1)+gamma3*dgam_vec(:,1)
        1017 41943834093 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)        !dynamic - line 2
        1018 41943834093 :           x   (:,2,iarea) = xtmp(:)                     !dynamic  - line 3
        1019 41943834093 :           dx  (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
        1020             :
        1021 13981278031 :           iarea=3
        1022 41943834093 :           xtmp(:        ) = x_start(:,4)+gamma3*dgam_vec(:,4)!
        1023 41943834093 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)        !dynamic - line 1
        1024 41943834093 :           x   (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1025 41943834093 :           dx  (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1026             :        case(3)
        1027 13958667257 :           iarea         = 3
        1028 41876001771 :           xtmp    (:  ) = x_start(:,5)+gamma3*dgam_vec(:,5)
        1029 41876001771 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
        1030 41876001771 :           x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1031 41876001771 :           dx       (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
        1032             :
        1033 13958667257 :           iarea         = 4
        1034 41876001771 :           xtmp    (:  ) = x_start(:,6)+gamma3*dgam_vec(:,6)
        1035 41876001771 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic - line 2
        1036 41876001771 :           x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1037 41876001771 :           dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1038             :        case(4)
        1039   701614480 :           iarea           = 1
        1040  2104843440 :           xtmp(:        ) = x_start(:,1)+gamma3*dgam_vec(:,1)
        1041  2104843440 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1042  2104843440 :           x (:,2,iarea) = xtmp(:)                      !dynamic
        1043  2104843440 :           dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:)  !dynamic
        1044             :
        1045   701614480 :           iarea         = 2
        1046  2104843440 :           xtmp    (:  ) = x_start(:,2)+gamma3*dgam_vec(:,2)
        1047  2104843440 :           xtmp2   (:  ) = x_start(:,3)+gamma3*dgam_vec(:,3)
        1048             :
        1049  2104843440 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic
        1050             :
        1051  2104843440 :           x (:,2,iarea) = xtmp (:)          !dynamic
        1052  2104843440 :           dx(:,2,iarea) = xtmp2(:)-xtmp(:)  !dynamic
        1053             :
        1054  2104843440 :           x (:,3,iarea) = xtmp2(:)                !dynamic
        1055  2104843440 :           dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:)   !dynamic
        1056             :
        1057   701614480 :           iarea            = 3
        1058  2104843440 :           xtmp (:        ) = x_start(:,4)+gamma3*dgam_vec(:,4)
        1059  2104843440 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic - line 1
        1060  2104843440 :           x    (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1061  2104843440 :           dx   (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1062             :        case(5)
        1063   705661840 :           iarea                = 3
        1064  2116985520 :           xtmp    (:  )        = x_start(:,5)+gamma3*dgam_vec(:,5)
        1065  2116985520 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
        1066  2116985520 :           x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1067  2116985520 :           dx       (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
        1068             :
        1069   705661840 :           iarea         = 4
        1070  2116985520 :           xtmp    (:  ) = x_start(:,6)+gamma3*dgam_vec(:,6)
        1071  2116985520 :           xtmp2   (:  ) = x_start(:,7)+gamma3*dgam_vec(:,7)
        1072             :
        1073  2116985520 :           dx(:,1,iarea) = xtmp(:)-x(:,1,iarea)   !dynamic - line 1
        1074  2116985520 :           x (:,2,iarea) = xtmp(:)          !dynamic -line 2
        1075  2116985520 :           dx       (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic - line 2
        1076  2116985520 :           x        (:,3,iarea) = xtmp2(:)               !dynamic  -line 1
        1077  2116985520 :           dx       (:,3,iarea) = x(:,1,iarea)-xtmp2(:)  !dynamic - line 1
        1078             :
        1079   705661840 :           iarea             = 5
        1080  2116985520 :           xtmp  (:  )       = x_start(:,8)+gamma3*dgam_vec(:,8)
        1081             :
        1082  2116985520 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)   !dynamic - line 1
        1083  2116985520 :           x        (:,2,iarea) = xtmp(:)                     !dynamic -line 2
        1084  2116985520 :           dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1085             :        case(6)
        1086    12502205 :           iarea = 1
        1087    37506615 :           xtmp(:  ) = x_start(:,1)+gamma3*dgam_vec(:,1)
        1088    37506615 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1089    37506615 :           x (:,2,iarea) = xtmp(:)                      !dynamic
        1090    37506615 :           dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:)  !dynamic
        1091             :
        1092    12502205 :           iarea         = 2
        1093    37506615 :           xtmp    (:  ) = x_start(:,2)+gamma3*dgam_vec(:,2)
        1094    37506615 :           xtmp2   (:  ) = x_start(:,3)+gamma3*dgam_vec(:,3)
        1095             :
        1096    37506615 :           dx(:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic
        1097    37506615 :           x (:,2,iarea) = xtmp (:)          !dynamic
        1098    37506615 :           dx(:,2,iarea) = xtmp2(:)-xtmp(:)  !dynamic
        1099    37506615 :           x (:,3,iarea) = xtmp2(:)                !dynamic
        1100    37506615 :           dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:)   !dynamic
        1101             :
        1102    12502205 :           iarea            = 3
        1103    37506615 :           xtmp (:  )       = x_start(:,4)+gamma3*dgam_vec(:,4)
        1104    37506615 :           xtmp2(:  )       = x_start(:,5)+gamma3*dgam_vec(:,5)
        1105    37506615 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1106    37506615 :           x    (:,2,iarea) = xtmp (:)         !dynamic
        1107    37506615 :           dx   (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
        1108    37506615 :           x    (:,3,iarea) = xtmp2(:)              !dynamic
        1109    37506615 :           dx   (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
        1110             :
        1111    12502205 :           iarea         = 4
        1112    37506615 :           xtmp    (:  ) = x_start(:,6)+gamma3*dgam_vec(:,6)
        1113    37506615 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic - line 2
        1114    37506615 :           x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1115    37506615 :           dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1116             :        case(7)
        1117    12446079 :           iarea=2
        1118    37338237 :           xtmp(:        ) = x_start(:,1)+gamma3*dgam_vec(:,1)
        1119    37338237 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)        !dynamic - line 2
        1120    37338237 :           x   (:,2,iarea) = xtmp(:)                     !dynamic  - line 3
        1121    37338237 :           dx  (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
        1122             :
        1123    12446079 :           iarea            = 3
        1124    37338237 :           xtmp (:  )       = x_start(:,4)+gamma3*dgam_vec(:,4)
        1125    37338237 :           xtmp2(:  )       = x_start(:,5)+gamma3*dgam_vec(:,5)
        1126    37338237 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1127    37338237 :           x    (:,2,iarea) = xtmp (:)         !dynamic
        1128    37338237 :           dx   (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
        1129    37338237 :           x    (:,3,iarea) = xtmp2(:)              !dynamic
        1130    37338237 :           dx   (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
        1131             :
        1132    12446079 :           iarea      = 4
        1133    37338237 :           xtmp    (:  ) = x_start(:,6)+gamma3*dgam_vec(:,6)
        1134    37338237 :           xtmp2   (:  ) = x_start(:,7)+gamma3*dgam_vec(:,7)
        1135             :
        1136    37338237 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
        1137    37338237 :           x        (:,2,iarea) = xtmp(:)              !dynamic
        1138    37338237 :           dx       (:,2,iarea) = xtmp2(:)-xtmp(:)     !dynamic
        1139    37338237 :           x        (:,3,iarea) = xtmp2(:)               !dynamic
        1140    37338237 :           dx       (:,3,iarea) = x(:,1,iarea)-xtmp2(:)  !dynamic
        1141             :
        1142    12446079 :           iarea      = 5
        1143    37338237 :           xtmp (:  ) = x_start(:,8)+gamma3*dgam_vec(:,8)
        1144    37338237 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)   !dynamic - line 1
        1145    37338237 :           x    (:,2,iarea) = xtmp(:)                     !dynamic -line 2
        1146    37338237 :           dx   (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1147             :        case(8)
        1148      327840 :           iarea = 1
        1149      983520 :           xtmp(:  ) = x_start(:,1)+gamma3*dgam_vec(:,1)
        1150      983520 :           dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1151      983520 :           x (:,2,iarea) = xtmp(:)                      !dynamic
        1152      983520 :           dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:)  !dynamic
        1153             :
        1154      327840 :           iarea         = 2
        1155      983520 :           xtmp    (:  ) = x_start(:,2)+gamma3*dgam_vec(:,2)
        1156      983520 :           xtmp2   (:  ) = x_start(:,3)+gamma3*dgam_vec(:,3)
        1157             :
        1158      983520 :           dx(:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic
        1159      983520 :           x (:,2,iarea) = xtmp (:)          !dynamic
        1160      983520 :           dx(:,2,iarea) = xtmp2(:)-xtmp(:)  !dynamic
        1161      983520 :           x (:,3,iarea) = xtmp2(:)                !dynamic
        1162      983520 :           dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:)   !dynamic
        1163             :
        1164      327840 :           iarea            = 3
        1165      983520 :           xtmp (:  )       = x_start(:,4)+gamma3*dgam_vec(:,4)
        1166      983520 :           xtmp2(:  )       = x_start(:,5)+gamma3*dgam_vec(:,5)
        1167      983520 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1168      983520 :           x    (:,2,iarea) = xtmp (:)         !dynamic
        1169      983520 :           dx   (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
        1170      983520 :           x    (:,3,iarea) = xtmp2(:)              !dynamic
        1171      983520 :           dx   (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
        1172             :
        1173      327840 :           iarea      = 4
        1174      983520 :           xtmp    (:  ) = x_start(:,6)+gamma3*dgam_vec(:,6)
        1175      983520 :           xtmp2   (:  ) = x_start(:,7)+gamma3*dgam_vec(:,7)
        1176             :
        1177      983520 :           dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic
        1178      983520 :           x        (:,2,iarea) = xtmp(:)              !dynamic
        1179      983520 :           dx       (:,2,iarea) = xtmp2(:)-xtmp(:)     !dynamic
        1180      983520 :           x        (:,3,iarea) = xtmp2(:)               !dynamic
        1181      983520 :           dx       (:,3,iarea) = x(:,1,iarea)-xtmp2(:)  !dynamic
        1182             :
        1183      327840 :           iarea      = 5
        1184      983520 :           xtmp (:  ) = x_start(:,8)+gamma3*dgam_vec(:,8)
        1185      983520 :           dx   (:,1,iarea) = xtmp(:)-x(:,1,iarea)   !dynamic - line 1
        1186      983520 :           x    (:,2,iarea) = xtmp(:)                     !dynamic -line 2
        1187      983520 :           dx   (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1188             :        case default
        1189 57430365546 :           call endrun('flow case not defined in get_flux_segments_area_iterate')
        1190             :        end select
        1191             :        !
        1192             :        ! compute flux integral
        1193             :        !
        1194 >23620*10^7 :        f2 = 0.0_r8
        1195 >23620*10^7 :        weight_area=w_static
        1196 >23620*10^7 :        do iarea=1,num_area
        1197 >34517*10^7 :          do iseg=1,num_seg(iarea)
        1198             : !rck vector directive needed here
        1199             : !DIR$ SIMD
        1200 >66558*10^7 :            do ipt=1,ngpc
        1201             :
        1202 >49918*10^7 :              xq(ipt) = x(1,iseg,iarea)+dx(1,iseg,iarea)*gspts(ipt)! create quadrature point locations
        1203 >49918*10^7 :              yq(ipt) = x(2,iseg,iarea)+dx(2,iseg,iarea)*gspts(ipt)
        1204             :
        1205 >49918*10^7 :              xq2      =  xq(ipt)*xq(ipt)
        1206 >49918*10^7 :              xq2i     =  1.0_r8/(1.0_r8+xq2)
        1207 >49918*10^7 :              rho      =  SQRT(1.0_r8+xq2+yq(ipt)*yq(ipt))
        1208 >49918*10^7 :              rhoi     =  1.0_r8/rho
        1209 >49918*10^7 :              yrh      =  yq(ipt)*rhoi
        1210 >66558*10^7 :              F(ipt,1) =  yrh*xq2i
        1211             :            end do
        1212 >84436*10^7 :            weight_area(iarea) = weight_area(iarea)+sum(gsweights(:)*F(:,1))*0.5_r8*dx(1,iseg,iarea)! integral
        1213             :          end do
        1214 >23620*10^7 :          f2 = f2+weight_area(iarea)*c(iarea)
        1215             :        end do
        1216 57430365546 :        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 77721254328 :        if (ABS(f2)<eps.or.lexit_after_one_more_iteration) then
        1224 20303098026 :          gamma=gamma3
        1225 20303098026 :          if (gamma>gamma_max) then
        1226    12209244 :            lexit_after_one_more_iteration=.true.
        1227    12209244 :            gamma=gamma_max
        1228    12209244 :            gamma3=gamma_max
        1229             :          else
        1230             :            exit
        1231             :          end if
        1232             :        else
        1233             :           !
        1234             :           ! Newton increment
        1235             :           !
        1236 37127267520 :          if (abs(f2-f1)<eps) then
        1237             :            !
        1238             :            ! if entering here abs(f2)>eps and abs(f1)>eps but abs(f2-f1)<eps
        1239             :            !
        1240          79 :            dgamma=-0.5_r8*(gamma2-gamma1)
        1241          79 :            lexit_after_one_more_iteration=.true.
        1242             :          else
        1243 37127267441 :            dgamma=(gamma2-gamma1)*f2/(f2-f1)
        1244             :          endif
        1245 37127267520 :          if (ABS(dgamma)>eps) then
        1246 37127267520 :            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 37127267520 :          gamma3=MAX(gamma3,gamma_min)
        1254             :          !
        1255             :          ! prepare for next iteration
        1256             :          !
        1257 37127267520 :          gamma1 = gamma2; f1 = f2; gamma2 = gamma3;
        1258             :        endif
        1259             :      end do
        1260 20290888782 :      if (iter>iter_max) write(iulog,*) "WARNING: iteration not converged",&
        1261           0 :           ABS(f2),flux,gamma1,gamma2,gamma3,ilev
        1262             :   end subroutine get_flux_segments_area_iterate
        1263             :
        1264   483116400 :   subroutine define_swept_areas(fvm,ilev,displ,base_vec,base_vtx,idx)
        1265             :     use control_mod, only : neast, nwest, seast, swest
        1266             :     implicit none
        1267             :     type (fvm_struct), intent(inout) :: fvm
        1268             :     integer          , intent(in)    :: ilev
        1269             :
        1270             :
        1271             :     integer, parameter :: num_area=5, num_sides=4, imin= 0, imax=nc+1
        1272             :     real (kind=r8)    , dimension(0:7       , imin:imax,imin:imax,num_sides), intent(out) :: displ
        1273             :     integer (kind=r8) , dimension(1:2,11    , imin:imax,imin:imax,num_sides), intent(out) :: base_vec
        1274             :     real (kind=r8)    , dimension(1:2, 6    , imin:imax,imin:imax,num_sides), intent(out) :: base_vtx
        1275             :     integer           , dimension(2,num_area, imin:imax,imin:imax,num_sides), intent(out) :: idx
        1276             :
        1277             :     real (kind=r8)        :: flux_sum     (0:nc+1,0:nc+1,2)
        1278             :     integer               :: degenerate   (1:nc+1,1:nc+1  )
        1279             :     integer               :: circular_flow(1:nc+1,1:nc+1  )
        1280             :     integer               :: illcond      (1:nc+1,1:nc+1)
        1281             :     integer               :: ib,i,j,sgn, iside, iarea
        1282             :
        1283             :     !
        1284             :     ! set where reconstruction function is as a function of area and side
        1285             :     !
        1286             :     integer, dimension(num_area*4), parameter :: idx_shift_tmp = (/-1,-1, 0, 1, 1,&  !iside=1
        1287             :                                                                     1, 0, 0, 0, 1,&  !iside=2
        1288             :                                                                     1, 1, 0,-1,-1,&  !iside=3
        1289             :                                                                    -1, 0, 0, 0,-1/)  !iside=4
        1290             :
        1291             :     integer, dimension(num_area*4), parameter :: idy_shift_tmp = (/-1, 0, 0, 0,-1,&  !iside=1
        1292             :                                                                    -1,-1, 0, 1, 1,&  !iside=2
        1293             :                                                                     1, 0, 0, 0, 1,&  !iside=3
        1294             :                                                                     1, 1, 0,-1,-1/)  !iside=4
        1295             :
        1296             :     integer, dimension(num_area,4), parameter :: idx_shift = RESHAPE(idx_shift_tmp,(/num_area,4/))
        1297             :     integer, dimension(num_area,4), parameter :: idy_shift = RESHAPE(idy_shift_tmp,(/num_area,4/))
        1298             :
        1299             :     integer, dimension(4), parameter :: iside_m1 = (/4,1,2,3/)
        1300             :     integer, dimension(4), parameter :: iside_p1 = (/2,3,4,1/)
        1301             :     integer, dimension(4), parameter :: iside_p2 = (/3,4,1,2/)
        1302             :     integer, dimension(4), parameter :: iside_p3 = (/4,1,2,3/)
        1303             :
        1304             :     integer, dimension(4), parameter :: imin_side = (/1   ,0   ,1   ,1   /)
        1305             :     integer, dimension(4), parameter :: imax_side = (/nc  ,nc  ,nc  ,nc+1/)
        1306             :     integer, dimension(4), parameter :: jmin_side = (/1   ,1   ,0   ,1   /)
        1307             :     integer, dimension(4), parameter :: jmax_side = (/nc+1,nc  ,nc  ,nc  /)
        1308             :
        1309             :
        1310             :
        1311             :     integer :: iur,jur,ilr,jlr,iul,jul,ill,jll
        1312             :
        1313   483116400 :     ib = fvm%cubeboundary
        1314 12077910000 :     flux_sum(0:nc+1,1:nc+1,1) = fvm%se_flux(0:nc+1,0:nc  ,3,ilev)-fvm%se_flux(0:nc+1,1:nc+1,1,ilev)
        1315 12561026400 :     flux_sum(1:nc+1,0:nc+1,2) = fvm%se_flux(0:nc  ,0:nc+1,2,ilev)-fvm%se_flux(1:nc+1,0:nc+1,4,ilev)
        1316             :
        1317             :     !
        1318             :     ! Degenerate case ("two departure points")
        1319             :     !
        1320             :     !           ||  |                        || no change in this situation ||  no change in this situation
        1321             :     !           ||  |                        ||                             ||
        1322             :     !           ||--------                   ||----------                   ||----------
        1323             :     !           ||  |                        ||                             ||
        1324             :     ! =======================      =======================         =====================
        1325             :     !       |   ||                       |   ||                             ||
        1326             :     !  -----|---||                 ------|---||                    ---------||
        1327             :     !       |   ||                       |   ||                             ||
        1328             :     !       |   ||                       |   ||                             ||
        1329             :     !
        1330             :     !
        1331 10145444400 :     where (flux_sum(0:nc,1:nc+1,1)*flux_sum(1:nc+1,1:nc+1,1)<0.0_r8.and.flux_sum(1:nc+1,0:nc,2)*flux_sum(1:nc+1,1:nc+1,2)<0.0_r8)
        1332             :        degenerate(:,:) = 0
        1333             :     elsewhere
        1334             :        degenerate(:,:) = 1
        1335             :     end where
        1336             :
        1337   483116400 :     if (ib>0) then
        1338    62268336 :        if (ib==swest) degenerate(1   ,1   ) = 1
        1339    61731540 :        if (ib==nwest) degenerate(1   ,nc+1) = 1
        1340    61731540 :        if (ib==neast) degenerate(nc+1,nc+1) = 1
        1341    61731540 :        if (ib==seast) degenerate(nc+1,1   ) = 1
        1342             :     end if
        1343             :
        1344  2415582000 :     do j=1,nc+1
        1345 10145444400 :        do i=1,nc+1
        1346 17392190400 :           do sgn=-1,1,2
        1347             :              if (&
        1348 61838899200 :                   sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
        1349 85028486400 :                   sgn*flux_sum(i  ,j,1)>0.0_r8.and.sgn*flux_sum(i,j  ,2)<0.0_r8) then
        1350    12213178 :                 circular_flow(i,j) = 0
        1351             :              else
        1352 15447511622 :                 circular_flow(i,j) = 1
        1353             :              end if
        1354             :           end do
        1355             :        end do
        1356             :     end do
        1357             :     !
        1358             :     ! wrap around corners
        1359             :     !
        1360   483116400 :     if (ib==nwest) then
        1361      536796 :        flux_sum(0,nc+1,1) = fvm%se_flux(0,nc,3,ilev)-fvm%se_flux(1,nc+1,4,ilev)
        1362      536796 :        flux_sum(1,nc+1,2) = fvm%se_flux(0,nc,3,ilev)-fvm%se_flux(1,nc+1,4,ilev)
        1363             :
        1364      536796 :        i=1;j=nc+1;
        1365      536796 :        circular_flow(i,j) = 1
        1366     1610388 :        do sgn=-1,1,2
        1367             :           if (&
        1368             :                sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
        1369     1610388 :                sgn*flux_sum(i  ,j,1)>0.0_r8.and.sgn*flux_sum(i,j  ,2)<0.0_r8) then
        1370        2579 :              circular_flow(i,j) = 0
        1371             :           end if
        1372             :        end do
        1373   482579604 :     else if (ib==swest) then
        1374      536796 :        flux_sum(0,1,1) = fvm%se_flux(1,0,4,ilev)-fvm%se_flux(0,1,1,ilev)
        1375      536796 :        flux_sum(1,0,2) = fvm%se_flux(0,1,1,ilev)-fvm%se_flux(1,0,4,ilev)
        1376      536796 :        i=1;j=1;
        1377      536796 :        circular_flow(i,j) = 1
        1378     1610388 :        do sgn=-1,1,2
        1379             :           if (&
        1380             :                sgn*flux_sum(i-1,j,1)<0.0_r8.and.&
        1381     1610388 :                sgn*flux_sum(i  ,j,1)>0.0_r8.and.sgn*flux_sum(i,j  ,2)<0.0_r8) then
        1382        1804 :              circular_flow(i,j) = 0
        1383             :           end if
        1384             :        end do
        1385   482042808 :     else if (ib==neast) then
        1386      536796 :        flux_sum(nc+1,nc+1,1) = fvm%se_flux(nc+1,nc,3,ilev)-fvm%se_flux(nc,nc+1,2,ilev)
        1387      536796 :        flux_sum(nc+1,nc+1,2) = fvm%se_flux(nc,nc+1,2,ilev)-fvm%se_flux(nc+1,nc,3,ilev)
        1388      536796 :        i=nc+1;j=nc+1;
        1389      536796 :        circular_flow(i,j) = 1
        1390     1610388 :        do sgn=-1,1,2
        1391             :           if (&
        1392     1073592 :                sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
        1393      536796 :                sgn*flux_sum(i,j  ,2)<0.0_r8) then
        1394        1303 :              circular_flow(i,j) = 0
        1395             :           end if
        1396             :        end do
        1397   481506012 :     else if (ib==seast) then
        1398      536796 :        flux_sum(nc+1,1   ,1) = fvm%se_flux(nc,0,2,ilev)-fvm%se_flux(nc+1,1,1,ilev)
        1399      536796 :        flux_sum(nc+1,0   ,2) = fvm%se_flux(nc,0,2,ilev)-fvm%se_flux(nc+1,1,1,ilev)
        1400      536796 :        i=nc+1;j=1;
        1401      536796 :        circular_flow(i,j) = 1
        1402     1610388 :        do sgn=-1,1,2
        1403             :           if (&
        1404     1073592 :                sgn*flux_sum(i-1,j,1)<0.0_r8.and.sgn*flux_sum(i,j-1,2)>0.0_r8.and.&
        1405      536796 :                sgn*flux_sum(i,j  ,2)<0.0_r8) then
        1406        1835 :              circular_flow(i,j) = 0
        1407             :           end if
        1408             :        end do
        1409             :     end if
        1410 10145444400 :     illcond = circular_flow*degenerate
        1411             :     !
        1412             :     !
        1413             :     !
        1414             :     !
        1415  2415582000 :     do iside=1,4
        1416  9179211600 :        do j=jmin_side(iside),jmax_side(iside)
        1417 31885682400 :           do i=imin_side(iside),imax_side(iside)
        1418 29953216800 :              if (fvm%se_flux(i,j,iside,ilev)>eps) then
        1419 11594793589 :                 iur = i+idx_shift(4,iside); jur = j+idy_shift(4,iside) !(i,j) index of upper right quadrant
        1420 11594793589 :                 ilr = i+idx_shift(5,iside); jlr = j+idy_shift(5,iside) !(i,j) index of lower left  quadrant
        1421 11594793589 :                 iul = i+idx_shift(2,iside); jul = j+idy_shift(2,iside) !(i,j) index of upper right quadrant
        1422 11594793589 :                 ill = i+idx_shift(1,iside); jll = j+idy_shift(1,iside) !(i,j) index of lower left  quadrant
        1423             :
        1424             :                 !iside=1
        1425 11594793589 :                 if (iside==1) then
        1426  2936257703 :                 displ(0,i,j,iside) = -flux_sum   (i  ,j  ,1)*illcond(i,j)     !center left
        1427  2936257703 :                 displ(1,i,j,iside) = -flux_sum   (i  ,j  ,1)*illcond(i+1,j)   !center right
        1428  2936257703 :                 displ(2,i,j,iside) =  flux_sum   (i+1,j  ,2)*illcond(i+1,j)   !c2
        1429  2936257703 :                 displ(3,i,j,iside) = -flux_sum   (i  ,j  ,2)*illcond(i  ,j)   !c3
        1430  2936257703 :                 displ(4,i,j,iside) = -flux_sum   (i+1,j  ,1)*illcond(i+1,j)   !r1
        1431  2936257703 :                 displ(5,i,j,iside) = -flux_sum   (i+1,j-1,2)*illcond(i+1,j)   !r2
        1432  2936257703 :                 displ(6,i,j,iside) = -flux_sum   (i-1,j  ,1)*illcond(i  ,j)   !l1
        1433  2936257703 :                 displ(7,i,j,iside) =  flux_sum   (i  ,j-1,2)*illcond(i  ,j)   !l2
        1434             :
        1435             :                 end if
        1436 11594793589 :                 if (iside==2) then
        1437             :                 !iside=2
        1438  2852956151 :                 displ(0,i,j,iside) =  flux_sum   (i+1,j  ,2)*illcond(i+1,j  )     !center left
        1439  2852956151 :                 displ(1,i,j,iside) =  flux_sum   (i+1,j  ,2)*illcond(i+1,j+1)   !center right
        1440  2852956151 :                 displ(2,i,j,iside) =  flux_sum   (i  ,j+1,1)*illcond(i+1,j+1)   !c2
        1441  2852956151 :                 displ(3,i,j,iside) = -flux_sum   (i  ,j  ,1)*illcond(i+1,j  )   !c3
        1442  2852956151 :                 displ(4,i,j,iside) =  flux_sum   (i+1,j+1,2)*illcond(i+1,j+1)   !r1
        1443  2852956151 :                 displ(5,i,j,iside) = -flux_sum   (i+1,j+1,1)*illcond(i+1,j+1)   !r2
        1444  2852956151 :                 displ(6,i,j,iside) =  flux_sum   (i+1,j-1,2)*illcond(i+1,j)   !l1
        1445  2852956151 :                 displ(7,i,j,iside) =  flux_sum   (i+1,j  ,1)*illcond(i+1,j)   !l2
        1446             :                 end if
        1447             :                 !iside=3
        1448 11594793589 :                 if (iside==3) then
        1449  2861139088 :                 displ(0,i,j,iside) =  flux_sum   (i  ,j+1,1)*illcond(i+1,j+1)     !center left
        1450  2861139088 :                 displ(1,i,j,iside) =  flux_sum   (i  ,j+1,1)*illcond(i  ,j+1)   !center right
        1451  2861139088 :                 displ(2,i,j,iside) = -flux_sum   (i  ,j  ,2)*illcond(i  ,j+1)   !c2
        1452  2861139088 :                 displ(3,i,j,iside) =  flux_sum   (i+1,j  ,2)*illcond(i+1,j+1)   !c3
        1453  2861139088 :                 displ(4,i,j,iside) =  flux_sum   (i-1,j+1,1)*illcond(i  ,j+1)   !r1
        1454  2861139088 :                 displ(5,i,j,iside) =  flux_sum   (i  ,j+1,2)*illcond(i  ,j+1)   !r2
        1455  2861139088 :                 displ(6,i,j,iside) =  flux_sum   (i+1,j+1,1)*illcond(i+1,j+1)   !l1
        1456  2861139088 :                 displ(7,i,j,iside) = -flux_sum   (i+1,j+1,2)*illcond(i+1,j+1)   !l2
        1457             :                 end if
        1458 11594793589 :                 if (iside==4) then
        1459             :                 !iside=4
        1460  2944440647 :                 displ(0,i,j,iside) = -flux_sum   (i  ,j  ,2)*illcond(i  ,j+1)     !center left
        1461  2944440647 :                 displ(1,i,j,iside) = -flux_sum   (i  ,j  ,2)*illcond(i  ,j  )   !center right
        1462  2944440647 :                 displ(2,i,j,iside) = -flux_sum   (i  ,j  ,1)*illcond(i  ,j  )   !c2
        1463  2944440647 :                 displ(3,i,j,iside) =  flux_sum   (i  ,j+1,1)*illcond(i  ,j+1)   !c3
        1464  2944440647 :                 displ(4,i,j,iside) = -flux_sum   (i  ,j-1,2)*illcond(i  ,j  )   !r1
        1465  2944440647 :                 displ(5,i,j,iside) =  flux_sum   (i-1,j  ,1)*illcond(i  ,j  )   !r2
        1466  2944440647 :                 displ(6,i,j,iside) = -flux_sum   (i  ,j+1,2)*illcond(i  ,j+1)   !l1
        1467  2944440647 :                 displ(7,i,j,iside) = -flux_sum   (i-1,j+1,1)*illcond(i  ,j+1)   !l2
        1468             :                 end if
        1469             :
        1470 34784380767 :                 base_vtx(:,1,i,j,iside) = fvm%vtx_cart(iside,:,i  ,j            )       !vertex center left
        1471 34784380767 :                 base_vtx(:,2,i,j,iside) = fvm%vtx_cart(iside_p1(iside),:,i  ,j  )       !vertex center right
        1472 34784380767 :                 base_vtx(:,3,i,j,iside) = fvm%vtx_cart(iside,:,iur,jur          )       !vertex upper right
        1473 34784380767 :                 base_vtx(:,4,i,j,iside) = fvm%vtx_cart(iside_p3(iside),:,ilr,jlr)       !vertex lower right
        1474 34784380767 :                 base_vtx(:,5,i,j,iside) = fvm%vtx_cart(iside_p1(iside),:,iul,jul)       !vertex upper left
        1475 34784380767 :                 base_vtx(:,6,i,j,iside) = fvm%vtx_cart(iside_p2(iside),:,ill,jll)       !vertex lower left
        1476             :
        1477 34784380767 :                 base_vec(:, 1,i,j,iside) = fvm%flux_vec    (:,i  ,j  ,iside          )      !vector center
        1478 34784380767 :                 base_vec(:, 2,i,j,iside) = fvm%flux_vec    (:,i  ,j  ,iside_p1(iside))      !vector center right
        1479 34784380767 :                 base_vec(:, 3,i,j,iside) = fvm%flux_vec    (:,i  ,j  ,iside_p3(iside))      !vector center left
        1480 34784380767 :                 base_vec(:, 4,i,j,iside) = fvm%flux_vec    (:,iur,jur,iside          )      !vector upper right 1
        1481 34784380767 :                 base_vec(:, 5,i,j,iside) = fvm%flux_vec    (:,iur,jur,iside_p3(iside))      !vector upper right 2
        1482 34784380767 :                 base_vec(:, 6,i,j,iside) = fvm%flux_vec    (:,ilr,jlr,iside_p3(iside))      !vector lower right 1
        1483 34784380767 :                 base_vec(:, 7,i,j,iside) = fvm%flux_vec    (:,ilr,jlr,iside_p2(iside))      !vector lower right 2
        1484 34784380767 :                 base_vec(:, 8,i,j,iside) = fvm%flux_vec    (:,iul,jul,iside          )      !vector upper left 1
        1485 34784380767 :                 base_vec(:, 9,i,j,iside) = fvm%flux_vec    (:,iul,jul,iside_p1(iside))      !vector upper left 2
        1486 34784380767 :                 base_vec(:,10,i,j,iside) = fvm%flux_vec    (:,ill,jll,iside_p1(iside))      !vector lower left 1
        1487 34784380767 :                 base_vec(:,11,i,j,iside) = fvm%flux_vec    (:,ill,jll,iside_p2(iside))      !vector lower left 2
        1488             :
        1489 69568761534 :                 do iarea=1,5
        1490 57973967945 :                    idx(1,iarea,i,j,iside) = i+idx_shift(iarea,iside)
        1491 69568761534 :                    idx(2,iarea,i,j,iside) = j+idy_shift(iarea,iside)
        1492             :                 end do
        1493             :              else
        1494 >10435*10^7 :                 displ(:,i,j,iside) = 9D99!for debugging
        1495             :              end if
        1496             :           end do
        1497             :        end do
        1498             :     end do
        1499             :     !
        1500             :     ! wrap around corners here
        1501             :     !
        1502             :
        1503   483116400 :   end subroutine define_swept_areas
        1504             :
        1505             :
        1506             :   !
        1507             :   ! Notation conventions used in define_area subroutines
        1508             :   !
        1509             :   !
        1510             :   !
        1511             :   !   ^    ||--->   ^   <---||    ^
        1512             :   !  /|\   || 3    /|\    2 ||   /|\
        1513             :   !   | 6  ||     1 |       ||    | 4
        1514             :   !   |    ||       |       ||    |
        1515             :   ! =================================
        1516             :   !        ||               ||
        1517             :   !        ||               ||
        1518             :   !      7 ||               || 5
        1519             :   !    <---||               ||--->
        1520             :   !
        1521             :
        1522   224442892 :   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   673328676 :     if (SUM(ABS(base_vec(:,9,i,j,iside))).NE.0) then
        1555   223985704 :        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   223985704 :        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      457188 :        gamma=displ(0,i,j,iside)
        1563             :     end if
        1564             :
        1565             :
        1566   673328676 :     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   673328676 :     x_start (:,1) = base_vtx(:, 6,i,j,iside)
        1568   673328676 :     dgam_vec(:,1) = base_vec(:,10,i,j,iside)*gamma
        1569             :
        1570   673328676 :     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   224442892 :     iarea                  = 1
        1573   224442892 :     num_seg       (iarea)  = 2
        1574   224442892 :     num_seg_static(iarea)  = 1
        1575             :
        1576   673328676 :     x_static (:,1,iarea) = base_vtx(:,6,i,j,iside)       !static
        1577   673328676 :     dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static
        1578             :
        1579   673328676 :     xtmp(:        ) = x_start(:,1)+dgam_vec(:,1)
        1580   673328676 :     x   (:,1,iarea) = xdep(:,1)                  !static
        1581   673328676 :     dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1582             :
        1583   673328676 :     x (:,2,iarea) = xtmp(:)                      !dynamic
        1584   673328676 :     dx(:,2,iarea) = x_static(:,1,iarea)-xtmp(:)  !dynamic
        1585             :     !
        1586             :     !
        1587             :     !
        1588   224442892 :     iarea                  = 2
        1589   224442892 :     num_seg       (iarea)  = 3
        1590             :
        1591   673328676 :     x_start (:,2) = base_vtx(:,5,i,j,iside)
        1592   673328676 :     dgam_vec(:,2) = base_vec(:,9,i,j,iside)*gamma
        1593   673328676 :     xtmp    (:  ) = x_start(:,2)+dgam_vec(:,2)
        1594             :
        1595   673328676 :     x_start (:,3) = base_vtx(:,5,i,j,iside)
        1596   673328676 :     dgam_vec(:,3) = base_vec(:,8,i,j,iside)*displ(0,i,j,iside)
        1597   673328676 :     xtmp2   (:  ) = x_start(:,3)+dgam_vec(:,3)
        1598             :
        1599   673328676 :     x   (:,1,iarea) = base_vtx(:,5,i,j,iside) !static
        1600   673328676 :     dx  (:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic
        1601             :
        1602   673328676 :     x (:,2,iarea) = xtmp (:)          !dynamic
        1603   673328676 :     dx(:,2,iarea) = xtmp2(:)-xtmp(:)  !dynamic
        1604             :
        1605   673328676 :     x (:,3,iarea) = xtmp2(:)                !dynamic
        1606   673328676 :     dx(:,3,iarea) = x(:,1,iarea)-xtmp2(:)   !dynamic
        1607   224442892 :   end subroutine define_area1_area2
        1608             :
        1609             :
        1610  5560574755 :   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 16681724265 :          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 16681724265 :     x_start (:,1) = base_vtx(:,5,i,j,iside)
        1647  5560574755 :     gamma         = displ(0,i,j,iside)
        1648 16681724265 :     dgam_vec(:,1) = base_vec(:,8,i,j,iside)*gamma
        1649  5560574755 :     if (present(displ_first_guess)) displ_first_guess = gamma
        1650             :
        1651  5560574755 :     iarea                  = 2
        1652  5560574755 :     num_seg       (iarea)  = 2
        1653  5560574755 :     num_seg_static(iarea)  = 1
        1654             :
        1655 16681724265 :     x_static (:,1,iarea) = base_vtx(:,5,i,j,iside)       !static  - line 1
        1656 16681724265 :     dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static  - line 1
        1657             :
        1658 16681724265 :     xtmp     (:        ) = x_start(:,1)+dgam_vec(:,1)
        1659 16681724265 :     x        (:,1,iarea) = xdep(:,1)                  !static  - line 2
        1660 16681724265 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic - line 2
        1661             :
        1662 16681724265 :     x        (:,2,iarea) = xtmp(:)                     !dynamic  - line 3
        1663 16681724265 :     dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 3
        1664  5560574755 :   end subroutine define_area2
        1665             :
        1666             :
        1667  5598169817 :   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  5598169817 :     xdep(:,2) = base_vtx(:,2,i,j,iside)+displ(1,i,j,iside)*base_vec(:,1,i,j,iside)&
        1700 22392679268 :          +MAX(0.0_r8,displ(2,i,j,iside))*base_vec(:,2,i,j,iside)
        1701 16794509451 :     x_start (:,4) = base_vtx(:,1,i,j,iside)
        1702  5598169817 :     gamma         = displ(0,i,j,iside)
        1703 16794509451 :     dgam_vec(:,4) = base_vec(:,1,i,j,iside)*gamma
        1704 16794509451 :     xtmp    (:  ) = x_start(:,4)+dgam_vec(:,4)
        1705             :
        1706  5598169817 :     if (present(displ_first_guess)) displ_first_guess = gamma
        1707             :
        1708  5598169817 :     iarea                  = 3
        1709  5598169817 :     num_seg       (iarea)  = 2
        1710  5598169817 :     num_seg_static(iarea)  = 2
        1711             :
        1712 16794509451 :     x_static (:,1,iarea) = xdep(:,2)                         !static  - line 3
        1713 16794509451 :     dx_static(:,1,iarea) = base_vtx(:,2,i,j,iside)-xdep(:,2) !static  - line 3
        1714             :
        1715 16794509451 :     x_static (:,2,iarea) = base_vtx(:,2,i,j,iside)                         !static  - line 4
        1716 16794509451 :     dx_static(:,2,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static  - line 4
        1717             :
        1718 16794509451 :     x        (:,1,iarea) = base_vtx(:,1,i,j,iside)    !static  - line 1
        1719 16794509451 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic - line 1
        1720             :
        1721 16794509451 :     x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1722 16794509451 :     dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1723  5598169817 :   end subroutine define_area3_left
        1724             :
        1725  5596792086 :   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  5596792086 :     xdep(:,1) = base_vtx(:,1,i,j,iside)+displ(0,i,j,iside)*base_vec(:,1,i,j,iside)&
        1755 22387168344 :          +MAX(0.0_r8,displ(3,i,j,iside))*base_vec(:,3,i,j,iside)
        1756 16790376258 :     x_start (:,5) = base_vtx(:,2,i,j,iside)
        1757  5596792086 :     gamma         = displ(1,i,j,iside)
        1758 16790376258 :     dgam_vec(:,5) = base_vec(:,1,i,j,iside)*gamma
        1759 16790376258 :     xtmp    (:  ) = x_start(:,5)+dgam_vec(:,5)
        1760             :
        1761  5596792086 :     iarea                  = 3
        1762  5596792086 :     num_seg       (iarea)  = 2
        1763  5596792086 :     num_seg_static(iarea)  = 2
        1764             :
        1765 16790376258 :     x_static (:,1,iarea) = base_vtx(:,1,i,j,iside)           !static  - line 1
        1766 16790376258 :     dx_static(:,1,iarea) = xdep(:,1)-base_vtx(:,1,i,j,iside) !static  - line 1
        1767             :
        1768 16790376258 :     x_static (:,2,iarea) = base_vtx(:,2,i,j,iside)                         !static  - line 4
        1769 16790376258 :     dx_static(:,2,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static  - line 4
        1770             :
        1771 16790376258 :     x        (:,1,iarea) = xdep(:,1)            !static  - line 2
        1772 16790376258 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea) !dynamic - line 2
        1773             :
        1774 16790376258 :     x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1775 16790376258 :     dx       (:,2,iarea) = x_static(:,2,iarea)-xtmp(:) !dynamic - line 2
        1776  5596792086 :   end subroutine define_area3_right
        1777             :
        1778             :
        1779   186847830 :   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   560543490 :     x_start (:,4) = base_vtx(:,1,i,j,iside)
        1805   560543490 :     x_start (:,5) = base_vtx(:,2,i,j,iside)
        1806   186847830 :     gamma         = displ(0,i,j,iside)
        1807   560543490 :     dgam_vec(:,4) = base_vec(:,1,i,j,iside)*gamma
        1808   560543490 :     dgam_vec(:,5) = base_vec(:,1,i,j,iside)*gamma
        1809   560543490 :     xtmp    (:  ) = x_start(:,4)+dgam_vec(:,4)
        1810   560543490 :     xtmp2   (:  ) = x_start(:,5)+dgam_vec(:,5)
        1811             :
        1812   186847830 :     iarea                  = 3
        1813   186847830 :     num_seg       (iarea)  = 3
        1814   186847830 :     num_seg_static(iarea)  = 1
        1815             :
        1816   560543490 :     x_static (:,1,iarea) = base_vtx(:,2,i,j,iside)                         !static
        1817   560543490 :     dx_static(:,1,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static
        1818             :
        1819   560543490 :     x        (:,1,iarea) = base_vtx(:,1,i,j,iside)    !static
        1820   560543490 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)       !dynamic
        1821             :
        1822   560543490 :     x        (:,2,iarea) = xtmp (:)         !dynamic
        1823   560543490 :     dx       (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic
        1824             :
        1825   560543490 :     x        (:,3,iarea) = xtmp2(:)              !dynamic
        1826   560543490 :     dx       (:,3,iarea) = x_start(:,5)-xtmp2(:) !dynamic
        1827   186847830 :   end subroutine define_area3_left_right
        1828             :
        1829   225596637 :   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   225596637 :     iarea                  = 4
        1862   225596637 :     num_seg       (iarea)  = 3
        1863             :
        1864   676789911 :     if (SUM(ABS(base_vec(:,5,i,j,iside))).NE.0) then
        1865   225340604 :        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   225340604 :        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      256033 :        gamma = displ(1,i,j,iside)
        1873             :     end if
        1874             :
        1875   225596637 :     if (present(displ_first_guess)) displ_first_guess = displ(1,i,j,iside)
        1876             :
        1877   676789911 :     x_start (:,6) = base_vtx(:,3,i,j,iside)
        1878   676789911 :     dgam_vec(:,6) = base_vec(:,4,i,j,iside)*displ(1,i,j,iside)
        1879   676789911 :     xtmp    (:  ) = x_start(:,6)+dgam_vec(:,6)
        1880   676789911 :     x_start (:,7) = base_vtx(:,3,i,j,iside)
        1881   676789911 :     dgam_vec(:,7) = base_vec(:,5,i,j,iside)*gamma
        1882   676789911 :     xtmp2   (:  ) = x_start(:,7)+dgam_vec(:,7)
        1883             :
        1884   676789911 :     x        (:,1,iarea) = base_vtx(:,3,i,j,iside)!static   -line 1
        1885   676789911 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)   !dynamic - line 1
        1886             :
        1887   676789911 :     x        (:,2,iarea) = xtmp(:)          !dynamic -line 2
        1888   676789911 :     dx       (:,2,iarea) = xtmp2(:)-xtmp(:) !dynamic - line 2
        1889             :
        1890   676789911 :     x        (:,3,iarea) = xtmp2(:)               !static   -line 1
        1891   676789911 :     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   676789911 :          -displ(4,i,j,iside)*base_vec(:,7,i,j,iside)
        1897   676789911 :     x_start (:,8) = base_vtx(:,4,i,j,iside)
        1898   676789911 :     dgam_vec(:,8) = base_vec(:,6,i,j,iside)*gamma
        1899   676789911 :     xtmp    (:  ) = x_start(:,8)+dgam_vec(:,8)
        1900             :
        1901   225596637 :     iarea                  = 5
        1902   225596637 :     num_seg       (iarea)  = 2
        1903   225596637 :     num_seg_static(iarea)  = 1
        1904             :
        1905   676789911 :     x        (:,1,iarea) = base_vtx(:,4,i,j,iside)!static   -line 1
        1906   676789911 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)   !dynamic - line 1
        1907             :
        1908   676789911 :     x_static (:,1,iarea) = xdep(:,1)                        !static - line 1
        1909   676789911 :     dx_static(:,1,iarea) = x(:,1,iarea)-x_static(:,1,iarea) !static - line 1
        1910             :
        1911   676789911 :     x        (:,2,iarea) = xtmp(:)                     !dynamic -line 2
        1912   676789911 :     dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1913   225596637 :   end subroutine define_area4_area5
        1914             :
        1915             :
        1916  5558043279 :   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  5558043279 :     iarea                  = 4
        1939  5558043279 :     num_seg       (iarea)  = 2
        1940  5558043279 :     num_seg_static(iarea)  = 1
        1941             :
        1942  5558043279 :     xdep(:,1) = base_vtx(:,3,i,j,iside)+MAX(0.0_r8,displ(4,i,j,iside))*base_vec(:,4,i,j,iside)&
        1943 22232173116 :          -displ(2,i,j,iside)*base_vec(:,5,i,j,iside)
        1944 16674129837 :     x_start (:,6) = base_vtx(:,3,i,j,iside)
        1945  5558043279 :     gamma         = displ(1,i,j,iside)
        1946 16674129837 :     dgam_vec(:,6) = base_vec(:,4,i,j,iside)*gamma
        1947 16674129837 :     xtmp    (:  ) = x_start(:,6)+dgam_vec(:,6)
        1948             :
        1949  5558043279 :     if (present(displ_first_guess)) displ_first_guess = gamma
        1950             :
        1951 16674129837 :     x_static (:,1,iarea) = xdep(:,1)                         !static
        1952 16674129837 :     dx_static(:,1,iarea) = base_vtx(:,3,i,j,iside)-xdep(:,1) !static
        1953             :
        1954 16674129837 :     x        (:,1,iarea) = base_vtx(:,3,i,j,iside) !static  - line 2
        1955 16674129837 :     dx       (:,1,iarea) = xtmp(:)-x(:,1,iarea)    !dynamic - line 2
        1956             :
        1957 16674129837 :     x        (:,2,iarea) = xtmp(:)                     !dynamic  -line 2
        1958 16674129837 :     dx       (:,2,iarea) = x_static(:,1,iarea)-xtmp(:) !dynamic - line 2
        1959  5558043279 :   end subroutine define_area4
        1960             :
        1961   212983856 :   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   638951568 :          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   638951568 :          displ(1,i,j,iside)*base_vec(:,1,i,j,iside)+displ(2,i,j,iside)*base_vec(:,2,i,j,iside)
        1998   638951568 :     xdep(:,2) = 0.5_r8*(xdep(:,1)+xdep(:,3))
        1999             :
        2000   212983856 :     gamma= se_flux_center
        2001             :     x_start(:,1) = ABS(base_vec(:,3,i,j,iside))*((xdep(:,2)-base_vtx(:,1,i,j,iside)))+&
        2002   638951568 :          base_vtx(:,1,i,j,iside) !xdep(2) - midway between departure points projected to side 1
        2003             :
        2004   638951568 :     dgam_vec(:,1) = gamma*base_vec(:,1,i,j,iside)
        2005             :
        2006   212983856 :     if (present(displ_first_guess)) displ_first_guess = gamma
        2007             :
        2008   638951568 :     xdep(:,2)     = x_start(:,1)+dgam_vec(:,1)
        2009   212983856 :     iarea                  = 3
        2010   212983856 :     num_seg       (iarea)  = 2
        2011   212983856 :     num_seg_static(iarea)  = 3
        2012             :
        2013             :     !                 ______X______
        2014             :     !        ||    2 /             \ 3    ||
        2015             :     !        ||  *--/               \--*  ||
        2016             :     !        || /                       \ ||
        2017             :     !        ||/ 1          5           4\||
        2018             :     !  ========================================
        2019             :     !        ||                           ||
        2020             :     !
        2021   638951568 :     x_static (:,1,iarea) = base_vtx(:,1,i,j,iside)       !static  - line 1
        2022   638951568 :     dx_static(:,1,iarea) = xdep(:,1)-x_static(:,1,iarea) !static  - line 1
        2023             :
        2024   638951568 :     x        (:,1,iarea) = xdep(:,1)                     !static  - line 2
        2025   638951568 :     dx       (:,1,iarea) = xdep(:,2)-x(:,1,iarea)        !dynamic - line 2
        2026             :
        2027   638951568 :     x        (:,2,iarea) = xdep(:,2)                     !dynamic - line 3
        2028   638951568 :     dx       (:,2,iarea) = xdep(:,3)-x(:,2,iarea)        !dynamic - line 3
        2029             :
        2030   638951568 :     x_static (:,2,iarea) = xdep(:,3)                                  !static  - line 4
        2031   638951568 :     dx_static(:,2,iarea) = base_vtx(:,2,i,j,iside)-x_static(:,2,iarea)!static  - line 4
        2032             :
        2033   638951568 :     x_static (:,3,iarea) = base_vtx(:,2,i,j,iside)                         !static - line 5
        2034   638951568 :     dx_static(:,3,iarea) = base_vtx(:,1,i,j,iside)-base_vtx(:,2,i,j,iside) !static - line 5
        2035             :
        2036   212983856 :   end subroutine define_area3_center
        2037             : end module fvm_consistent_se_cslam

Generated by: LCOV version 1.14