!    This file is part of ELPA.
!
!    The ELPA library was originally created by the ELPA consortium,
!    consisting of the following organizations:
!
!    - Max Planck Computing and Data Facility (MPCDF), formerly known as
!      Rechenzentrum Garching der Max-Planck-Gesellschaft (RZG),
!    - Bergische Universität Wuppertal, Lehrstuhl für angewandte
!      Informatik,
!    - Technische Universität München, Lehrstuhl für Informatik mit
!      Schwerpunkt Wissenschaftliches Rechnen ,
!    - Fritz-Haber-Institut, Berlin, Abt. Theorie,
!    - Max-Plack-Institut für Mathematik in den Naturwissenschaften,
!      Leipzig, Abt. Komplexe Strukutren in Biologie und Kognition,
!      and
!    - IBM Deutschland GmbH
!
!
!    More information can be found here:
!    http://elpa.mpcdf.mpg.de/
!
!    ELPA is free software: you can redistribute it and/or modify
!    it under the terms of the version 3 of the license of the
!    GNU Lesser General Public License as published by the Free
!    Software Foundation.
!
!    ELPA is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public License
!    along with ELPA.  If not, see <http://www.gnu.org/licenses/>
!
!    ELPA reflects a substantial effort on the part of the original
!    ELPA consortium, and we ask you to respect the spirit of the
!    license that we chose: i.e., please contribute any changes you
!    may have back to the original ELPA library distribution, and keep
!    any derivatives of ELPA under the same license that we chose for
!    the original distribution, the GNU Lesser General Public License.
!
!
! Copyright of the original code rests with the authors inside the ELPA
! consortium. The copyright of any additional modifications shall rest
! with their original authors, but shall adhere to the licensing terms
! distributed along with the original code in the file "COPYING".
!
! Author: Andreas Marek, MPCDF

module ELPA_utilities

   use, intrinsic :: iso_fortran_env, only : output_unit, error_unit
   use, intrinsic :: iso_c_binding
   implicit none

   private ! By default, all routines contained are private

   public :: output_unit, error_unit
   public :: check_alloc, check_alloc_CUDA_f, check_memcpy_CUDA_f, check_dealloc_CUDA_f
   public :: check_host_alloc_CUDA_f, check_host_dealloc_CUDA_f, check_host_register_CUDA_f, check_host_unregister_CUDA_f
   public :: check_memset_cuda_f
   public :: check_allocate_f, check_deallocate_f
   public :: map_global_array_index_to_local_index
   public :: pcol, prow
   public :: local_index                ! Get local index of a block cyclic distributed matrix
   public :: least_common_multiple      ! Get least common multiple

   !******
contains

   !Processor col for global col number
   pure function pcol(global_col, nblk, np_cols) result(local_col)
      use, intrinsic :: iso_c_binding, only : c_int
      implicit none
      integer(kind=c_int), intent(in) :: global_col, nblk, np_cols
      integer(kind=c_int)             :: local_col
      local_col = MOD((global_col-1)/nblk,np_cols)
   end function

   !Processor row for global row number
   pure function prow(global_row, nblk, np_rows) result(local_row)
      use, intrinsic :: iso_c_binding, only : c_int
      implicit none
      integer(kind=c_int), intent(in) :: global_row, nblk, np_rows
      integer(kind=c_int)             :: local_row
      local_row = MOD((global_row-1)/nblk,np_rows)
   end function

!-------------------------------------------------------------------------------

   function map_global_array_index_to_local_index(iGLobal, jGlobal, iLocal, jLocal , nblk, np_rows, np_cols, my_prow, my_pcol) &
      result(possible)
      use, intrinsic :: iso_c_binding, only : c_int
      implicit none

      integer(kind=c_int)              :: pi, pj, li, lj, xi, xj
      integer(kind=c_int), intent(in)  :: iGlobal, jGlobal, nblk, np_rows, np_cols, my_prow, my_pcol
      integer(kind=c_int), intent(out) :: iLocal, jLocal
      logical                       :: possible

      possible = .true.
      iLocal = 0
      jLocal = 0

      pi = prow(iGlobal, nblk, np_rows)

      if (my_prow .ne. pi) then
         possible = .false.
         return
      endif

      pj = pcol(jGlobal, nblk, np_cols)

      if (my_pcol .ne. pj) then
         possible = .false.
         return
      endif
      li = (iGlobal-1)/(np_rows*nblk) ! block number for rows
      lj = (jGlobal-1)/(np_cols*nblk) ! block number for columns

      xi = mod( (iGlobal-1),nblk)+1   ! offset in block li
      xj = mod( (jGlobal-1),nblk)+1   ! offset in block lj

      iLocal = li * nblk + xi
      jLocal = lj * nblk + xj

   end function

   integer function local_index(idx, my_proc, num_procs, nblk, iflag)

!-------------------------------------------------------------------------------
!  local_index: returns the local index for a given global index
!               If the global index has no local index on the
!               processor my_proc behaviour is defined by iflag
!
!  Parameters
!
!  idx         Global index
!
!  my_proc     Processor row/column for which to calculate the local index
!
!  num_procs   Total number of processors along row/column
!
!  nblk        Blocksize
!
!  iflag       Controls the behaviour if idx is not on local processor
!              iflag< 0 : Return last local index before that row/col
!              iflag==0 : Return 0
!              iflag> 0 : Return next local index after that row/col
!-------------------------------------------------------------------------------
      implicit none

      integer(kind=c_int) :: idx, my_proc, num_procs, nblk, iflag

      integer(kind=c_int) :: iblk

      iblk = (idx-1)/nblk  ! global block number, 0 based

      if (mod(iblk,num_procs) == my_proc) then

         ! block is local, always return local row/col number

         local_index = (iblk/num_procs)*nblk + mod(idx-1,nblk) + 1

      else

         ! non local block

         if (iflag == 0) then

            local_index = 0

         else

            local_index = (iblk/num_procs)*nblk

            if (mod(iblk,num_procs) > my_proc) local_index = local_index + nblk

            if (iflag>0) local_index = local_index + 1
         endif
      endif

   end function local_index

   integer function least_common_multiple(a, b)

      ! Returns the least common multiple of a and b
      ! There may be more efficient ways to do this, we use the most simple approach
      implicit none
      integer(kind=c_int), intent(in) :: a, b

      do least_common_multiple = a, a*(b-1), a
         if(mod(least_common_multiple,b)==0) exit
      enddo
      ! if the loop is left regularly, least_common_multiple = a*b

   end function least_common_multiple

   subroutine check_alloc(function_name, variable_name, istat, errorMessage)

      implicit none

      character(len=*), intent(in)    :: function_name
      character(len=*), intent(in)    :: variable_name
      integer(kind=c_int), intent(in)    :: istat
      character(len=*), intent(in)    :: errorMessage

      if (istat .ne. 0) then
         print *, function_name, ": error when allocating ", variable_name, " ", errorMessage
         stop 1
      endif
   end subroutine

   subroutine check_alloc_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_malloc when allocating "
         stop 1
      endif
   end subroutine

   subroutine check_dealloc_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_free when deallocating "
         stop 1
      endif
   end subroutine

   subroutine check_memcpy_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_memcpy when copying "
         stop 1
      endif
   end subroutine

   subroutine check_host_alloc_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_alloc_host when allocating "
         stop 1
      endif
   end subroutine

   subroutine check_host_dealloc_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_free_host when deallocating "
         stop 1
      endif
   end subroutine

   subroutine check_host_register_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_host_register when registering "
         stop 1
      endif
   end subroutine

   subroutine check_host_unregister_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_host_unregister when unregistering "
         stop 1
      endif
   end subroutine

   subroutine check_memset_CUDA_f(file_name, line, successCUDA)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in)    :: line
      logical                         :: successCUDA

      if (.not.(successCUDA)) then
         print *, file_name, ":", line,  " error in cuda_memset "
         stop 1
      endif
   end subroutine

   subroutine check_allocate_f(file_name, line, success, errorMessage)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in) :: line
      integer(kind=c_int)             :: success
      character(len=*)                :: errorMessage

      if ( success .ne. 0) then
         print *, file_name, ":", line,  " error in allocate: " // errorMessage
         stop 1
      endif
   end subroutine

   subroutine check_deallocate_f(file_name, line, success, errorMessage)

      implicit none

      character(len=*), intent(in)    :: file_name
      integer(kind=c_int), intent(in) :: line
      integer(kind=c_int)             :: success
      character(len=*)                :: errorMessage

      if ( success .ne. 0) then
         print *, file_name, ":", line,  " error in deallocate: " // errorMessage
         stop 1
      endif
   end subroutine
end module ELPA_utilities


text:
                  subroutine moto_z(RHS_x)

!-------------------------------------------------------------!
!             F. Califano, 2006                               !
!                                                             !
!            Eq. MOTO componente z:                           !
!                                                             !
!          Div[ BB - n(UU + de^2 jj) - P_tot * Ident ]        !
!                                                             !
! Axz = - Den * (Uix * Uiz + de2 * Uex * Uez) + Bx Bz         !
! Ayz = - Den * (Uiy * Uiz + de2 * Uey * Uez) + By Bz         !
! Azz = - Den * (Uiz * Uiz + de2 * Uez * Uez) + Bz Bz - P_tot !
!-------------------------------------------------------------!

!**************************************************
!  MPI PARALLEL VERSION: VALENTINI-FAGANELLO 2009
!   3D PARALLE VERSION: FAGANELLO 2010
!  LANDAU-FLUID VERSION  FAGANELLO 2011
!
! FLR-Landau-Fluid VERSION: Cerri 2011
!
!**************************************************

use parameter_mod
use box_mod
use deriv_mod
use fields_UJ_mod
use fields_DP_mod
use fields_EB_mod
use dom_distr_mod
use parallel_mod

IMPLICIT NONE

integer :: ix, iy, iz

REAL(dp),DIMENSION (nxl,nyl,nzl) :: RHS_x

REAL(dp), ALLOCATABLE :: zy(:),zx(:), zz(:)
REAL(dp), ALLOCATABLE :: AA(:,:,:)
REAL(dp), ALLOCATABLE :: At(:,:,:), B_inv(:,:,:)
REAL(dp), ALLOCATABLE :: Ptot(:,:,:)

ALLOCATE( At( nx, nyt, nz ) )  ! nyt ~ ny/nproc
allocate(zy(ny))
allocate(zx(nx))
allocate(zz(nz))
allocate(AA( nxl, ny, nz ))
allocate(Ptot( nxl, ny, nz ))
allocate(B_inv(nxl, ny, nz))

Ptot = Bx**2 + By**2 + Bz**2
B_inv = 1.0d0 / Ptot

 !*************************************************************!
 !  Calcolo del primo termine MOTO_z da derivare rispetto a x  !
 !*************************************************************!

AA = Bx * Bz * ( 1.0d0 - B_inv * ( pe_para + pi_para - pe_perp - pi_perp ) ) &
       - Den * (Uix * Uiz + de2 * Uex * Uez) - Gi_xz ! = Axz

 CALL traspdist( AA, At, 1 )

do iz = 1, nz
 do iy = 1, nyt
  call  derx_1( At(:,iy,iz), zx )
  At(:,iy,iz) = zx
 enddo
enddo

 CALL traspdist( RHS_x, At, -1 )

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! eventuali correz. su divB nel bulk
! CALL traspdist( Bx, At, 1 )
!do iz = 1, nz
! do iy = 1, nyt
!  call  derx_1( At(:,iy,iz), zx )
!  At(:,iy,iz) = zx
! enddo
!enddo
! CALL traspdist( AA, At, -1 )
!do iz = 1, nz
! do ix = 1, nxl
!  call  dery_1(By(ix,:,iz), zy)
!  AA(ix,:,iz) = AA(ix,:,iz) + zy
! enddo
!enddo
!RHS_x = RHS_x - Bz * AA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! N.B per ora tralascio ogni termime de2....corretto solo MHD
! N.B. in questo caso ho messo, nei punti 1,nx U e non Ui/Ue

IF (ibc .EQ. 4) then
 if ( mpime == 0 ) then

 RHS_x(1,:,:) = 0.5 * Den(1,:,:) * ( Uz(1,:,:) * (2. * L_0(1,:,:) *c_s_inv(1,:,:) + &
                alpha_2(1,:,:) * c_s_inv(1,:,:) *&
               (L_s_p(1,:,:) +  L_s_m(1,:,:)) - alpha_1(1,:,:) * f_inv(1,:,:) *&
                 (L_f_p(1,:,:) +  L_f_m(1,:,:))) &
     + a_inv(1,:,:) * ( - ay(1,:,:) * (L_a_p(1,:,:) + L_a_m(1,:,:)) + az(1,:,:) * &
   (segno(1,:,:) * alpha_1(1,:,:)* (L_s_p(1,:,:) - L_s_m(1,:,:)) &
    + ax(1,:,:) * alpha_2(1,:,:) * f_inv(1,:,:) * (L_f_p(1,:,:) - L_f_m(1,:,:)))))

!call derx_BC(Bz,grad)
!RHS_x(1,:,:) =  Bx(1,:,:) * grad(1,:,:)
do iz = 1, nz
 call dery_1(By(1,:,iz),zy)
 RHS_x(1,:,iz) = RHS_x(1,:,iz) - Bz(1,:,iz) * zy
enddo

 endif
 if ( mpime == nprow -1 ) then

 RHS_x(nxl,:,:) = 0.5 * Den(nxl,:,:) * ( Uz(nxl,:,:) * (2. * L_0(2,:,:) *c_s_inv(2,:,:) + &
                alpha_2(2,:,:) * c_s_inv(2,:,:) *&
               (L_s_p(2,:,:) +  L_s_m(2,:,:)) - alpha_1(2,:,:) * f_inv(2,:,:) *&
                 (L_f_p(2,:,:) +  L_f_m(2,:,:))) &
     + a_inv(2,:,:) * ( - ay(2,:,:) * (L_a_p(2,:,:) + L_a_m(2,:,:)) + az(2,:,:) * &
   (segno(2,:,:) * alpha_1(2,:,:)* (L_s_p(2,:,:) - L_s_m(2,:,:)) &
    + ax(2,:,:) * alpha_2(2,:,:) * f_inv(2,:,:) * (L_f_p(2,:,:) - L_f_m(2,:,:)))))

!call derx_BC(Bz,grad)
!RHS_x(nxl,:,:) = Bx(nxl,:,:) * grad(nxl,:,:)
do iz = 1, nz
 call dery_1(By(nxl,:,iz),zy)
 RHS_x(nxl,:,iz) = RHS_x(nxl,:,iz) - Bz(nxl,:,iz) * zy
enddo

 endif
ENDIF

 !*************************************************************!
 ! Calcolo del secondo termine MOTO_z da derivare rispetto a y !
 !*************************************************************!

AA = By * Bz * ( 1.0d0 - B_inv * (pe_para + pi_para - pe_perp - pi_perp ) ) &
        - Den * (Uiy * Uiz + de2 * Uey * Uez) - Gi_yz ! = Ayz

do iz = 1, nz
 do ix = 1, nxl
  call  dery_1(AA(ix,:,iz), zy)
  RHS_x(ix,:,iz) = RHS_x(ix,:,iz) + zy
 enddo
enddo

 !*************************************************************!
 ! Calcolo del terzo termine MOTO_z da derivare rispetto a z   !
 !                                                             !
 ! Ptot = pe_perp + pi_perp + B^2 / 2                          !
 !                                                             !
 ! PAI_zz = 0 per ipotesi (c.magnetico principale lungo z)     !
 ! (vedi commenti in flr_i.f90 per maggiori dettagli)          !
 !*************************************************************!

RHS_x = RHS_x * de12inv

deallocate(zx)
deallocate(zy)
deallocate(zz)
deallocate(AA)
deallocate(At)
deallocate(Ptot)
deallocate(B_inv)

end subroutine

max_stars_repo_name: sscerr/eTF

max_issues_repo_name: sscerr/eTF

max_forks_repo_name: sscerr/eTF


text:
!######################################################################
      module lsm
!######################################################################
	  SAVE
      real,allocatable,dimension(:,:):: f,g,hj_phi,hj_phi1
      real :: dt_reinit

      end module lsm
!######################################################################
      subroutine initial_lsm_3d_channel
!######################################################################
      use vars
      use lsm
      use multidata
	use mpi

      implicit none
      integer :: i,j,k,ib,tti,ttj,ttk,sn,sn1
      integer :: is,ie,js,je,ks,ke
      real :: b,dummy
      character*8 :: chb,chb1
      character*31 :: gridfile
      character *80 dummyline

	if (myrank.eq.0) then
	  write(*,'(a)') '**********************************************'
	  write(*,'(a)') '*'
	  write(*,'(a)') '*            TWO-PHASE SIMULATION'
	  if (l_lsmbase) then
	    write(*,'(a,f8.3)') '* Base case with rigid lid at z=',length
	  else
	    write(*,'(a,f8.3)') '* Initial water level: z=',length
	  endif
	  write(*,'(a)') '*'
	  write(*,'(a)') '**********************************************'
	endif

      do ib=1,nbp
        tti=dom(ib)%ttc_i;  ttj=dom(ib)%ttc_j;  ttk=dom(ib)%ttc_k
        is=dom(ib)%isp; ie=dom(ib)%iep
        js=dom(ib)%jsp; je=dom(ib)%jep
        ks=dom(ib)%ksp; ke=dom(ib)%kep

        allocate (dom(ib)%phi(tti,ttj,ttk),
     &dom(ib)%phi_init(tti,ttj,ttk),dom(ib)%phi_reinit(tti,ttj,ttk),
     &dom(ib)%phi_new(tti,ttj,ttk),dom(ib)%dphi_dx(tti,ttj,ttk),
     &dom(ib)%dphi_dy(tti,ttj,ttk),dom(ib)%dphi_dz(tti,ttj,ttk),
     &dom(ib)%s_phi0(tti,ttj,ttk),dom(ib)%h_phi(tti,ttj,ttk),
     &dom(ib)%phim(tti,ttj,ttk))
!
! Read in phi field if restarting from previous solution
!
        if (L_LSM .and. lrestart) then
          write(chb,'(i8)') dom_id(ib)
          write(chb1,'(i8)') itime_start
          sn=len(trim(adjustl(chb)))
          sn1=len(trim(adjustl(chb1)))
          chb=repeat('0',(4-sn))//trim(adjustl(chb))
          chb1=repeat('0',(6-sn1))//trim(adjustl(chb1))
          gridfile='tecout_phi_'//trim(adjustl(chb))//'_'//
     & 'initial'//'.dat'
          open (unit=703, file=gridfile)
          read (703,*) dummyline
          read (703,*) dummyline
          read (703,*) dummyline
          read (703,*) dummyline
          read (703,*) dummyline
          do k=ks-1,ke
            do j=js-1,je
              do i=is-1,ie
                read (703,73) dummy,dummy,dummy,
     & dom(ib)%phi(i,j,k),dom(ib)%phim(i,j,k),
     & dom(ib)%dens(i,j,k),dom(ib)%mu(i,j,k)
              end do
            end do
          end do
 73       format (10e25.8)
          close(703)
          dom(ib)%phi_init = 0.0
          dom(ib)%phi_new = 0.0
          dom(ib)%phi_reinit = 0.0
        else
!
! Initialise uniform phi field, if not restarting
!
          dom(ib)%phi=1.0
          dom(ib)%phi_init = 0.0
          dom(ib)%phi_new = 0.0
          dom(ib)%phi_reinit = 0.0

          if (trim(keyword).eq.'channel') then      ! Channel flow case
            do k=1,ttk
              do j=1,ttj
                do i=1,tti
!
! Set initial free surface profile, if not uniform (this is case dependent)
! Initialise length (i.e. distance of free surface from bed)
!========== Richard cube test =========================================
!          if ((dom(ib)%xc(i).ge.0).and.(dom(ib)%xc(i).le.0.035)) then
!            length=0.025
!          else if ((dom(ib)%xc(i).gt.0.035).and.
!     &             (dom(ib)%xc(i).le.0.065)) then
!            length=-0.266666667*dom(ib)%xc(i)+0.0343
!          else if ((dom(ib)%xc(i).gt.0.065).and.
!     &             (dom(ib)%xc(i).le.0.25088)) then
!            length=0.017
!          end if
!========== Sibel constriction test ===================================
!          if ((dom(ib)%xc(i).ge.0).and.(dom(ib)%xc(i).le.0.59)) then
!            length=0.076
!          else if ((dom(ib)%xc(i).gt.0.59).and.
!     &             (dom(ib)%xc(i).le.0.885)) then
!            length=-0.213559322*dom(ib)%xc(i)+0.202
!          else if ((dom(ib)%xc(i).gt.0.885).and.
!     &             (dom(ib)%xc(i).le.1.475)) then
!            length=0.013
!          end if
!======================================================================
!
! Initialise phi (free surface defined by phi=0, phi=-ve above, +ve below)
!
                  if (dom(ib)%zc(k).lt.length)   then
                    dom(ib)%phi(i,j,k) = 1.0*abs(dom(ib)%zc(k)-length)
                  else if (dom(ib)%zc(k).gt.length)   then
                    dom(ib)%phi(i,j,k) = -1.0*abs(dom(ib)%zc(k)-length)
                  else if (dom(ib)%zc(k).eq.length)  then
                    dom(ib)%phi(i,j,k) = 0.0
                  end if
                  dom(ib)%s_phi0(i,j,k)=dom(ib)%phi(i,j,k)
                  dom(ib)%phi_reinit(i,j,k)=dom(ib)%phi(i,j,k)
                end do
              end do
            end do

          else if (trim(keyword).eq.'wave') then   ! Solitary wave case

            do k=1,ttk
              do j=1,ttj
                do i=1,tti
                  b=length/(cosh(sqrt(3.*length)/2.*(dom(ib)%xc(i))))**2
                  if (dom(ib)%zc(k).lt.(b+1))   then
                    dom(ib)%phi(i,j,k) = 1.0*abs(dom(ib)%zc(k)-(b+1))
                  else if (dom(ib)%zc(k).gt.(b+1))   then
                    dom(ib)%phi(i,j,k) = -1.0*abs(dom(ib)%zc(k)-(b+1))
                  else if (dom(ib)%zc(k).eq.(b+1))  then
                    dom(ib)%phi(i,j,k) = 0.0
                  end if
                  dom(ib)%s_phi0(i,j,k)=dom(ib)%phi(i,j,k)
                end do
              end do
            end do

          end if

        end if

        dt_reinit = cfl_lsm*max(dom(ib)%dx,dom(ib)%dy,dom(ib)%dz)

      end do
!
! Set level set function phi to a signed distance function
!
      if (reinit)  call tvd_rk_reinit
!
! Define density and viscosity above, below and across the surface
!
      call heaviside

      return
      end subroutine initial_lsm_3d_channel

!######################################################################
      subroutine lsm_3d
!######################################################################
      use vars
      use lsm
      use multidata

      implicit none
      integer :: i,j,k,ib

      call tvd_rk_3step

      return
      end subroutine lsm_3d
!######################################################################
      subroutine  tvd_rk_3step
!######################################################################
!**********************************************************************
! 3 step runge kutta routine to return convected phi field
!**********************************************************************
      use vars
      use lsm
      use multidata
      use mpi

      implicit none
      integer :: i,j,k,ib
      real :: uijk,vijk,wijk,h1,h2,h3
!
! First RK step (phi --> phi_new)
!
      call dphi_a_v_3d(14)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

           if (i.eq.dom(ib)%isp) then
           uijk=1.0/16.0*(5.0*dom(ib)%u(i-1,j,k)+15.0*dom(ib)%u(i,j,k)-
     &                   5.0*dom(ib)%u(i+1,j,k)+dom(ib)%u(i+2,j,k))
           else if (i.eq.dom(ib)%iep) then
           uijk=1.0/16.0*(dom(ib)%u(i-3,j,k)-5.0*dom(ib)%u(i-2,j,k)+
     &              15.0*dom(ib)%u(i-1,j,k)+5.0*dom(ib)%u(i,j,k))
           else
           uijk=1.0/16.0*(-dom(ib)%u(i-2,j,k)+9.0*dom(ib)%u(i-1,j,k)+
     &                9.0*dom(ib)%u(i,j,k)-dom(ib)%u(i+1,j,k))
           end if

           if (j.eq.dom(ib)%jsp) then
           vijk=1.0/16.0*(5.0*dom(ib)%v(i,j-1,k)+15.0*dom(ib)%v(i,j,k)-
     &                   5.0*dom(ib)%v(i,j+1,k)+dom(ib)%v(i,j+2,k))
           else if (j.eq.dom(ib)%jep) then
           vijk=1.0/16.0*(dom(ib)%v(i,j-3,k)-5.0*dom(ib)%v(i,j-2,k)+
     &              15.0*dom(ib)%v(i,j-1,k)+5.0*dom(ib)%v(i,j,k))
           else
           vijk=1.0/16.0*(-dom(ib)%v(i,j-2,k)+9.0*dom(ib)%v(i,j-1,k)+
     &                9.0*dom(ib)%v(i,j,k)-dom(ib)%v(i,j+1,k))
           end if

           if (k.eq.dom(ib)%ksp) then
           wijk=1.0/16.0*(5.0*dom(ib)%w(i,j,k-1)+15.0*dom(ib)%w(i,j,k)-
     &                   5.0*dom(ib)%w(i,j,k+1)+dom(ib)%w(i,j,k+2))
           else if (k.eq.dom(ib)%kep) then
           wijk=1.0/16.0*(dom(ib)%w(i,j,k-3)-5.0*dom(ib)%w(i,j,k-2)+
     &              15.0*dom(ib)%w(i,j,k-1)+5.0*dom(ib)%w(i,j,k))
           else
           wijk=1.0/16.0*(-dom(ib)%w(i,j,k-2)+9.0*dom(ib)%w(i,j,k-1)+
     &                9.0*dom(ib)%w(i,j,k)-dom(ib)%w(i,j,k+1))
           end if
!
! Convection
!
           h1 = uijk*dom(ib)%dphi_dx(i,j,k)
           h2 = vijk*dom(ib)%dphi_dy(i,j,k)
           h3 = wijk*dom(ib)%dphi_dz(i,j,k)

           dom(ib)%phi_new(i,j,k) = dom(ib)%phi(i,j,k)-dt*(h1+h2+h3)

          end do
         end do
        end do

      end do

      call bound_lsm(16)    !(phi_new)
!
! 2nd RK step (phi_new --> pih_init)
!
      call dphi_a_v_3d(16)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

           if (i.eq.dom(ib)%isp) then
           uijk=1.0/16.0*(5.0*dom(ib)%u(i-1,j,k)+15.0*dom(ib)%u(i,j,k)-
     &                   5.0*dom(ib)%u(i+1,j,k)+dom(ib)%u(i+2,j,k))
           else if (i.eq.dom(ib)%iep) then
           uijk=1.0/16.0*(dom(ib)%u(i-3,j,k)-5.0*dom(ib)%u(i-2,j,k)+
     &              15.0*dom(ib)%u(i-1,j,k)+5.0*dom(ib)%u(i,j,k))
           else
           uijk=1.0/16.0*(-dom(ib)%u(i-2,j,k)+9.0*dom(ib)%u(i-1,j,k)+
     &                9.0*dom(ib)%u(i,j,k)-dom(ib)%u(i+1,j,k))
           end if

           if (j.eq.dom(ib)%jsp) then
           vijk=1.0/16.0*(5.0*dom(ib)%v(i,j-1,k)+15.0*dom(ib)%v(i,j,k)-
     &                   5.0*dom(ib)%v(i,j+1,k)+dom(ib)%v(i,j+2,k))
           else if (j.eq.dom(ib)%jep) then
           vijk=1.0/16.0*(dom(ib)%v(i,j-3,k)-5.0*dom(ib)%v(i,j-2,k)+
     &              15.0*dom(ib)%v(i,j-1,k)+5.0*dom(ib)%v(i,j,k))
           else
           vijk=1.0/16.0*(-dom(ib)%v(i,j-2,k)+9.0*dom(ib)%v(i,j-1,k)+
     &                9.0*dom(ib)%v(i,j,k)-dom(ib)%v(i,j+1,k))
           end if

           if (k.eq.dom(ib)%ksp) then
           wijk=1.0/16.0*(5.0*dom(ib)%w(i,j,k-1)+15.0*dom(ib)%w(i,j,k)-
     &                   5.0*dom(ib)%w(i,j,k+1)+dom(ib)%w(i,j,k+2))
           else if (k.eq.dom(ib)%kep) then
           wijk=1.0/16.0*(dom(ib)%w(i,j,k-3)-5.0*dom(ib)%w(i,j,k-2)+
     &              15.0*dom(ib)%w(i,j,k-1)+5.0*dom(ib)%w(i,j,k))
           else
           wijk=1.0/16.0*(-dom(ib)%w(i,j,k-2)+9.0*dom(ib)%w(i,j,k-1)+
     &                9.0*dom(ib)%w(i,j,k)-dom(ib)%w(i,j,k+1))
           end if
!
! Convection
!
           h1 = uijk*dom(ib)%dphi_dx(i,j,k)
           h2 = vijk*dom(ib)%dphi_dy(i,j,k)
           h3 = wijk*dom(ib)%dphi_dz(i,j,k)

           dom(ib)%phi_init(i,j,k) = 0.75*dom(ib)%phi(i,j,k) +
     & 0.25*dom(ib)%phi_new(i,j,k) - 0.25*dt*(h1+h2+h3)

          end do
         end do
        end do

      end do

      call bound_lsm(17)   !(phi_init)
!
! 3rd RK step (phi_init --> phi)
!
      call dphi_a_v_3d(17)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

           if (i.eq.dom(ib)%isp) then
           uijk=1.0/16.0*(5.0*dom(ib)%u(i-1,j,k)+15.0*dom(ib)%u(i,j,k)-
     &                   5.0*dom(ib)%u(i+1,j,k)+dom(ib)%u(i+2,j,k))
           else if (i.eq.dom(ib)%iep) then
           uijk=1.0/16.0*(dom(ib)%u(i-3,j,k)-5.0*dom(ib)%u(i-2,j,k)+
     &              15.0*dom(ib)%u(i-1,j,k)+5.0*dom(ib)%u(i,j,k))
           else
           uijk=1.0/16.0*(-dom(ib)%u(i-2,j,k)+9.0*dom(ib)%u(i-1,j,k)+
     &                9.0*dom(ib)%u(i,j,k)-dom(ib)%u(i+1,j,k))
           end if

           if (j.eq.dom(ib)%jsp) then
           vijk=1.0/16.0*(5.0*dom(ib)%v(i,j-1,k)+15.0*dom(ib)%v(i,j,k)-
     &                   5.0*dom(ib)%v(i,j+1,k)+dom(ib)%v(i,j+2,k))
           else if (j.eq.dom(ib)%jep) then
           vijk=1.0/16.0*(dom(ib)%v(i,j-3,k)-5.0*dom(ib)%v(i,j-2,k)+
     &              15.0*dom(ib)%v(i,j-1,k)+5.0*dom(ib)%v(i,j,k))
           else
           vijk=1.0/16.0*(-dom(ib)%v(i,j-2,k)+9.0*dom(ib)%v(i,j-1,k)+
     &                9.0*dom(ib)%v(i,j,k)-dom(ib)%v(i,j+1,k))
           end if

           if (k.eq.dom(ib)%ksp) then
           wijk=1.0/16.0*(5.0*dom(ib)%w(i,j,k-1)+15.0*dom(ib)%w(i,j,k)-
     &                   5.0*dom(ib)%w(i,j,k+1)+dom(ib)%w(i,j,k+2))
           else if (k.eq.dom(ib)%kep) then
           wijk=1.0/16.0*(dom(ib)%w(i,j,k-3)-5.0*dom(ib)%w(i,j,k-2)+
     &              15.0*dom(ib)%w(i,j,k-1)+5.0*dom(ib)%w(i,j,k))
           else
           wijk=1.0/16.0*(-dom(ib)%w(i,j,k-2)+9.0*dom(ib)%w(i,j,k-1)+
     &                9.0*dom(ib)%w(i,j,k)-dom(ib)%w(i,j,k+1))
           end if
!
! Convection
!
           h1 = uijk*dom(ib)%dphi_dx(i,j,k)
           h2 = vijk*dom(ib)%dphi_dy(i,j,k)
           h3 = wijk*dom(ib)%dphi_dz(i,j,k)

           dom(ib)%phi(i,j,k) = 1.0/3.0*dom(ib)%phi(i,j,k) +
     & 2.0/3.0*dom(ib)%phi_init(i,j,k) - 2.0/3.0*dt*(h1+h2+h3)   ! We now have updated phi field

!
! Hold level to be held constant at inflow and outflow (if required - may help with stability in inflow-outflow sims)
!
           if (trim(keyword).eq.'channel' .and. lends) then
             if (dom(ib)%iprev.lt.0) then
               if ((i.ge.dom(ib)%isu).and.(i.le.dom(ib)%isu+5)) then
                 if (dom(ib)%zc(k).lt.length)   then
                   dom(ib)%phi(i,j,k) = 1.0*abs(dom(ib)%zc(k)-length)
                 else if (dom(ib)%zc(k).gt.length)   then
                   dom(ib)%phi(i,j,k) = -1.0*abs(dom(ib)%zc(k)-length)
                 else if (dom(ib)%zc(k).eq.length)  then
                   dom(ib)%phi(i,j,k) = 0.0
                 end if
               end if
             end if
!          length=0.0354
             if  (dom(ib)%inext.lt.0) then
               if ((i.le.dom(ib)%ieu).and.(i.ge.dom(ib)%ieu-5)) then
                 if (dom(ib)%zc(k).lt.length)   then
                   dom(ib)%phi(i,j,k) = 1.0*abs(dom(ib)%zc(k)-length)
                 else if (dom(ib)%zc(k).gt.length)   then
                   dom(ib)%phi(i,j,k) = -1.0*abs(dom(ib)%zc(k)-length)
                 else if (dom(ib)%zc(k).eq.length)  then
                   dom(ib)%phi(i,j,k) = 0.0
                 end if
               end if
             end if
           end if

         end do
        end do
       end do

      end do

      call bound_lsm(14)  !(phi)
!
! Reinitialise phi field to a signed distance function
!
      if (reinit) then
        call tvd_rk_reinit     !(phi)
      end if

      return
      end
!#######################################################################
      subroutine dphi_a_v_3d(op)     !(fi)
!#######################################################################
      use vars
      use lsm
      use mpi
      use multidata
      implicit none

      integer :: i,j,k,ifi,op,ib
      real :: uijk,vijk,wijk,phi_minus_ip12j
      real :: phi_plus_ip12j,phi_minus_ijp12
      real :: phi_plus_ijp12,superbee

      call hj_weno_dxplus_3d(op)
      call hj_weno_dyplus_3d(op)
      call hj_weno_dzplus_3d(op)
      call hj_weno_dxminus_3d(op)
      call hj_weno_dyminus_3d(op)
      call hj_weno_dzminus_3d(op)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

          if (i.eq.dom(ib)%isp) then
          uijk=1.0/16.0*(5.0*dom(ib)%u(i-1,j,k)+15.0*dom(ib)%u(i,j,k)-
     &                   5.0*dom(ib)%u(i+1,j,k)+dom(ib)%u(i+2,j,k))
          else if (i.eq.dom(ib)%iep) then
          uijk=1.0/16.0*(dom(ib)%u(i-3,j,k)-5.0*dom(ib)%u(i-2,j,k)+
     &              15.0*dom(ib)%u(i-1,j,k)+5.0*dom(ib)%u(i,j,k))
          else
          uijk=1.0/16.0*(-dom(ib)%u(i-2,j,k)+9.0*dom(ib)%u(i-1,j,k)+
     &                9.0*dom(ib)%u(i,j,k)-dom(ib)%u(i+1,j,k))
          end if

          if (j.eq.dom(ib)%jsp) then
          vijk=1.0/16.0*(5.0*dom(ib)%v(i,j-1,k)+15.0*dom(ib)%v(i,j,k)-
     &                   5.0*dom(ib)%v(i,j+1,k)+dom(ib)%v(i,j+2,k))
          else if (j.eq.dom(ib)%jep) then
          vijk=1.0/16.0*(dom(ib)%v(i,j-3,k)-5.0*dom(ib)%v(i,j-2,k)+
     &              15.0*dom(ib)%v(i,j-1,k)+5.0*dom(ib)%v(i,j,k))
          else
          vijk=1.0/16.0*(-dom(ib)%v(i,j-2,k)+9.0*dom(ib)%v(i,j-1,k)+
     &                9.0*dom(ib)%v(i,j,k)-dom(ib)%v(i,j+1,k))
          end if

          if (k.eq.dom(ib)%ksp) then
          wijk=1.0/16.0*(5.0*dom(ib)%w(i,j,k-1)+15.0*dom(ib)%w(i,j,k)-
     &                   5.0*dom(ib)%w(i,j,k+1)+dom(ib)%w(i,j,k+2))
          else if (k.eq.dom(ib)%kep) then
          wijk=1.0/16.0*(dom(ib)%w(i,j,k-3)-5.0*dom(ib)%w(i,j,k-2)+
     &              15.0*dom(ib)%w(i,j,k-1)+5.0*dom(ib)%w(i,j,k))
          else
          wijk=1.0/16.0*(-dom(ib)%w(i,j,k-2)+9.0*dom(ib)%w(i,j,k-1)+
     &                9.0*dom(ib)%w(i,j,k)-dom(ib)%w(i,j,k+1))
          end if

           if (uijk.gt.0.0)
     &        dom(ib)%dphi_dx(i,j,k) = dom(ib)%dphi_dxminus(i,j,k)
           if (uijk.lt.0.0)
     &        dom(ib)%dphi_dx(i,j,k) = dom(ib)%dphi_dxplus(i,j,k)
           if (uijk.eq.0.0)
     &        dom(ib)%dphi_dx(i,j,k) = 0.0

           if (vijk.gt.0.0)
     &        dom(ib)%dphi_dy(i,j,k) = dom(ib)%dphi_dyminus(i,j,k)
           if (vijk.lt.0.0)
     &        dom(ib)%dphi_dy(i,j,k) = dom(ib)%dphi_dyplus(i,j,k)
           if (vijk.eq.0.0)
     &        dom(ib)%dphi_dy(i,j,k) = 0.0

           if (wijk.gt.0.0)
     &        dom(ib)%dphi_dz(i,j,k) = dom(ib)%dphi_dzminus(i,j,k)
           if (wijk.lt.0.0)
     &        dom(ib)%dphi_dz(i,j,k) = dom(ib)%dphi_dzplus(i,j,k)
           if (wijk.eq.0.0)
     &        dom(ib)%dphi_dz(i,j,k) = 0.0

         end do
       end do
      end do

      end do

      return
      end subroutine dphi_a_v_3d
!######################################################################
      subroutine tvd_rk_reinit
!######################################################################
!**********************************************************************
! 3 step runge kutta routine which returns re-initialised phi field
!**********************************************************************
      use vars
      use lsm
      use multidata
      use mpi

      implicit none
      integer :: i,j,k,it,ib
      logical :: bool
      real :: max_abs,abs_dphi,abs_phidiff,max_phidiff
      real :: local_max_abs,local_max_phidiff

      bool=.false.
      it=0

      do while (bool.eq..false.)
!
! 1st step
!
      call dphi_for_reinit(14)      ! (phi)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

        abs_dphi = sqrt(dom(ib)%dphi_dx(i,j,k)**2+
     & dom(ib)%dphi_dy(i,j,k)**2+dom(ib)%dphi_dz(i,j,k)**2)

!	  dom(ib)%abs_dphi_check(i,j,k)=abs_dphi

        dom(ib)%s_phi0(i,j,k) = dom(ib)%phi(i,j,k)/
     &      sqrt(dom(ib)%phi(i,j,k)**2+(abs_dphi**2)*dom(ib)%dx**2)

        dom(ib)%phi_reinit(i,j,k) =dom(ib)%phi(i,j,k) +
     &  dt_reinit*(dom(ib)%s_phi0(i,j,k)-
     &   dom(ib)%s_phi0(i,j,k)*abs_dphi)

        end do
       end do
      end do

      end do

      call bound_lsm(15) ! (phi_reinit)
!
! 2nd step
!
      call dphi_for_reinit(15) ! (phi_reinit)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

        abs_dphi = sqrt(dom(ib)%dphi_dx(i,j,k)**2+
     & dom(ib)%dphi_dy(i,j,k)**2+dom(ib)%dphi_dz(i,j,k)**2)

        dom(ib)%s_phi0(i,j,k) = dom(ib)%phi_reinit(i,j,k)/
     &sqrt(dom(ib)%phi_reinit(i,j,k)**2+(abs_dphi**2)*dom(ib)%dx**2)

        dom(ib)%phi_reinit(i,j,k) = 0.75*dom(ib)%phi(i,j,k)+
     &0.25*dom(ib)%phi_reinit(i,j,k)+
     &0.25*dt_reinit*(dom(ib)%s_phi0(i,j,k)-
     &dom(ib)%s_phi0(i,j,k)*abs_dphi)

          end do
         end do
        end do

      end do

      call bound_lsm(15) ! (phi_reinit)
!
! 3rd step
!
      call dphi_for_reinit(15) ! (phi_reinit)

      do ib=1,nbp

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

        abs_dphi = sqrt(dom(ib)%dphi_dx(i,j,k)**2+
     & dom(ib)%dphi_dy(i,j,k)**2+dom(ib)%dphi_dz(i,j,k)**2)

        dom(ib)%s_phi0(i,j,k) = dom(ib)%phi_reinit(i,j,k)/
     &sqrt(dom(ib)%phi_reinit(i,j,k)**2+(abs_dphi**2)*dom(ib)%dx**2)

        dom(ib)%phi_reinit(i,j,k) =1.0/3.0*dom(ib)%phi(i,j,k)+
     & 2.0/3.0*dom(ib)%phi_reinit(i,j,k)+
     & 2.0/3.0*dt_reinit*(dom(ib)%s_phi0(i,j,k)-
     & dom(ib)%s_phi0(i,j,k)*abs_dphi)

          end do
         end do
        end do

      end do

      call bound_lsm(15) ! (phi_reinit)

!
! Check convergence
!
      call dphi_for_reinit(15) ! (phi_reinit)

      max_abs = 0.0
      max_phidiff = 0.0

      do ib=1,nbp

       do k=dom(ib)%ksp,dom(ib)%kep
        do i=dom(ib)%isp,dom(ib)%iep
         do j=dom(ib)%jsp,dom(ib)%jep

          abs_dphi = sqrt(dom(ib)%dphi_dx(i,j,k)**2+
     & dom(ib)%dphi_dy(i,j,k)**2+dom(ib)%dphi_dz(i,j,k)**2)

          max_abs  = max(max_abs,abs_dphi)
          abs_phidiff=abs(dom(ib)%phi(i,j,k)-dom(ib)%phi_reinit(i,j,k))
          max_phidiff = max(max_phidiff,abs_phidiff)

         end do
        end do
       end do

      end do

      local_max_abs = max_abs

      call mpi_allreduce(local_max_abs,max_abs,1,mpi_flt,mpi_max,
     &mpi_comm_world,ierr)

      local_max_phidiff = max_phidiff

      call mpi_allreduce(local_max_phidiff,max_phidiff,1,mpi_flt,
     &mpi_max,mpi_comm_world,ierr)

      if ((max_phidiff.lt.reldif_lsm).and.(it.ge.1)) then
        bool=.true.
      else
        if (it.ge.ntime_reinit) then
          bool=.true.
        else
          it=it+1
        end if
      end if

      do ib=1,nbp
        do k=dom(ib)%ksp,dom(ib)%kep
          do i=dom(ib)%isp,dom(ib)%iep
            do j=dom(ib)%jsp,dom(ib)%jep
              dom(ib)%phi(i,j,k) = dom(ib)%phi_reinit(i,j,k)
            end do
          end do
        end do
      end do

      call bound_lsm(14) ! (phi)

      end do

      if (myrank.eq.0) then
        write(*,*) 'norm v (reinit)', max_abs, 'needed steps', it
        write(numfile3,'(i8,f18.8,i8,2f18.8)') ntime,max_abs,it,ctime,
     & dt
      endif

      return
      end subroutine tvd_rk_reinit
!######################################################################
      subroutine dphi_for_reinit(op)
!######################################################################
      use vars
      use lsm
      use multidata

      implicit none

      integer :: i,j,k,ifi,op,ib
      real :: s_phi012
      double precision, pointer, dimension(:,:,:)::fi
      real :: lsv,lssig,xm,xp,ym,yp,zm,zp

          call hj_weno_dxplus_3d(op)
          call hj_weno_dyplus_3d(op)
          call hj_weno_dzplus_3d(op)
          call hj_weno_dxminus_3d(op)
          call hj_weno_dyminus_3d(op)
          call hj_weno_dzminus_3d(op)

      do ib=1,nbp

       select case (op)
        case (14) ! (phi)
          fi => dom(ib)%phi
        case (15) ! (phi_reinit)
          fi => dom(ib)%phi_reinit
       end select

        do k=dom(ib)%ksp,dom(ib)%kep
         do i=dom(ib)%isp,dom(ib)%iep
          do j=dom(ib)%jsp,dom(ib)%jep

              lsv = fi(i,j,k)
              lssig = dom(ib)%s_phi0(i,j,k) ! s_phi0=phi initially

          xm=dom(ib)%dphi_dxminus(i,j,k)
          xp=dom(ib)%dphi_dxplus(i,j,k)
          ym=dom(ib)%dphi_dyminus(i,j,k)
          yp=dom(ib)%dphi_dyplus(i,j,k)
          zm=dom(ib)%dphi_dzminus(i,j,k)
          zp=dom(ib)%dphi_dzplus(i,j,k)

      if ((xm*lssig.gt.0.0).and.(xp*lssig.gt.-xm*lssig)) then
            dom(ib)%dphi_dx(i,j,k)=dom(ib)%dphi_dxminus(i,j,k)
      end if

      if((xp*lssig.lt.0.0).and.(xm*lssig.lt.-xp*lssig)) then
            dom(ib)%dphi_dx(i,j,k)=dom(ib)%dphi_dxplus(i,j,k)
      end if

      if ((xp*lssig.gt.0.0).and.(xm*lssig.lt.0.0)) then
            dom(ib)%dphi_dx(i,j,k)=0.0
      end if

      if ((ym*lssig.gt.0.0).and.(yp*lssig.gt.-ym*lssig)) then
            dom(ib)%dphi_dy(i,j,k)=dom(ib)%dphi_dyminus(i,j,k)
      end if

      if ((yp*lssig.lt.0.0).and.(ym*lssig.lt.-yp*lssig)) then
           dom(ib)%dphi_dy(i,j,k)=dom(ib)%dphi_dyplus(i,j,k)
      end if

      if ((yp*lssig.gt.0.0).and.(ym*lssig.lt.0.0)) then
          dom(ib)%dphi_dy(i,j,k) = 0.0
      end if

      if ((zm*lssig.gt.0.0).and.(zp*lssig.gt.-zm*lssig)) then
            dom(ib)%dphi_dz(i,j,k)=dom(ib)%dphi_dzminus(i,j,k)
      end if

      if((zp*lssig.lt.0.0).and.(zm*lssig.lt.-zp*lssig)) then
            dom(ib)%dphi_dz(i,j,k)=dom(ib)%dphi_dzplus(i,j,k)
      end if

      if ((zp*lssig.gt.0.0).and.(zm*lssig.lt.0.0)) then
            dom(ib)%dphi_dz(i,j,k)=0.0
      end if

          end do
         end do
        end do

       end do

      return
      end
!######################################################################
      subroutine hj_weno_dxplus_3d(op)
!######################################################################
      use vars
      use multidata

      implicit none

      integer :: i,j,k,l,k_star,np,nq,nr,op,ib
      double precision, pointer, dimension(:,:,:)::fi
      real :: q1,q2,q3,c_star,v1,v2,v3,v4,v5,v11,v22
	real :: v33,v44,v55,s1,s2,s3,a1,a2,a3,w1,w2,w3
	real :: e

      do ib=1,nbp

       select case (op)
        case (1)
          fi => dom(ib)%u
          np = dom(ib)%niul;  nq = dom(ib)%njul;  nr = dom(ib)%nkul
        case (2)
          fi => dom(ib)%v
          np = dom(ib)%nivl;  nq = dom(ib)%njvl;  nr = dom(ib)%nkvl
        case (3)
          fi => dom(ib)%w
          np = dom(ib)%niwl;  nq = dom(ib)%njwl;  nr = dom(ib)%nkwl
        case (14) ! (phi)
          fi => dom(ib)%phi
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (15) ! (phi_reinit)
          fi => dom(ib)%phi_reinit
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (16) ! (phi_new)
          fi => dom(ib)%phi_new
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (17) ! (phi_init)
          fi => dom(ib)%phi_init
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
       end select

      if (dom(ib)%iprev.lt.0 .and. dom(ib)%bc_west.ne.5) then
       do k=1,nr
        do j=1,nq
            fi(2,j,k)    = fi(3,j,k)
            fi(1,j,k)    = fi(2,j,k)
        end do
       end do
      else if (dom(ib)%inext.lt.0 .and. dom(ib)%bc_east.ne.5) then
       do k=1,nr
        do j=1,nq
            fi(np-1,j,k) = fi(np-2,j,k)
            fi(np,j,k)   = fi(np-1,j,k)
        end do
       end do
      end if
!
! Compute 1st order spatial derivative
!
      do k=1,nr
       do i=1,np-1
        do j=1,nq
           dom(ib)%d1(i,j,k) = (fi(i+1,j,k) - fi(i,j,k))/dom(ib)%dx
        end do
       end do
      end do
!
! Compute 5th order derivative (+ve direction)
!
      dom(ib)%dphi_dxplus = 0.0

      do k=1,nr
       do i=1,np-6
        do j=1,nq

            l=i

           v1 = dom(ib)%d1(l+5,j,k)
           v2 = dom(ib)%d1(l+4,j,k)
           v3 = dom(ib)%d1(l+3,j,k)
           v4 = dom(ib)%d1(l+2,j,k)
           v5 = dom(ib)%d1(l+1,j,k)

           v11 = v1**2
           v22 = v2**2
           v33 = v3**2
           v44 = v4**2
           v55 = v5**2

           e = (1.0e-6) * max(v11,v22,v33,v44,v55) + 1.0d-99

      s1 = (13.0/12.0)*(v1-2.0*v2+v3)**2 + 0.25*(v1-4.0*v2+3.0*v3)**2
      s2 = (13.0/12.0)*(v2-2.0*v3+v4)**2 + 0.25*(v2-v4)**2
      s3 = (13.0/12.0)*(v3-2.0*v4+v5)**2 + 0.25*(3.0*v3-4.0*v4+v5)**2

           a1 = 0.1/(e+s1)**2
           a2 = 0.6/(e+s2)**2
           a3 = 0.3/(e+s3)**2

           w1 = a1/(a1+a2+a3)
           w2 = a2/(a1+a2+a3)
           w3 = a3/(a1+a2+a3)

      dom(ib)%dphi_dxplus(i+3,j,k) = w1*(v1*(1.0/3.0)-v2*(7.0/6.0)+
     &v3*(11.0/6.0)) + w2*(-v2*(1.0/6.0)+v3*(5.0/6.0)+v4*(1.0/3.0))+
     &w3*(v3*(1.0/3.0)+v4*(5.0/6.0)-v5*(1.0/6.0))

        end do
       end do
      end do

      end do

      return
      end subroutine hj_weno_dxplus_3d
!######################################################################
      subroutine hj_weno_dxminus_3d(op)
!######################################################################
      use vars
      use multidata

      implicit none

      integer :: i,j,k,l,k_star,np,nq,nr,op,ib
      double precision, pointer, dimension(:,:,:)::fi
      real :: q1,q2,q3,c_star,v1,v2,v3,v4,v5,v11,v22
	real :: v33,v44,v55,s1,s2,s3,a1,a2,a3,w1,w2,w3
	real :: e

      do ib=1,nbp

       select case (op)
        case (1)
          fi => dom(ib)%u
          np = dom(ib)%niul;  nq = dom(ib)%njul;  nr = dom(ib)%nkul
        case (2)
          fi => dom(ib)%v
          np = dom(ib)%nivl;  nq = dom(ib)%njvl;  nr = dom(ib)%nkvl
        case (3)
          fi => dom(ib)%w
          np = dom(ib)%niwl;  nq = dom(ib)%njwl;  nr = dom(ib)%nkwl
        case (14)  ! (phi)
          fi => dom(ib)%phi
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (15)  ! (phi_reinit)
          fi => dom(ib)%phi_reinit
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (16)  ! (phi_new)
          fi => dom(ib)%phi_new
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (17)  ! (phi_init)
          fi => dom(ib)%phi_init
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
       end select

      if (dom(ib)%iprev.lt.0 .and. dom(ib)%bc_west.ne.5) then
       do k=1,nr
        do j=1,nq
            fi(2,j,k)    = fi(3,j,k)
            fi(1,j,k)    = fi(2,j,k)
        end do
       end do
      else if (dom(ib)%inext.lt.0 .and. dom(ib)%bc_east.ne.5) then
       do k=1,nr
        do j=1,nq
            fi(np-1,j,k) = fi(np-2,j,k)
            fi(np,j,k)   = fi(np-1,j,k)
        end do
       end do
      end if
!
! Compute 1st order spatial derivative
!
      do k=1,nr
       do i=1,np-1
        do j=1,nq
           dom(ib)%d1(i,j,k) = (fi(i+1,j,k) - fi(i,j,k))/dom(ib)%dx
        end do
       end do
      end do
!
! Compute 5th order derivative (-ve direction)
!
      dom(ib)%dphi_dxminus = 0.0

      do k=1,nr
       do i=1,np-6
        do j=1,nq

           l=i-1

           v1 = dom(ib)%d1(l+1,j,k)
           v2 = dom(ib)%d1(l+2,j,k)
           v3 = dom(ib)%d1(l+3,j,k)
           v4 = dom(ib)%d1(l+4,j,k)
           v5 = dom(ib)%d1(l+5,j,k)

           v11 = v1**2
           v22 = v2**2
           v33 = v3**2
           v44 = v4**2
           v55 = v5**2

           e = (1.0e-6) * max(v11,v22,v33,v44,v55) + 1.0d-99

      s1 = (13.0/12.0)*(v1-2.0*v2+v3)**2 + 0.25*(v1-4.0*v2+3.0*v3)**2
      s2 = (13.0/12.0)*(v2-2.0*v3+v4)**2 + 0.25*(v2-v4)**2
      s3 = (13.0/12.0)*(v3-2.0*v4+v5)**2 + 0.25*(3.0*v3-4.0*v4+v5)**2

           a1 = 0.1/(e+s1)**2
           a2 = 0.6/(e+s2)**2
           a3 = 0.3/(e+s3)**2

           w1 = a1/(a1+a2+a3)
           w2 = a2/(a1+a2+a3)
           w3 = a3/(a1+a2+a3)

      dom(ib)%dphi_dxminus(i+3,j,k) = w1*(v1*(1.0/3.0)-v2*(7.0/6.0)+
     &v3*(11.0/6.0)) + w2*(-v2*(1.0/6.0)+v3*(5.0/6.0)+v4*(1.0/3.0))+
     &w3*(v3*(1.0/3.0)+v4*(5.0/6.0)-v5*(1.0/6.0))

        end do
       end do
      end do

      end do

      return
      end subroutine hj_weno_dxminus_3d
!######################################################################
      subroutine hj_weno_dyplus_3d(op)
!######################################################################
      use vars
      use multidata

      implicit none

      integer :: i,j,k,l,k_star,np,nq,nr,op,ib
      double precision, pointer, dimension(:,:,:)::fi
      real :: q1,q2,q3,c_star,v1,v2,v3,v4,v5,v11,v22
	real :: v33,v44,v55,s1,s2,s3,a1,a2,a3,w1,w2,w3
	real :: e

      do ib=1,nbp

       select case (op)
        case (1)
          fi => dom(ib)%u
          np = dom(ib)%niul;  nq = dom(ib)%njul;  nr = dom(ib)%nkul
        case (2)
          fi => dom(ib)%v
          np = dom(ib)%nivl;  nq = dom(ib)%njvl;  nr = dom(ib)%nkvl
        case (3)
          fi => dom(ib)%w
          np = dom(ib)%niwl;  nq = dom(ib)%njwl;  nr = dom(ib)%nkwl
        case (14)  ! (phi)
          fi => dom(ib)%phi
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (15)  ! (phi_reinit)
          fi => dom(ib)%phi_reinit
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (16)  ! (phi_new)
          fi => dom(ib)%phi_new
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (17)  ! (phi_init)
          fi => dom(ib)%phi_init
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
       end select

      if (dom(ib)%jprev.lt.0 .and. dom(ib)%bc_south.ne.5) then
       do k=1,nr
        do i=1,np
            fi(i,2,k)    = fi(i,3,k)
            fi(i,1,k)    = fi(i,2,k)
        end do
       end do
      else if (dom(ib)%jnext.lt.0 .and. dom(ib)%bc_north.ne.5) then
       do k=1,nr
        do i=1,np
            fi(i,nq-1,k) = fi(i,nq-2,k)
            fi(i,nq,k)   = fi(i,nq-1,k)
        end do
       end do
      end if
!
! Compute 1st order spatial derivative
!
      do k=1,nr
       do i=1,np
        do j=1,nq-1
           dom(ib)%d1(i,j,k) = (fi(i,j+1,k) - fi(i,j,k))/dom(ib)%dy
        end do
       end do
      end do
!
! Compute 5th order derivative (+ve direction)
!
      dom(ib)%dphi_dyplus = 0.0

      do k=1,nr
       do i=1,np
        do j=1,nq-6
           l=j

           v1 = dom(ib)%d1(i,l+5,k)
           v2 = dom(ib)%d1(i,l+4,k)
           v3 = dom(ib)%d1(i,l+3,k)
           v4 = dom(ib)%d1(i,l+2,k)
           v5 = dom(ib)%d1(i,l+1,k)

           v11 = v1**2
           v22 = v2**2
           v33 = v3**2
           v44 = v4**2
           v55 = v5**2

           e = (1.0e-6) * max(v11,v22,v33,v44,v55) + 1.0d-99

      s1 = (13.0/12.0)*(v1-2.0*v2+v3)**2 + 0.25*(v1-4.0*v2+3.0*v3)**2
      s2 = (13.0/12.0)*(v2-2.0*v3+v4)**2 + 0.25*(v2-v4)**2
      s3 = (13.0/12.0)*(v3-2.0*v4+v5)**2 + 0.25*(3.0*v3-4.0*v4+v5)**2

           a1 = 0.1/(e+s1)**2
           a2 = 0.6/(e+s2)**2
           a3 = 0.3/(e+s3)**2

           w1 = a1/(a1+a2+a3)
           w2 = a2/(a1+a2+a3)
           w3 = a3/(a1+a2+a3)

      dom(ib)%dphi_dyplus(i,j+3,k) = w1*(v1*(1.0/3.0)-v2*(7.0/6.0)+
     &v3*(11.0/6.0)) + w2*(-v2*(1.0/6.0)+v3*(5.0/6.0)+v4*(1.0/3.0))+
     &w3*(v3*(1.0/3.0)+v4*(5.0/6.0)-v5*(1.0/6.0))

        end do
       end do
      end do

      end do

      return
      end subroutine hj_weno_dyplus_3d
!######################################################################
      subroutine hj_weno_dyminus_3d(op)
!######################################################################
      use vars
      use multidata

      implicit none

      integer :: i,j,k,l,k_star,np,nq,nr,op,ib
      double precision, pointer, dimension(:,:,:)::fi
      real :: q1,q2,q3,c_star,v1,v2,v3,v4,v5,v11,v22
	real :: v33,v44,v55,s1,s2,s3,a1,a2,a3,w1,w2,w3
	real :: e

      do ib=1,nbp

       select case (op)
        case (1)
          fi => dom(ib)%u
          np = dom(ib)%niul;  nq = dom(ib)%njul;  nr = dom(ib)%nkul
        case (2)
          fi => dom(ib)%v
          np = dom(ib)%nivl;  nq = dom(ib)%njvl;  nr = dom(ib)%nkvl
        case (3)
          fi => dom(ib)%w
          np = dom(ib)%niwl;  nq = dom(ib)%njwl;  nr = dom(ib)%nkwl
        case (14) ! (phi)
          fi => dom(ib)%phi
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (15) ! (phi_reinit)
          fi => dom(ib)%phi_reinit
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (16) ! (phi_new)
          fi => dom(ib)%phi_new
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (17) ! (phi_init)
          fi => dom(ib)%phi_init
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
       end select

      if (dom(ib)%jprev.lt.0 .and. dom(ib)%bc_south.ne.5) then
       do k=1,nr
        do i=1,np
            fi(i,2,k)    = fi(i,3,k)
            fi(i,1,k)    = fi(i,2,k)
        end do
       end do
      else if (dom(ib)%jnext.lt.0 .and. dom(ib)%bc_north.ne.5) then
       do k=1,nr
        do i=1,np
            fi(i,nq-1,k) = fi(i,nq-2,k)
            fi(i,nq,k)   = fi(i,nq-1,k)
        end do
       end do
      end if
!
! Compute 1st order spatial derivative
!
      do k=1,nr
       do i=1,np
        do j=1,nq-1
           dom(ib)%d1(i,j,k) = (fi(i,j+1,k) - fi(i,j,k))/dom(ib)%dy
        end do
       end do
      end do
!
! Compute 5th order derivative (-ve direction)
!
      dom(ib)%dphi_dyminus = 0.0

      do k=1,nr
       do i=1,np
        do j=1,nq-6

           l=j-1

           v1 = dom(ib)%d1(i,l+1,k)
           v2 = dom(ib)%d1(i,l+2,k)
           v3 = dom(ib)%d1(i,l+3,k)
           v4 = dom(ib)%d1(i,l+4,k)
           v5 = dom(ib)%d1(i,l+5,k)

           v11 = v1**2
           v22 = v2**2
           v33 = v3**2
           v44 = v4**2
           v55 = v5**2

           e = (1.0e-6) * max(v11,v22,v33,v44,v55) + 1.0d-99

      s1 = (13.0/12.0)*(v1-2.0*v2+v3)**2 + 0.25*(v1-4.0*v2+3.0*v3)**2
      s2 = (13.0/12.0)*(v2-2.0*v3+v4)**2 + 0.25*(v2-v4)**2
      s3 = (13.0/12.0)*(v3-2.0*v4+v5)**2 + 0.25*(3.0*v3-4.0*v4+v5)**2

           a1 = 0.1/(e+s1)**2
           a2 = 0.6/(e+s2)**2
           a3 = 0.3/(e+s3)**2

           w1 = a1/(a1+a2+a3)
           w2 = a2/(a1+a2+a3)
           w3 = a3/(a1+a2+a3)

      dom(ib)%dphi_dyminus(i,j+3,k) = w1*(v1*(1.0/3.0)-v2*(7.0/6.0)+
     &v3*(11.0/6.0)) + w2*(-v2*(1.0/6.0)+v3*(5.0/6.0)+v4*(1.0/3.0))+
     &w3*(v3*(1.0/3.0)+v4*(5.0/6.0)-v5*(1.0/6.0))

        end do
       end do
      end do

      end do

      return
      end subroutine hj_weno_dyminus_3d
!######################################################################
      subroutine hj_weno_dzplus_3d(op)
!######################################################################
      use vars
      use multidata

      implicit none

      integer :: i,j,k,l,k_star,np,nq,nr,op,ib
      double precision, pointer, dimension(:,:,:)::fi
      real :: e,q1,q2,q3,c_star,v1,v2,v3,v4,v5,v11,v22
	real :: v33,v44,v55,s1,s2,s3,a1,a2,a3,w1,w2,w3

      do ib=1,nbp

       select case (op)
        case (1)
          fi => dom(ib)%u
          np = dom(ib)%niul;  nq = dom(ib)%njul;  nr = dom(ib)%nkul
        case (2)
          fi => dom(ib)%v
          np = dom(ib)%nivl;  nq = dom(ib)%njvl;  nr = dom(ib)%nkvl
        case (3)
          fi => dom(ib)%w
          np = dom(ib)%niwl;  nq = dom(ib)%njwl;  nr = dom(ib)%nkwl
        case (14) ! (phi)
          fi => dom(ib)%phi
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (15) ! (phi_reinit)
          fi => dom(ib)%phi_reinit
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (16) ! (phi_new)
          fi => dom(ib)%phi_new
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (17) ! (phi_init)
          fi => dom(ib)%phi_init
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
       end select

      if (dom(ib)%kprev.lt.0 .and. dom(ib)%bc_bottom.ne.5) then
       do i=1,np
        do j=1,nq
            fi(i,j,2)    = fi(i,j,3)
            fi(i,j,1)    = fi(i,j,2)
        end do
       end do
      else if (dom(ib)%knext.lt.0 .and. dom(ib)%bc_top.ne.5) then
       do i=1,np
        do j=1,nq
            fi(i,j,nr-1) = fi(i,j,nr-2)
            fi(i,j,nr)   = fi(i,j,nr-1)
        end do
       end do
      end if
!
! Compute 1st order spatial derivative
!
      do k=1,nr-1
       do i=1,np
        do j=1,nq
           dom(ib)%d1(i,j,k) = (fi(i,j,k+1) - fi(i,j,k))/dom(ib)%dz
        end do
       end do
      end do
!
! Compute 5th order derivative (-ve direction)
!
      dom(ib)%dphi_dzplus = 0.0

      do k=1,nr-6
       do i=1,np
        do j=1,nq

            l=k

           v1 = dom(ib)%d1(i,j,l+5)
           v2 = dom(ib)%d1(i,j,l+4)
           v3 = dom(ib)%d1(i,j,l+3)
           v4 = dom(ib)%d1(i,j,l+2)
           v5 = dom(ib)%d1(i,j,l+1)

           v11 = v1**2
           v22 = v2**2
           v33 = v3**2
           v44 = v4**2
           v55 = v5**2

           e = (1.0e-6) * max(v11,v22,v33,v44,v55) + 1.0d-99

      s1 = (13.0/12.0)*(v1-2.0*v2+v3)**2 + 0.25*(v1-4.0*v2+3.0*v3)**2
      s2 = (13.0/12.0)*(v2-2.0*v3+v4)**2 + 0.25*(v2-v4)**2
      s3 = (13.0/12.0)*(v3-2.0*v4+v5)**2 + 0.25*(3.0*v3-4.0*v4+v5)**2

           a1 = 0.1/(e+s1)**2
           a2 = 0.6/(e+s2)**2
           a3 = 0.3/(e+s3)**2

           w1 = a1/(a1+a2+a3)
           w2 = a2/(a1+a2+a3)
           w3 = a3/(a1+a2+a3)

      dom(ib)%dphi_dzplus(i,j,k+3) = w1*(v1*(1.0/3.0)-v2*(7.0/6.0)+
     &v3*(11.0/6.0)) + w2*(-v2*(1.0/6.0)+v3*(5.0/6.0)+v4*(1.0/3.0))+
     &w3*(v3*(1.0/3.0)+v4*(5.0/6.0)-v5*(1.0/6.0))

        end do
       end do
      end do

      end do

      return
      end subroutine hj_weno_dzplus_3d
!######################################################################
      subroutine hj_weno_dzminus_3d(op)
!######################################################################
      use vars
      use multidata

      implicit none

      integer :: i,j,k,l,k_star,np,nq,nr,op,ib
      double precision, pointer, dimension(:,:,:)::fi
      real :: q1,q2,q3,c_star,v1,v2,v3,v4,v5,v11,v22
	real :: v33,v44,v55,s1,s2,s3,a1,a2,a3,w1,w2,w3
	real :: e

      do ib=1,nbp

       select case (op)
        case (1)
          fi => dom(ib)%u
          np = dom(ib)%niul;  nq = dom(ib)%njul;  nr = dom(ib)%nkul
        case (2)
          fi => dom(ib)%v
          np = dom(ib)%nivl;  nq = dom(ib)%njvl;  nr = dom(ib)%nkvl
        case (3)
          fi => dom(ib)%w
          np = dom(ib)%niwl;  nq = dom(ib)%njwl;  nr = dom(ib)%nkwl
        case (14)  ! (phi)
          fi => dom(ib)%phi
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (15)  ! (phi_reinit)
          fi => dom(ib)%phi_reinit
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (16)  ! (phi_new)
          fi => dom(ib)%phi_new
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
        case (17)  ! (phi_init)
          fi => dom(ib)%phi_init
          np = dom(ib)%ttc_i; nq = dom(ib)%ttc_j; nr = dom(ib)%ttc_k
       end select

      if (dom(ib)%kprev.lt.0 .and. dom(ib)%bc_bottom.ne.5) then
       do i=1,np
        do j=1,nq
            fi(i,j,2)    = fi(i,j,3)
            fi(i,j,1)    = fi(i,j,2)
        end do
       end do
      else if (dom(ib)%knext.lt.0 .and. dom(ib)%bc_top.ne.5) then
       do i=1,np
        do j=1,nq
            fi(i,j,nr-1) = fi(i,j,nr-2)
            fi(i,j,nr)   = fi(i,j,nr-1)
        end do
       end do
      end if
!
! Compute 1st order spatial derivative
!
      do k=1,nr-1
       do i=1,np
        do j=1,nq
           dom(ib)%d1(i,j,k) = (fi(i,j,k+1) - fi(i,j,k))/dom(ib)%dz
        end do
       end do
      end do
!
! Compute 5th order derivative (-ve direction)
!
      dom(ib)%dphi_dzminus = 0.0

      do k=1,nr-6
       do i=1,np
        do j=1,nq

           l=k-1

           v1 = dom(ib)%d1(i,j,l+1)
           v2 = dom(ib)%d1(i,j,l+2)
           v3 = dom(ib)%d1(i,j,l+3)
           v4 = dom(ib)%d1(i,j,l+4)
           v5 = dom(ib)%d1(i,j,l+5)

           v11 = v1**2
           v22 = v2**2
           v33 = v3**2
           v44 = v4**2
           v55 = v5**2

           e = (1.0e-6) * max(v11,v22,v33,v44,v55) + 1.0d-99

      s1 = (13.0/12.0)*(v1-2.0*v2+v3)**2 + 0.25*(v1-4.0*v2+3.0*v3)**2
      s2 = (13.0/12.0)*(v2-2.0*v3+v4)**2 + 0.25*(v2-v4)**2
      s3 = (13.0/12.0)*(v3-2.0*v4+v5)**2 + 0.25*(3.0*v3-4.0*v4+v5)**2

           a1 = 0.1/(e+s1)**2
           a2 = 0.6/(e+s2)**2
           a3 = 0.3/(e+s3)**2

           w1 = a1/(a1+a2+a3)
           w2 = a2/(a1+a2+a3)
           w3 = a3/(a1+a2+a3)

      dom(ib)%dphi_dzminus(i,j,k+3) = w1*(v1*(1.0/3.0)-v2*(7.0/6.0)+
     &v3*(11.0/6.0)) + w2*(-v2*(1.0/6.0)+v3*(5.0/6.0)+v4*(1.0/3.0))+
     &w3*(v3*(1.0/3.0)+v4*(5.0/6.0)-v5*(1.0/6.0))

        end do
       end do
      end do

      end do

      return
      end subroutine hj_weno_dzminus_3d
!#####################################################################
      subroutine heaviside
!######################################################################
      use vars
      use lsm
      use mpi
      use multidata

      implicit none
      integer :: i,j,k,n_epsl,ib,tti,ttj,ttk
      real :: epsl
      real, parameter :: pi = 3.14159265359
!
! Define an infinitely differentiable smoothed heaviside function h_phi
!
      n_epsl = 2

      do ib=1,nbp

       epsl = n_epsl*max(dom(ib)%dx,dom(ib)%dy,dom(ib)%dz)

        do k=dom(ib)%ksp-pl,dom(ib)%kep+pl
         do i=dom(ib)%isp-pl,dom(ib)%iep+pl
          do j=dom(ib)%jsp-pl,dom(ib)%jep+pl

          if (dom(ib)%phi(i,j,k).lt.(-1.0*epsl))
     &  dom(ib)%h_phi(i,j,k) = 0.0  ! h_phi=0 above free surface

          if (dom(ib)%phi(i,j,k).gt.(epsl)) dom(ib)%h_phi(i,j,k) = 1.0  ! h_phi=1.0 below free surface

          if (abs(dom(ib)%phi(i,j,k)).le.epsl) then
!
! Transition zone across free surface (2 grid cells width either side)
!
          dom(ib)%h_phi(i,j,k)=0.5*(1.0+dom(ib)%phi(i,j,k)/
     &epsl+1.0/pi*sin(pi*dom(ib)%phi(i,j,k)/epsl))

          end if

        dom(ib)%dens(i,j,k)=densg+(densl-densg)*dom(ib)%h_phi(i,j,k) !dens = densg above free surface, densl below
        dom(ib)%mu(i,j,k)=mug+(mul-mug)*dom(ib)%h_phi(i,j,k)  !mu = mug above free surface, mul below

          end do
         end do

        end do

      end do

      call bound_lsm(18) !(dens)
      call bound_lsm(19) !(mu)

      return
      end subroutine heaviside
!#####################################################################

ext: for


      subroutine res_control (jres)

      use basin_module
      use reservoir_data_module
      use time_module
      use reservoir_module
      use climate_module
      use hydrograph_module
      use conditional_module
      use water_body_module

      implicit none

      integer :: ii                   !none          |counter
      integer :: jres                 !none          |reservoir number
      integer :: idat                 !              |
      integer :: ihyd                 !none          |counter
      integer :: ised                 !none          |counter
      integer :: irel                 !              |
      integer :: inut                 !none          |counter
      integer :: iob                  !none          |counter
      real :: pvol_m3
      real :: evol_m3

      iob = res_ob(jres)%ob

      !! set water body pointer to res
      wbody => res(jres)
      wbody_wb => res_wat_d(jres)

      ht1 = ob(icmd)%hin    !! set incoming flow
      ht2 = resz            !! zero outgoing flow

      !! add incoming flow to reservoir
      res(jres) = res(jres) + ht1

      if (time%yrc > res_hyd(jres)%iyres .or. (time%mo >= res_hyd(jres)%mores   &
                                   .and. time%yrc == res_hyd(jres)%iyres)) then
        !! perform reservoir water/sediment balance
        idat = res_ob(jres)%props
        ihyd = res_dat(idat)%hyd
        ised = res_dat(idat)%sed
        if(time%step == 0) then
          !! determine reservoir outflow
          irel = res_dat(idat)%release
          d_tbl => dtbl_res(irel)
          pvol_m3 = res_ob(jres)%pvol
          evol_m3 = res_ob(jres)%evol
          call conditions (ihyd)
          call res_hydro (jres, irel, ihyd, pvol_m3, evol_m3)
          call res_sediment (jres, ihyd, ised)
	    else
	      !call res_hourly
        endif

      !! calculate water balance for day
      iwst = ob(iob)%wst
      res_wat_d(jres)%evap = 10. * res_hyd(ihyd)%evrsv * wst(iwst)%weat%pet * res_wat_d(jres)%area_ha
      res_wat_d(jres)%seep = 240. * res_hyd(ihyd)%k * res_wat_d(jres)%area_ha
      res_wat_d(jres)%precip = 10. * wst(iwst)%weat%precip * res_wat_d(jres)%area_ha

      !! add precip to reservoir storage
      res(jres)%flo = res(jres)%flo + res_wat_d(jres)%precip

      !! subtract outflow from reservoir storage
      res(jres)%flo = res(jres)%flo - ht2%flo
      if (res(jres)%flo < 0.) then
        ht2%flo = ht2%flo + res(jres)%flo
        res(jres)%flo = 0.
      end if

      !! subtract evaporation from reservoir storage
      res(jres)%flo = res(jres)%flo - res_wat_d(jres)%evap
      if (res(jres)%flo < 0.) then
        res_wat_d(jres)%evap = res_wat_d(jres)%evap + res(jres)%flo
        res(jres)%flo = 0.
      end if

      !! subtract seepage from reservoir storage
      res(jres)%flo = res(jres)%flo - res_wat_d(jres)%seep
      if (res(jres)%flo < 0.) then
        res_wat_d(jres)%seep = res_wat_d(jres)%seep + res(jres)%flo
        res(jres)%flo = 0.
      end if

        !! update surface area
        if (res(jres)%flo > 0.) then
          res_wat_d(jres)%area_ha = res_ob(jres)%br1 * res(jres)%flo ** res_ob(jres)%br2
        else
          res_wat_d(jres)%area_ha = 0.
        end if

        !! perform reservoir nutrient balance
        inut = res_dat(idat)%nut
        call res_nutrient (jres, inut, iob)

        !! perform reservoir pesticide transformations
        call res_pest (jres)

        !! set values for outflow variables
        ob(icmd)%hd(1) = ht2

        if (time%step > 0) then
          do ii = 1, time%step
            ob(icmd)%ts(1,ii) = ht2 / real(time%step)
          end do
        end if

        !! set inflow and outflow variables for reservoir_output
        if (time%yrs > pco%nyskip) then
          res_in_d(jres) = ht1
          res_out_d(jres) = ht2
          res_in_d(jres)%flo = res(jres)%flo / 10000.               !m^3 -> ha-m
          res_out_d(jres)%flo = res(jres)%flo / 10000.              !m^3 -> ha-m
          res_wat_d(jres)%evap = res_wat_d(jres)%evap / 10000.      !m^3 -> ha-m
          res_wat_d(jres)%seep = res_wat_d(jres)%seep / 10000.      !m^3 -> ha-m
          res_wat_d(jres)%precip = res_wat_d(jres)%precip / 10000.    !m^3 -> ha-m
        end if

      else
        !! reservoir has not been constructed yet
        ob(icmd)%hd(1) = ob(icmd)%hin
      end if

      return
      end subroutine res_control


text:
module superfluid_def
    use const_def, only: dp

    ! computes superfluid critical temperatures as a function of density
    ! these come courtesy of Dany Page

    ! types. these are pointers to locations in sf_tables
    integer, parameter :: undefined_gap = -1
    integer, parameter :: proton_1S0 = 1
    integer, parameter :: neutron_1S0 = 2
    integer, parameter :: neutron_3P2 = 3
    integer, parameter :: max_number_sf_types = 3

    type sf_table_type
        integer :: which_gap
        logical :: is_loaded
        integer :: nv
        real(dp) :: kF_min, kF_max
        real(dp), dimension(:), allocatable :: kF
        real(dp), dimension(:), pointer :: f=>null()        ! for storing the interpolation coefficients
    end type sf_table_type

    type (sf_table_type), dimension(max_number_sf_types), target :: sf_tables

    logical :: sf_is_initialized = .false.
    character(len=256), save :: sf_datadir

    ! you might want to scale the gaps by some factor; this is normally set to one
    real(dp), dimension(max_number_sf_types), save :: sf_scale

end module superfluid_def

ext: f

max_stars_repo_name: LBJ-Wade/dStar

max_issues_repo_name: LBJ-Wade/dStar

max_forks_repo_name: LBJ-Wade/dStar


text:
      subroutine DSVA(A,MDA,M,N,MDATA,B,SING,KPVEC,NAMES,ISCALE,D,WORK)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
c>> 1996-05-10 DSVA Krogh  Only "*" output on C version.
c>> 1995-11-15 DSVA Krogh  Moved formats up for C conversion.
C>> 1994-10-20 DSVA Krogh  Changes to use M77CON
C>> 1992-04-23 DSVA CLL Replaced Hollerith with '...' in a FORMAT.
C>> 1992-03-18 CLL New treatment of WIDTH and length of labeling names.
C>> 1989-11-21 CLL Fixed bug that caused divide by 0 when M < N.
C>> 1989-03-07 DSVA   CLL Changed logical PRINT to integer vector KPVEC.
C>> 1988-12-01 DSVA   CLL Added PRINT into arg list.
C>> 1987-11-24 DSVA   Lawson  Initial code.
c
C     Singular Value Analysis.  Computes the singular value
c     decomposition of the matrix of a least squares problem, and
c     produces a printed report.
c     ------------------------------------------------------------------
c                         Subroutine Arguments
c     A(,)     [inout]  On entry, contains the M x N matrix of the least
c              squares problem to be analyzed.  This could be a matrix
c              obtained by preliminary orthogonal transformations
c              applied to the actual problem matrix which may have had
c              more rows (See MDATA below.)
c
c     MDA   [in]  First dimensioning parameter for A(,).  Require
c              MDA .ge. max(M, N).
c
c     M,N   [in]  No. of rows and columns, respectively, in the
c              matrix, A.  Either M > N or M .le. N is permitted.
c              Require M > 0 and N > 0.
c
c     MDATA [in]  No. of rows in actual least squares problem.
c              Generally MDATA .ge. M.  MDATA is used only in computing
c              statistics for the report and is not used as a loop
c              count or array dimension.
c
c     B()   [inout]  On entry, contains the right-side vector, b, of the
c              least squares problem.  This vector is of length, M.
c              On return, contains the vector, g = (U**t)*b, where U
c              comes from the singular value decomposition of A.  The
c              vector , g, is also of length M.
c
c     SING()   [out] On return, contains the singular values of A, in
c              descending order, in locations indexed 1 thru min(M,N).
c
c     KPVEC()  [integer array, in]  Array of integers to select print
c              options.  KPVEC(1) determines whether the rest of
c              the array is to be used or ignored.
c              If KPVEC(1) = 1, the contents of (KPVEC(I), I=2,4)
c              will be used to set internal variables as follows:
c                       PRBLK = KPVEC(2)
c                       UNIT  = KPVEC(3)
c                       WIDTH = KPVEC(4)
c              If KPVEC(1) = 0 default settings will be used.  The user
c              need not dimension KPVEC() greater than 1.  The subr will
c              set PRBLK = 111111, UNIT = -1, and WIDTH = 69.
c
c              The internal variables PRBLK, UNIT, and WIDTH are
c              interpreted as follows:
c
c              PRBLK    The decimal representation of PRBLK must be
c              representable as at most 6 digits, each being 0 or 1.
c              The decimal digits will be interpreted as independant
c              on/off flags for the 6 possible blocks of printed output.
c              Examples:  111111 selects all blocks, 0 suppresses all
c              printing,  101010 selects the 1st, 3rd, and 5th blocks,
c              etc.
c              The six blocks are:
c              1. Header, with size and scaling option parameters.
c              2. V-matrix.  Amount of output depends on M and N.
c              3. Singular values and related quantities.  Amount of
c                 output depends on N.
c              4. Listing of YNORM and RNORM and their logarithms.
c                 Amount of output depends on N.
c              5. Levenberg-Marquart analysis.
c              6. Candidate solutions.  Amount of output depends on
c                 M and N.
c
c              UNIT     Selects the output unit.  If UNIT .ge. 0,
c              UNIT will be used as the output unit number.
c              If UNIT = -1, output will be written to the "*" output
c              unit, i.e., the standard system output unit.
c              The calling program unit is responsible for opening
c              and/or closing the selected output unit if the host
c              system requires these actions.
c
c        WIDTH    Default value is 79.  Determines the width of
c        blocks 2, 3, and 6 of the output report.
c        Block 3 will use 95(+1) cols if WIDTH .ge. 95, and otherwise
c        69(+1) cols.
c        Blocks 2 and 6 are printed by subroutine DPRTSV.  These blocks
c        generally use at most WIDTH(+1) cols, but will use more if
c        the names are so long that more space is needed to print one
c        name and one numeric column.  The (+1)'s above are reminders
c        that in all cases there is one extra initial column for Fortran
c        "carriage control".  The carriage control character will always
c        be a blank.
c        Blocks 1, 4, and 5 have fixed widths of 63(+1), 66(+1) and
c        66(+1), respectively.
c
c     NAMES()  [in]  NAMES(j), for j = 1, ..., N, may contain a
c              name for the jth component of the solution
c              vector.  The declared length of the elements of the
c              NAMES() array is not specifically limited, but a
c              greater length reduces the space available for columns
c              of the matris to be printed.
c              If NAMES(1) contains only blank characters,
c              it will be assumed that no names have been provided,
c              and this subr will not access the NAMES() array beyond
c              the first element.
C
C     ISCALE   [in]  Set by the user to 1, 2, or 3 to select the column
c              scaling option.
C                  1   SUBR WILL USE IDENTITY SCALING AND IGNORE THE D()
C                      ARRAY.
C                  2   SUBR WILL SCALE NONZERO COLS TO HAVE UNIT EUCLID-
C                      EAN LENGTH AND WILL STORE RECIPROCAL LENGTHS OF
C                      ORIGINAL NONZERO COLS IN D().
C                  3   USER SUPPLIES COL SCALE FACTORS IN D(). SUBR
C                      WILL MULT COL J BY D(J) AND REMOVE THE SCALING
C                      FROM THE SOLN AT THE END.
c
c     D()   [ignored or out or in]  Usage of D() depends on ISCALE as
c              described above.  When used, its length must be
c              at least N.
c
c     WORK()   [scratch]   Work space of length at least 2*N.  Used
c              directly in this subr and also in _SVDRS.
c     ------------------------------------------------------------------
c--D replaces "?": ?SVA, ?SVDRS, ?PRTSV
c     ------------------------------------------------------------------
c          This code was originally developed by Charles L. Lawson and
c     Richard J. Hanson at Jet Propulsion Laboratory in 1973.  The
c     original code was described and listed in the book,
c
c                  Solving Least Squares Problems
c                  C. L. Lawson and R. J. Hanson
c                  Prentice-Hall, 1974
c
c     Feb 1985, Mar 1987, C. L. Lawson & S. Y. Chiu, JPL.
c     Adapted code from the Lawson & Hanson book to Fortran 77 for use
c     in the JPL MATH77 library.
c     Prefixing subprogram names with S or D for s.p. or d.p. versions.
c     Using generic names for intrinsic functions.
c     Adding calls to BLAS and MATH77 error processing subrs in some
c     program units.
c     Dec 1988 CLL.  For MATH77 this subr was intended to have PRINT in
c     the arg list, and was described that way in the MATH77 manual, but
c     PRINT was omitted in the code.  Now fixing this at Rel. 2.2 by
c     adding PRINT to the arg list.
c     Dec 1988 CLL.  Also for MATH77 the WORK() array was introduced.
c     Previous version used a long SING() array for work space.  WORK()
c     was introduced incorrectly in the S.P. version.  Now correcting
c     this and using the processor CHGTYP to keep the S.P. and D.P.
c     versions similar.
c     1989-03-07 CLL>  Changing the logical PRINT to an integer vector
c     KPVEC().  Giving user control of quantity, destination, and
c     line width of output.
c     ------------------------------------------------------------------
      integer I, IE, IPASS, ISCALE, J, K, KPVEC(4), M, MDA, MDATA
      integer MINMN, MINMN1, MPASS, N, NSOL
      integer PRBLK, UNIT, WIDTH
      double precision A(MDA,N), A1, A2, A3, A4, ALAMB, ALN10, B(M)
      double precision D(N), DEL, DENOM, EL, EL2
      double precision ONE, PCOEF, RL, RNORM, RS
      double precision SB, SING(N), SL, TEMP, TEN, THOU, TWENTY
      double precision WORK(2*N), YL, YNORM, YS, YSQ, ZERO
      character*(*) NAMES(N)
      logical BLK(6), NARROW
c++ CODE for ~.C. is active
      logical STAR
c++ END
C
      parameter( ZERO = 0.0D0, ONE = 1.0D0)
      parameter( TEN = 10.0D0, TWENTY = 20.0D0, THOU = 1000.0D0)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  220 FORMAT (1X/' INDEX  SING. VAL.    P COEF   ',
     *        '   RECIPROCAL    G COEF         G**2   ',
     *        '   CUMULATIVE  SCALED SQRT'/
     *    31x,'   SING. VAL.',26x,
     *        '  SUM of SQRS  of CUM.S.S.')
  221 FORMAT (1X/' INDEX  SING. VAL.    P COEF   ',
     *        '   RECIPROCAL    G COEF     SCALED SQRT'/
     *    31x,'   SING. VAL.',13x,'  of CUM.S.S.')
  222 FORMAT (1X/' INDEX  SING. VAL.    G COEF         G**2   ',
     *        '   CUMULATIVE  SCALED SQRT'/
     *    44x,'  SUM of SQRS  of CUM.S.S.')

  230 FORMAT (' ',4X,'0',64X,2g13.4)
  231 FORMAT (' ',4X,'0',51X, g13.4)
  232 FORMAT (' ',4X,'0',38X,2g13.4)

  240 FORMAT (' ',i5,g12.4,6g13.4)

  260 FORMAT (1X,' M  = ',I6,',   N  =',I4,',   MDATA  =',I8)
  270 FORMAT (1X/' Singular Value Analysis of the least squares',
     *        ' problem,  A*X = B,'/
     *        ' scaled as (A*D)*Y = B.')
  280 FORMAT (1X/' Scaling option No.',I2,'.  D is a diagonal',
     * ' matrix with the following diagonal elements..'/(5X,10E12.4))
  290 FORMAT (1X/' Scaling option No. 1.   D is the identity matrix.'/
     * 1X)
  300 FORMAT (1X/' INDEX',12X,'YNORM      RNORM',11X,
     * '      LOG10      LOG10'/
     * 45X,'      YNORM      RNORM'/1X)
  310 FORMAT (' ',I5,6X,2E11.3,11X,2F11.3)
  320 FORMAT (1X/
     *' Norms of solution and residual vectors for a range of values'/
     *' of the Levenberg-Marquardt parameter, LAMBDA.'//
     *        '      LAMBDA      YNORM      RNORM',
     *         '      LOG10      LOG10      LOG10'/
     *     34X,'     LAMBDA      YNORM      RNORM')
  330 FORMAT (1X, 3E11.3, 3F11.3)
c     ------------------------------------------------------------------
      IF (M.LE.0 .OR. N.LE.0) RETURN
      MINMN = MIN(M,N)
      MINMN1 = MINMN + 1
      if(KPVEC(1) .eq. 0) then
         PRBLK = 111111
         UNIT = -1
         WIDTH = 79
      else
         PRBLK = KPVEC(2)
         UNIT = KPVEC(3)
         WIDTH = KPVEC(4)
      endif
c++ CODE for ~.C. is active
      STAR = UNIT .lt. 0
c++ END
c                           Build logical array BLK() by testing
c                           decimal digits of PRBLK.
      do 20 I=6, 1, -1
         J = PRBLK/10
         BLK(I) = (PRBLK - 10*J) .gt. 0
         PRBLK = J
   20 continue
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c                         Optionally print header and M, N, MDATA
      if(BLK(1)) then
c++ CODE for ~.C. is active
         if(STAR) then
c++ END
            write (*,270)
            write (*,260) M,N,MDATA
c++ CODE for ~.C. is active
         else
            write (UNIT,270)
            write (UNIT,260) M,N,MDATA
         endif
c++ END
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c                                  Handle scaling as selected by ISCALE.
      if( ISCALE .eq. 1) then
         if(BLK(1)) then
c++ CODE for ~.C. is active
            if(STAR) then
c++ END
               write (*,290)
c++ CODE for ~.C. is active
            else
               write (UNIT,290)
            endif
c++ END
         endif
      else
C
C                                  Apply column scaling to A.
C
         DO 52 J = 1,N
            A1 = D(J)
            if( ISCALE .le. 2) then
               SB = ZERO
               DO 30 I = 1,M
   30             SB = SB + A(I,J)**2
               A1 = sqrt(SB)
               IF (A1.EQ.ZERO) A1 = ONE
               A1 = ONE/A1
               D(J) = A1
            endif
            DO 50 I = 1,M
               A(I,J) = A(I,J)*A1
   50       continue
   52    continue
         if(BLK(1)) then
c++ CODE for ~.C. is active
            if(STAR) then
c++ END
               write (*,280) ISCALE,(D(J),J = 1,N)
c++ CODE for ~.C. is active
            else
               write (UNIT,280) ISCALE,(D(J),J = 1,N)
            endif
c++ END
         endif
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C         Compute the Singular Value Decomposition of the scaled matrix.
C
      call DSVDRS (A,MDA,M,N,B,M,1,SING,WORK)
c
c                                               Determine NSOL.
      NSOL = MINMN
      do 60 J = 1,MINMN
         if(SING(J) .eq. ZERO) then
            NSOL = J-1
            go to 65
         endif
   60 continue
   65 continue
C
c              The array B() contains the vector G.
C              Compute cumulative sums of squares of components of
C              G and store them in WORK(I), I = 1,...,MINMN+1
C
      SB = ZERO
      DO 70 I = MINMN1,M
         SB = SB + B(I)**2
   70 CONTINUE
      WORK(MINMN+1) = SB
      DO  75 J = MINMN, 1, -1
         SB = SB + B(J)**2
         WORK(J) = SB
   75 CONTINUE
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C                                          PRINT THE V MATRIX.
C
      if(BLK(2)) CALL DPRTSV (A,MDA,N,N,NAMES,1, UNIT, WIDTH)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C                                REPLACE V BY D*V IN THE ARRAY A()
      if (ISCALE .gt.1) then
         do 82 I = 1,N
            do 80 J = 1,N
               A(I,J) = D(I)*A(I,J)
   80       continue
   82    continue
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(BLK(3)) then
c                     Print singular values and other summary results.
c
c     Output will be done using one of two layouts.  The narrow
c     layout uses 69 cols + 1 for carriage control, and makes two passes
c     through the computation.
c     The wide layout uses 95 cols + 1 for carriage control, and makes
c     only one pass through the computation.
c
C                       G  NOW IN  B() ARRAY.  V NOW IN A(,) ARRAY.
C
      NARROW = WIDTH .lt. 95
      MPASS = 1
      if(NARROW) MPASS = 2
      do 170 IPASS = 1, MPASS
c++ CODE for ~.C. is active
         if(STAR) then
c++ END
            if(NARROW) then
               if(IPASS .eq. 1) then
                  write(*,221)
               else
                  write(*,222)
               endif
            else
               write (*,220)
            endif
c++ CODE for ~.C. is active
         else
            if(NARROW) then
               if(IPASS .eq. 1) then
                  write(UNIT,221)
               else
                  write(UNIT,222)
               endif
            else
               write (UNIT,220)
            endif
         endif
c++ END
c                               The following stmt converts from
c                               integer to floating-point.
      DENOM = max(1,MDATA)
      A3 = WORK(1)
      A4 = sqrt(A3/ DENOM  )
c++ CODE for ~.C. is active
      if(STAR) then
c++ END
         if(NARROW) then
            if(IPASS .eq. 1) then
               write(*,231) A4
            else
               write(*,232) A3, A4
            endif
         else
            write (*,230) A3,A4
         endif
c++ CODE for ~.C. is active
      else
         if(NARROW) then
            if(IPASS .eq. 1) then
               write(UNIT,231) A4
            else
               write(UNIT,232) A3, A4
            endif
         else
            write (UNIT,230) A3,A4
         endif
      endif
c++ END
C
      DO 160 K = 1,MINMN
         if (SING(K).EQ.ZERO) then
            PCOEF = ZERO
c++ CODE for ~.C. is active
            if(STAR) then
c++ END
               write (*,240) K,SING(K)
c++ CODE for ~.C. is active
            else
               write (UNIT,240) K,SING(K)
            endif
c++ END
         else
            PCOEF = B(K) / SING(K)
            A1 = ONE / SING(K)
            A2 = B(K)**2
            A3 = WORK(K+1)
c                                 The following stmt converts from
c                                 integer to floating-point.
            TEMP = MAX(1,MDATA-K)
            A4 = sqrt(A3/TEMP )
c++ CODE for ~.C. is active
            if(STAR) then
c++ END
               if(NARROW) then
                  if(IPASS .eq. 1) then
                     write(*,240) K,SING(K),PCOEF,A1,B(K),      A4
                  else
                     write(*,240) K,SING(K),         B(K),A2,A3,A4
                  endif
               else
                  write (*,240) K,SING(K),PCOEF,A1,B(K),A2,A3,A4
               endif
c++ CODE for ~.C. is active
            else
               if(NARROW) then
                  if(IPASS .eq. 1) then
                     write(UNIT,240) K,SING(K),PCOEF,A1,B(K),      A4
                  else
                     write(UNIT,240) K,SING(K),         B(K),A2,A3,A4
                  endif
               else
                  write (UNIT,240) K,SING(K),PCOEF,A1,B(K),A2,A3,A4
               endif

            endif
c++ END
         endif
  160 continue
  170 continue
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if( BLK(4) ) then
C
C         Compute and print values of YNORM, RNORM and their logarithms.
C
c++ CODE for ~.C. is active
      if(STAR) then
c++ END
         write (*,300)
c++ CODE for ~.C. is active
      else
         write (UNIT,300)
      endif
c++ END
      YSQ = ZERO
      do 180 J = 0, NSOL
         if(J .ne. 0) YSQ = YSQ + (B(J) / SING(J))**2
         YNORM = sqrt(YSQ)
      RNORM = sqrt(WORK(J+1))
      YL = -THOU
      IF (YNORM .GT. ZERO) YL = LOG10(YNORM)
      RL = -THOU
      IF (RNORM .GT. ZERO) RL = LOG10(RNORM)
c++ CODE for ~.C. is active
      if(STAR) then
c++ END
         write (*,310) J,YNORM,RNORM,YL,RL
c++ CODE for ~.C. is active
      else
         write (UNIT,310) J,YNORM,RNORM,YL,RL
      endif
c++ END
  180 continue
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if( BLK(5) .and. SING(1) .ne. ZERO ) then
C
C     COMPUTE VALUES OF XNORM AND RNORM FOR A SEQUENCE OF VALUES OF
C     THE LEVENBERG-MARQUARDT PARAMETER.
C
         EL = LOG10(SING(1)) + ONE
         EL2 = LOG10(SING(NSOL)) - ONE
         DEL = (EL2-EL) / TWENTY
         ALN10 = LOG(TEN)
c++ CODE for ~.C. is active
         if(STAR) then
c++ END
            write (*,320)
c++ CODE for ~.C. is active
         else
            write (UNIT,320)
         endif
c++ END
         DO 200 IE = 1,21
C                                    COMPUTE        ALAMB = 10.0**EL
            ALAMB = EXP(ALN10*EL)
            YS = ZERO
            RS = WORK(NSOL+1)
            DO 190 I = 1,MINMN
               SL = SING(I)**2 + ALAMB**2
               YS = YS + (B(I)*SING(I)/SL)**2
               RS = RS + (B(I)*(ALAMB**2)/SL)**2
  190       CONTINUE
            YNORM = sqrt(YS)
            RNORM = sqrt(RS)
            RL = -THOU
            IF (RNORM.GT.ZERO) RL = LOG10(RNORM)
            YL = -THOU
            IF (YNORM.GT.ZERO) YL = LOG10(YNORM)
c++ CODE for ~.C. is active
            if(STAR) then
c++ END
               write (*,330) ALAMB,YNORM,RNORM,EL,YL,RL
c++ CODE for ~.C. is active
            else
               write (UNIT,330) ALAMB,YNORM,RNORM,EL,YL,RL
            endif
c++ END
            EL = EL + DEL
  200    CONTINUE
      endif
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C                   Compute and optionally print candidate solutions.
C
      do 215 K = 1,NSOL
         PCOEF = B(K) / SING(K)
         DO 210 I = 1,N
                A(I,K) = A(I,K) * PCOEF
  210           IF (K.GT.1) A(I,K) = A(I,K) + A(I,K-1)
  215 continue
      if (BLK(6) .and. NSOL.GE.1)
     *       CALL DPRTSV (A,MDA,N,NSOL,NAMES,2,UNIT,WIDTH)
      return
      END


text:
c  bpsplinb -- dmc 30 May 1996
c
c  set up coefficients for bicubic spline with following BC's:
c  * LHS and RHS BC under user control (see comments)
c  * derivatives periodic in second coordinate (use pspline.for)
c
      subroutine r8bpsplinb(x,inx,th,inth,fspl,inf3,
     >                    ibcxmin,bcxmin,ibcxmax,bcxmax,
     >                    wk,nwk,ilinx,ilinth,ier)
c
      IMPLICIT NONE
      INTEGER, PARAMETER :: R8=SELECTED_REAL_KIND(12,100)
      integer inx,inth,inf3,nwk
      REAL*8 x(inx),th(inth),fspl(4,4,inf3,inth),wk(nwk)
      integer ibcxmin,ibcxmax
      REAL*8 bcxmin(inth),bcxmax(inth)
c
c  input:
c    x(1...inx) -- abscissae, first dimension of data
c   th(1...inth) -- abscissae, second (periodic) dimension of data
c   fspl(1,1,1..inx,1..inth) -- function values
c   inf3 -- fspl dimensioning, inf3.ge.inx required.
c
c  boundary conditions input:
c   ibcxmin -- indicator for boundary condition at x(1):
c    bcxmin(...) -- boundary condition data
c     =0 -- use "not a knot", bcxmin(...) ignored
c     =1 -- match slope, specified at x(1),th(ith) by bcxmin(ith)
c     =2 -- match 2nd derivative, specified at x(1),th(ith) by bcxmin(ith)
c     =3 -- boundary condition is slope=0 (df/dx=0) at x(1), all th(j)
c     =4 -- boundary condition is d2f/dx2=0 at x(1), all th(j)
c     =5 -- match 1st derivative to 1st divided difference
c     =6 -- match 2nd derivative to 2nd divided difference
c     =7 -- match 3rd derivative to 3rd divided difference
c           (for more detailed definition of BCs 5-7, see the
c           comments of subroutine mkspline)
c   NOTE bcxmin(...) referenced ONLY if ibcxmin=1 or ibcxmin=2
c
c   ibcxmax -- indicator for boundary condition at x(nx):
c    bcxmax(...) -- boundary condition data
c     (interpolation as with ibcxmin, bcxmin)
c
c  output:
c   fspl(*,*,1..inx,1..inth) -- bicubic spline coeffs (4x4)
c   ...fspl(1,1,*,*) is not replaced.
c
      integer ilinx,ilinth,ier
c
c   ilinx -- =1 on output if x(inx) pts are nearly evenly spaced (tol=1e-3)
c   ilinth-- =1 on output if th(inth) evenly spaced (tol=1e-3)
c
c   ier -- completion code, 0 for normal
c
c  workspace:
c   wk -- must be at least 5*max(inx,inth) large
c  nwk -- size of workspace
c
c---------------------------------
c  compute bicubic spline of 2d function, given values at the grid
c  grid crossing points, f(1,1,i,j)=f(x(i),th(j)).
c
c  on evaluation:  for point x btw x(i) and x(i+1), dx=x-x(i)
c                       and th btw th(j) and th(j+1), dt=th-th(j),
c
c      spline =
c        f(1,1,i,j) + dx*f(2,1,i,j) + dx**2*f(3,1,i,j) + dx**3*f(4,1,i,j)
c   +dt*(f(1,2,i,j) + dx*f(2,2,i,j) + dx**2*f(3,2,i,j) + dx**3*f(4,2,i,j))
c   +d2*(f(1,3,i,j) + dx*f(2,3,i,j) + dx**2*f(3,3,i,j) + dx**3*f(4,3,i,j))
c   +d3*(f(1,4,i,j) + dx*f(2,4,i,j) + dx**2*f(3,4,i,j) + dx**3*f(4,4,i,j))
c
c      where d2=dt**2 and d3=dt**3.
c
c---------------------------------
c
      integer ierx,ierth
      REAL*8 zdum(1)
c
c---------------------------------
c
      ier=0
      if(nwk.lt.5*max(inx,inth)) then
         write(6,'('' ?bpsplinb:  workspace too small.'')')
         ier=1
      endif
      if(inx.lt.2) then
         write(6,'('' ?bpsplinb:  at least 2 x points required.'')')
         ier=1
      endif
      if(inth.lt.2) then
         write(6,'('' ?bpsplinb:  need at least 2 theta points.'')')
         ier=1
      endif
c
      call ibc_ck(ibcxmin,'bcspline','xmin',0,7,ier)
      call ibc_ck(ibcxmax,'bcspline','xmax',0,7,ier)
c
c  check ilinx & x vector
c
      call r8splinck(x,inx,ilinx,1.0E-3_r8,ierx)
      if(ierx.ne.0) ier=2
c
      if(ier.eq.2) then
         write(6,'('' ?bpsplinb:  x axis not strict ascending'')')
      endif
c
c  check ilinth & th vector
c
      call r8splinck(th,inth,ilinth,1.0E-3_r8,ierth)
      if(ierth.ne.0) ier=3
c
      if(ier.eq.3) then
         write(6,'('' ?bpsplinb:  th axis not strict ascending'')')
      endif
c
      if(ier.ne.0) return
c
c------------------------------------
      zdum=0.0_r8
c
      call r8bcspline(x,inx,th,inth,fspl,inf3,
     >   ibcxmin,bcxmin,ibcxmax,bcxmax,
     >   -1,zdum,-1,zdum,
     >   wk,nwk,ilinx,ilinth,ier)
c
      end

max_stars_repo_name: ORNL-Fusion/RFSciDAC-testing

max_issues_repo_name: ORNL-Fusion/RFSciDAC-testing

max_forks_repo_name: ORNL-Fusion/RFSciDAC-testing


text:
subroutine AspirationUpdate(nodeUpdate, payoffUpdate, randUpdate)
    use global
    implicit none
    integer,intent(in)::nodeUpdate
    real,intent(in)::payoffUpdate, randUpdate

    real::switchProb
    real::AspirationMinusPayoff
    !!

    AspirationMinusPayoff = aspirationValueNode(nodeUpdate) - payoffUpdate
    switchProb = 1.0 / (1 + exp(-1.0 * selectionIntensity * AspirationMinusPayoff))
    !call RANDOM_NUMBER(randUpdate)
    !call RNUN(randUpdateArray)
    !randUpdate = randUpdateArray(1)

    if(randUpdate < switchProb)then
        strategyPGG(nodeUpdate) = 1 - strategyPGG(nodeUpdate)
    end if

    !write(*,*) "Aspiration Value"
    !write(*,*) aspirationValueNode(nodeUpdate)
    !write(*,*) "Payoff Update Node"
    !write(*,*) payoffUpdate
    !write(*,*) "randUpdate"
    !write(*,*) randUpdate
    !write(*,*) "Switching Probability"
    !write(*,*) switchProb
    !pause
    return
end subroutine AspirationUpdate

max_stars_repo_name: leizhougetbetter/IAD

max_issues_repo_name: leizhougetbetter/IAD

max_forks_repo_name: leizhougetbetter/IAD