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

Generated by: LCOV version 1.14