!
!  Copyright 2019-2020 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130

! for single-scale Maxwell-TDDFT method
module fdtd_coulomb_gauge
  implicit none

contains

subroutine fdtd_singlescale(itt,lg,mg,system,info,rho,Vh,j_e,srg_scalar,Ac,div_Ac,fw)
  use structures
  use math_constants,only : zi,pi
  use phys_constants, only: cspeed_au
  use salmon_global, only: dt,method_singlescale,yn_symmetrized_stencil
  use sendrecv_grid, only: update_overlap_real8
  use stencil_sub, only: calc_gradient_field
  use communication, only: comm_is_root, comm_summation
  use inputoutput, only: t_unit_time
  use timer
  implicit none
  integer                 ,intent(in) :: itt
  type(s_rgrid)           ,intent(in) :: lg,mg
  type(s_dft_system)      ,intent(in) :: system
  type(s_parallel_info)   ,intent(in) :: info
  type(s_scalar)          ,intent(in) :: rho,Vh ! electron number density & Hartree potential
  type(s_vector)          ,intent(in) :: j_e    ! electron number current density (without rho*A/c)
  type(s_sendrecv_grid)               :: srg_scalar
  type(s_vector)                      :: Ac     ! A/c, A: vector potential, c: speed of light
  type(s_scalar)                      :: div_Ac ! div(A/c)
  type(s_singlescale)                 :: fw     ! FDTD working arrays, etc.
  !
  integer,parameter :: mstep=100
  integer :: ix,iy,iz,i1,ii,krd(3,3),lcs(3,3,3),dr(3)
  real(8) :: Hvol,hgs(3),dt_m,tm,coef,lap_A,Energy_em,diff_A,coef2 &
  & ,e_em,e_em_wrk,e_joule,e_joule_wrk,e_poynting(2),e_poynting_wrk(2),rho_t
  real(8),dimension(3) :: out_Aext,out_Ab1,out_Ab2,wrk,wrk2,wrk3,wrk4,vec_je,Aext0,Aext1,Aext0_old,Aext1_old
  real(8) :: e_poy1,e_poy2,rtmp1(6),rtmp2(6)

  call timer_begin(LOG_SS_FDTD_CALC)

  krd = 0
  krd(1,1) = 1; krd(2,2) = 1; krd(3,3) = 1

  lcs = 0
  lcs(1,2,3) = 1; lcs(3,1,2) = 1; lcs(2,3,1) = 1
  lcs(1,3,2) = -1; lcs(2,1,3) = -1; lcs(3,2,1) = -1

  hgs = system%hgs
  Hvol = system%Hvol
  dt_m = dt / dble(mstep)

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

  fw%box = 0d0
  wrk = 0d0
!$OMP parallel do collapse(2) private(ix,iy,iz,vec_je,rho_t) reduction(+:wrk)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    vec_je = ( j_e%v(1:3,ix,iy,iz) + fw%vec_je_old(1:3,ix,iy,iz) )*0.5d0 ! j_e(t) = ( j_e(t+dt/2) + j_e(t-dt/2) )/2
    rho_t  = ( rho%f(ix,iy,iz)     + fw%rho_old(ix,iy,iz)        )*0.5d0 ! rho(t) = ( rho(t+dt/2) + rho(t-dt/2) )/2
    fw%curr(ix,iy,iz,1:3) = vec_je + rho_t * fw%vec_Ac_m(1,ix,iy,iz,1:3) ! curr(t): electron number current density
    wrk = wrk + fw%curr(ix,iy,iz,1:3)

    fw%box(ix,iy,iz) = Vh%f(ix,iy,iz) ! Vh(t+dt/2)
  end do
  end do
  end do
  wrk = wrk/dble(lg%num(1)*lg%num(2)*lg%num(3))
  call timer_end(LOG_SS_FDTD_CALC)

  call timer_begin(LOG_SS_FDTD_COMM_COLL)
  call comm_summation(wrk,fw%curr_ave,3,info%icomm_r)
  call timer_end(LOG_SS_FDTD_COMM_COLL)

! gradient of d(Vh)/dt (Vh: Hartree potential)
  call timer_begin(LOG_SS_FDTD_COMM)
  if(info%if_divide_rspace) call update_overlap_real8(srg_scalar, mg, fw%box1)
  call timer_end(LOG_SS_FDTD_COMM)

  call timer_begin(LOG_SS_FDTD_CALC)
  call calc_gradient_field(mg,fw%coef_nab,system%rmatrix_B,fw%box1,fw%grad_Vh) ! grad[Vh(t+dt/2)]

  !$OMP parallel do collapse(2) private(ix,iy,iz)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    fw%current4pi(ix,iy,iz,1:3) = ( fw%grad_Vh(1:3,ix,iy,iz) - fw%grad_Vh_old(1:3,ix,iy,iz) ) /dt & ! d(grad(Vh))/dt
                                & - 4d0*pi * fw%curr(ix,iy,iz,1:3) ! 4*pi* (electric current density)
  end do
  end do
  end do
  call timer_end(LOG_SS_FDTD_CALC)

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

! FDTD loop: Ac(t) --> Ac(t+dt)

  if(method_singlescale=='3d') then
    call fdtd
  else
    call fdtd_gbp
  end if

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

! divergence & rotation of Ac(t+dt)

  call timer_begin(LOG_SS_FDTD_CALC)
  fw%div_Ac = 0d0
  fw%rot_Ac = 0d0
  do i1=1,3

!$OMP parallel do collapse(2) private(ix,iy,iz)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
      fw%box(ix,iy,iz) = fw%vec_Ac_m(1,ix,iy,iz,i1)
    end do
    end do
    end do
    call timer_begin(LOG_SS_FDTD_COMM)
    call update_overlap_real8(fw%srg_eg, fw%eg, fw%box)
    call timer_end(LOG_SS_FDTD_COMM)
    if(mg%is(3)==lg%is(3))then
!$OMP parallel do collapse(2) private(ix,iy,iz)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        fw%box(ix,iy,lg%is(3)-1) = fw%vec_Ac_boundary_bottom(ix,iy,i1)
      end do
      end do
    end if
    if(mg%ie(3)==lg%ie(3))then
!$OMP parallel do collapse(2) private(ix,iy,iz)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        fw%box(ix,iy,lg%ie(3)+1) = fw%vec_Ac_boundary_top(ix,iy,i1)
      end do
      end do
    end if

  !$OMP parallel do collapse(2) private(ix,iy,iz,wrk,dr,diff_A)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)
    ! rot(A)
      wrk(1) = ( fw%box(ix+1,iy,iz) - fw%box(ix-1,iy,iz) ) / ( 2d0* Hgs(1) )
      wrk(2) = ( fw%box(ix,iy+1,iz) - fw%box(ix,iy-1,iz) ) / ( 2d0* Hgs(2) )
      wrk(3) = ( fw%box(ix,iy,iz+1) - fw%box(ix,iy,iz-1) ) / ( 2d0* Hgs(3) )
      fw%rot_Ac(:,ix,iy,iz) = fw%rot_Ac(:,ix,iy,iz) + lcs(:,1,i1) * wrk(1) + lcs(:,2,i1) * wrk(2) + lcs(:,3,i1) * wrk(3)
    ! div(A)
      dr = krd(:,i1)
      diff_A = ( fw%box(ix+dr(1),iy+dr(2),iz+dr(3)) - fw%box(ix-dr(1),iy-dr(2),iz-dr(3)) ) / ( 2d0* Hgs(i1) )
      fw%div_Ac(ix,iy,iz) = fw%div_Ac(ix,iy,iz) + diff_A
    end do
    end do
    end do

  end do

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

! Ac & div_Ac for TDDFT

!$OMP parallel do collapse(2) private(ix,iy,iz)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)

    Ac%v(1:3,ix,iy,iz) = ( fw%vec_Ac_m(1,ix,iy,iz,1:3) + fw%vec_Ac_old(1:3,ix,iy,iz) ) * 0.5d0 ! Ac(t+dt/2) = ( A(t+dt) + A(t) )/2
    div_Ac%f(ix,iy,iz)   = ( fw%div_Ac(ix,iy,iz) + fw%div_Ac_old(ix,iy,iz) ) * 0.5d0 ! div ( A(t+dt) + A(t) )/2

  end do
  end do
  end do

  if(yn_symmetrized_stencil=='y' .and. info%if_divide_rspace) then
    call update_overlap_real8(srg_scalar, mg, Ac%v(1,:,:,:))
    call update_overlap_real8(srg_scalar, mg, Ac%v(2,:,:,:))
    call update_overlap_real8(srg_scalar, mg, Ac%v(3,:,:,:))
  end if

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

! Poynting vector, energy, etc.

  coef = cspeed_au / (4d0*pi)
  coef2 = Hvol / (8d0*pi)
  e_em_wrk = 0d0 ! for Electro-Magnetic energy
  e_joule_wrk = 0d0 ! for Joule dissipated power

!$OMP parallel do collapse(2) private(ix,iy,iz,wrk,wrk2,wrk3,wrk4) reduction(+:e_em_wrk,e_joule_wrk)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)

    wrk4 = ( fw%vec_Ac_m(1,ix,iy,iz,:) - fw%vec_Ac_old(:,ix,iy,iz) ) / dt ! (A(t+dt)-A(t))/dt
    wrk  = - (-1d0)*fw%grad_Vh(:,ix,iy,iz) - wrk4 ! E
    wrk2 = cspeed_au * fw%rot_Ac(:,ix,iy,iz)      ! B
    wrk3 = - fw%curr(ix,iy,iz,1:3)                ! j
    fw%poynting_vector(:,ix,iy,iz) = coef * ( lcs(:,1,2) * wrk(1) * wrk2(2) + lcs(:,1,3) * wrk(1) * wrk2(3) &
                                            + lcs(:,2,1) * wrk(2) * wrk2(1) + lcs(:,2,3) * wrk(2) * wrk2(3) &
                                            + lcs(:,3,1) * wrk(3) * wrk2(1) + lcs(:,3,2) * wrk(3) * wrk2(2) ) ! E x B
    e_em_wrk = e_em_wrk + coef2 * ( sum(wrk**2) + sum(wrk2**2) ) ! ( E^2 + B^2 )/(8*pi)
    e_joule_wrk = e_joule_wrk + sum(wrk3*wrk) * Hvol ! j*E

  end do
  end do
  end do
  call timer_end(LOG_SS_FDTD_CALC)

  call timer_begin(LOG_SS_FDTD_COMM_COLL)
  call comm_summation(e_em_wrk,e_em,info%icomm_r)
  call comm_summation(e_joule_wrk,e_joule,info%icomm_r)
  call timer_end(LOG_SS_FDTD_COMM_COLL)

  call timer_begin(LOG_SS_FDTD_CALC)

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

  ! integral(A) @ z = 0 (bottom boundary)
  wrk = 0d0
  if (mg%is(3) == lg%is(3)) then
!$omp parallel do private(ix,iy,iz) reduction(+:wrk)
    do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        wrk(1) = wrk(1) + Ac%v(1,ix,iy,mg%is(3))
        wrk(2) = wrk(2) + Ac%v(2,ix,iy,mg%is(3))
        wrk(3) = wrk(3) + Ac%v(3,ix,iy,mg%is(3))
      end do
    end do
  end if

  ! integral(A) @ z = az (top boundary)
  wrk3 = 0d0
  if (mg%ie(3) == lg%ie(3)) then
!$omp parallel do private(ix,iy,iz) reduction(+:wrk3)
    do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        wrk3(1) = wrk3(1) + Ac%v(1,ix,iy,mg%ie(3))
        wrk3(2) = wrk3(2) + Ac%v(2,ix,iy,mg%ie(3))
        wrk3(3) = wrk3(3) + Ac%v(3,ix,iy,mg%ie(3))
      end do
    end do
  end if
  call timer_end(LOG_SS_FDTD_CALC)

  call timer_begin(LOG_SS_FDTD_COMM_COLL)
  rtmp1 = [wrk, wrk3]
  call comm_summation(rtmp1,rtmp2,6,info%icomm_r)
  wrk (1:3) = rtmp2(1:3)
  wrk3(1:3) = rtmp2(4:6)
  call timer_end(LOG_SS_FDTD_COMM_COLL)

  call timer_begin(LOG_SS_FDTD_CALC)
  out_Ab1 = wrk  / dble(lg%num(1)*lg%num(2))
  out_Ab2 = wrk3 / dble(lg%num(1)*lg%num(2))

! Surface integral of the poynting vector S
  e_poy1 = 0d0
  e_poy2 = 0d0
  coef = Hgs(1)*Hgs(2)
  if(mg%is(3)==lg%is(3)) then ! integral(S) @ z = 0 (bottom boundary)
!$omp parallel do private(ix,iy,iz) reduction(+:e_poy1)
    do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        e_poy1 = e_poy1 + fw%poynting_vector(3,ix,iy,lg%is(3)) * coef
      end do
    end do
  end if
  if(mg%ie(3)==lg%ie(3)) then ! integral(S) @ z = az (top boundary)
!$omp parallel do private(ix,iy,iz) reduction(+:e_poy2)
    do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        e_poy2 = e_poy2 + fw%poynting_vector(3,ix,iy,lg%ie(3)) * coef
      end do
    end do
  end if
  call timer_end(LOG_SS_FDTD_CALC)

  call timer_begin(LOG_SS_FDTD_COMM_COLL)
  e_poynting_wrk = [e_poy1, e_poy2]
  call comm_summation(e_poynting_wrk,e_poynting,2,info%icomm_r)
  call timer_end(LOG_SS_FDTD_COMM_COLL)

  call timer_begin(LOG_SS_FDTD_CALC)
  Energy_em = e_em
  fw%Energy_joule = fw%Energy_joule + dt*e_joule
  fw%Energy_poynting = fw%Energy_poynting + dt*e_poynting

  if(comm_is_root(info%id_rko)) write(fw%fh_rt_micro,'(99(1X,E23.15E3))') &
    dble(itt)*dt*t_unit_time%conv,out_Ab1,out_Ab2,out_Aext,fw%curr_ave,fw%E_electron,fw%Energy_poynting,Energy_em,fw%Energy_joule

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

! FIXME: is following result unused?
! for spatial distribution of excitation energy
!  coef = Hgs(1)*Hgs(2)
!  fw%integral_poynting_tmp = 0d0
!!$omp parallel do collapse(2) private(iz,iy,ix)
!  do iy=mg%is(2),mg%ie(2)
!  do ix=mg%is(1),mg%ie(1)
!    do iz=mg%is(3),mg%ie(3)
!      fw%integral_poynting_tmp(iz) = fw%integral_poynting_tmp(iz) + fw%poynting_vector(3,ix,iy,iz) * coef
!    end do
!  end do
!  end do
!  call timer_end(LOG_SS_FDTD_CALC)
!
!  call timer_begin(LOG_SS_FDTD_COMM_COLL)
!  call comm_summation(fw%integral_poynting_tmp,fw%integral_poynting_tmp2,lg_num(3),info%icomm_r)
!  call timer_end(LOG_SS_FDTD_COMM_COLL)
!
!  call timer_begin(LOG_SS_FDTD_CALC)
!  fw%integral_poynting = fw%integral_poynting + dt * fw%integral_poynting_tmp2

!  fw%Ac_zt_t = 0d0
!! for the vector potential Ax(z,t)
!  do iz=mg%is(3),mg%ie(3)
!!$omp parallel do collapse(2) private(iy,ix)
!    do iy=mg%is(2),mg%ie(2)
!    do ix=mg%is(1),mg%ie(1)
!        fw%Ac_zt_t(1,iz) = fw%Ac_zt_t(1,iz) + Ac%v(1,ix,iy,iz) / (lg%num(1)*lg%num(2))
!        fw%Ac_zt_t(2,iz) = fw%Ac_zt_t(2,iz) + Ac%v(2,ix,iy,iz) / (lg%num(1)*lg%num(2))
!        fw%Ac_zt_t(3,iz) = fw%Ac_zt_t(3,iz) + Ac%v(3,ix,iy,iz) / (lg%num(1)*lg%num(2))
!    end do
!    end do
!  end do
!  call timer_end(LOG_SS_FDTD_CALC)
!
!  call timer_begin(LOG_SS_FDTD_COMM_COLL)
!  call comm_summation(fw%Ac_zt_t,fw%Ac_zt,size(fw%Ac_zt),info%icomm_r)
!  call timer_end(LOG_SS_FDTD_COMM_COLL)
!
!  if(comm_is_root(info%id_rko)) then
!    do iz=lg%is(3),lg%ie(3)
!      write(fw%fh_Ac_zt,fmt='(99(1X,E23.15E3))',advance='no') dble(iz)*hgs(3),fw%Ac_zt(1,iz),fw%Ac_zt(2,iz),fw%Ac_zt(3,iz)
!    end do
!    write(fw%fh_Ac_zt,'()')
!  end if

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

! stock old Ac, grad_Vh, j_e, & rho

!$OMP parallel do collapse(2) private(ix,iy,iz)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    fw%vec_Ac_old(:,ix,iy,iz)    = fw%vec_Ac_m(1,ix,iy,iz,1:3) ! Ac(t+dt) --> Ac(t) of next step
    fw%div_Ac_old(ix,iy,iz)      = fw%div_Ac(ix,iy,iz)      ! div Ac(t+dt) --> div Ac(t) of next step
    fw%grad_Vh_old(1:3,ix,iy,iz) = fw%grad_Vh(1:3,ix,iy,iz) ! grad[Vh(t-dt/2)] of next step
    fw%vec_je_old(1:3,ix,iy,iz)  = j_e%v(1:3,ix,iy,iz) ! j_e(t-dt/2) of next step
    fw%rho_old(ix,iy,iz)         = rho%f(ix,iy,iz)     ! rho(t-dt/2) of next step
  end do
  end do
  end do

  call timer_end(LOG_SS_FDTD_CALC)

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

  return

contains

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

  subroutine pulse(t,r,A_ext)
    use em_field, only: calc_Ac_ext
    implicit none
    real(8),intent(in)  :: t,r
    real(8),intent(out) :: A_ext(3)
    !
    real(8) :: tt

    tt = t - r/cspeed_au
    call calc_Ac_ext(tt,A_ext)

    return
  end subroutine pulse

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

  subroutine fdtd
    implicit none

    call timer_begin(LOG_SS_FDTD_CALC)
    do ii=1,mstep

    !$OMP parallel do collapse(2) private(ix,iy,iz)
      do iz=mg%is(3),mg%ie(3)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)

        fw%vec_Ac_m(-1,ix,iy,iz,1:3) = fw%vec_Ac_m(0,ix,iy,iz,1:3)
        fw%vec_Ac_m(0 ,ix,iy,iz,1:3) = fw%vec_Ac_m(1,ix,iy,iz,1:3)

      end do
      end do
      end do

      do i1=1,3

    !$OMP parallel do collapse(2) private(ix,iy,iz)
        do iz=mg%is(3),mg%ie(3)
        do iy=mg%is(2),mg%ie(2)
        do ix=mg%is(1),mg%ie(1)
          fw%box(ix,iy,iz) = fw%vec_Ac_m(0,ix,iy,iz,i1)
        end do
        end do
        end do
        call timer_end(LOG_SS_FDTD_CALC)

        call timer_begin(LOG_SS_FDTD_COMM)
        call update_overlap_real8(fw%srg_eg, fw%eg, fw%box)
        call timer_end(LOG_SS_FDTD_COMM)

        call timer_begin(LOG_SS_FDTD_CALC)
        if(mg%is(3)==lg%is(3))then
    !$OMP parallel do collapse(2) private(ix,iy,iz)
          do iy=mg%is(2),mg%ie(2)
          do ix=mg%is(1),mg%ie(1)
            fw%box(ix,iy,lg%is(3)-1) = fw%vec_Ac_boundary_bottom(ix,iy,i1)
          end do
          end do
        end if
        if(mg%ie(3)==lg%ie(3))then
    !$OMP parallel do collapse(2) private(ix,iy,iz)
          do iy=mg%is(2),mg%ie(2)
          do ix=mg%is(1),mg%ie(1)
            fw%box(ix,iy,lg%ie(3)+1) = fw%vec_Ac_boundary_top(ix,iy,i1)
          end do
          end do
        end if

    !$OMP parallel do collapse(2) private(ix,iy,iz,lap_A)
        do iz=mg%is(3),mg%ie(3)
        do iy=mg%is(2),mg%ie(2)
        do ix=mg%is(1),mg%ie(1)
          lap_A = ( - 2d0* fw%box(ix,iy,iz) + fw%box(ix-1,iy,iz) + fw%box(ix+1,iy,iz) ) / Hgs(1)**2 &
                + ( - 2d0* fw%box(ix,iy,iz) + fw%box(ix,iy-1,iz) + fw%box(ix,iy+1,iz) ) / Hgs(2)**2 &
                + ( - 2d0* fw%box(ix,iy,iz) + fw%box(ix,iy,iz-1) + fw%box(ix,iy,iz+1) ) / Hgs(3)**2
          fw%vec_Ac_m(1,ix,iy,iz,i1) = ( cspeed_au * dt_m )**2 * lap_A &
                                    + 2.d0* fw%box(ix,iy,iz) - fw%vec_Ac_m(-1,ix,iy,iz,i1) &
                                    + dt_m**2 * fw%current4pi(ix,iy,iz,i1)
        end do
        end do
        end do

      end do ! i1 (spacial )

    ! external field
      tm = ( dble(itt-1) + dble(ii-1)/dble(mstep) ) *dt
      call pulse(tm,0d0,   Aext0_old)
      call pulse(tm,Hgs(3),Aext1_old)
      call pulse(tm+dt_m,0d0,   Aext0)
      call pulse(tm+dt_m,Hgs(3),Aext1)
      out_Aext = Aext1

    ! z axis: Mur absorbing boundary condition
      coef = ( cspeed_au * dt_m - Hgs(3) ) / ( cspeed_au * dt_m + Hgs(3) )
      if(mg%is(3)==lg%is(3))then
    !$OMP parallel do collapse(2) private(ix,iy,iz)
        do iy=mg%is(2),mg%ie(2)
        do ix=mg%is(1),mg%ie(1)
        ! absorbing boundary condition with the incident field vec_Ac_ext
          fw%vec_Ac_boundary_bottom(ix,iy,1:3) = Aext0 &
                                          + ( fw%vec_Ac_m(0,ix,iy,lg%is(3),1:3) - Aext1_old )  &
                                          + coef* ( ( fw%vec_Ac_m(1,ix,iy,lg%is(3),1:3) - Aext1 ) &
                                                  - ( fw%vec_Ac_boundary_bottom_old(ix,iy,1:3) - Aext0_old ) )
        end do
        end do
      end if
      if(mg%ie(3)==lg%ie(3))then
    !$OMP parallel do collapse(2) private(ix,iy,iz)
        do iy=mg%is(2),mg%ie(2)
        do ix=mg%is(1),mg%ie(1)
          fw%vec_Ac_boundary_top(ix,iy,1:3) = fw%vec_Ac_m(0,ix,iy,lg%ie(3),1:3)   &
                                      + coef* ( fw%vec_Ac_m(1,ix,iy,lg%ie(3),1:3) - fw%vec_Ac_boundary_top_old(ix,iy,1:3) )
        end do
        end do
      end if

    !$OMP parallel do collapse(2) private(ix,iy,iz)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        fw%vec_Ac_boundary_bottom_old(ix,iy,1:3) = fw%vec_Ac_boundary_bottom(ix,iy,1:3)
        fw%vec_Ac_boundary_top_old   (ix,iy,1:3) = fw%vec_Ac_boundary_top   (ix,iy,1:3)
      end do
      end do

    end do ! ii=1,mstep

    call timer_end(LOG_SS_FDTD_CALC)

    return
  end subroutine fdtd

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

  subroutine fdtd_gbp
    implicit none

    call timer_begin(LOG_SS_FDTD_CALC)
    fw%tmp_zt = 0d0
    do iz=mg%is(3),mg%ie(3)
  !$omp parallel do collapse(2) private(iy,ix)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        fw%tmp_zt(iz,1:3) = fw%tmp_zt(iz,1:3) + fw%current4pi(ix,iy,iz,1:3) / (lg%num(1)*lg%num(2))
      end do
      end do
    end do
    call timer_end(LOG_SS_FDTD_CALC)

    call timer_begin(LOG_SS_FDTD_COMM_COLL)
    call comm_summation(fw%tmp_zt,fw%curr4pi_zt,size(fw%curr4pi_zt),info%icomm_r)
    call timer_end(LOG_SS_FDTD_COMM_COLL)

    call timer_begin(LOG_SS_FDTD_CALC)
    do ii=1,mstep

    !$OMP parallel do
      do iz=lg%is(3),lg%ie(3)

        fw%Ac_zt_m(iz,-1,1:3) = fw%Ac_zt_m(iz,0,1:3)
        fw%Ac_zt_m(iz,0 ,1:3) = fw%Ac_zt_m(iz,1,1:3)

      end do

      fw%Ac_zt_m(lg%is(3)-1,0,1:3) = fw%Ac_zt_boundary_bottom(1:3)
      fw%Ac_zt_m(lg%ie(3)+1,0,1:3) = fw%Ac_zt_boundary_top   (1:3)

!$OMP parallel do collapse(2) private(i1,iz,lap_A)
      do i1=1,3
        do iz=lg%is(3),lg%ie(3)
          lap_A = ( - 2d0* fw%Ac_zt_m(iz,0,i1) + fw%Ac_zt_m(iz-1,0,i1) + fw%Ac_zt_m(iz+1,0,i1) ) / Hgs(3)**2
          fw%Ac_zt_m(iz,1,i1) = ( cspeed_au * dt_m )**2 * lap_A &
                                    + 2.d0* fw%Ac_zt_m(iz,0,i1) - fw%Ac_zt_m(iz,-1,i1) &
                                    + dt_m**2 * fw%curr4pi_zt(iz,i1)
        end do
      end do ! i1 (spacial )

    ! external field
      tm = ( dble(itt-1) + dble(ii-1)/dble(mstep) ) *dt
      call pulse(tm,0d0,   Aext0_old)
      call pulse(tm,Hgs(3),Aext1_old)
      call pulse(tm+dt_m,0d0,   Aext0)
      call pulse(tm+dt_m,Hgs(3),Aext1)
      out_Aext = Aext1

    ! z axis: Mur absorbing boundary condition
      coef = ( cspeed_au * dt_m - Hgs(3) ) / ( cspeed_au * dt_m + Hgs(3) )

    ! absorbing boundary condition with the incident field vec_Ac_ext
      fw%Ac_zt_boundary_bottom = Aext0 &
                                      + ( fw%Ac_zt_m(lg%is(3),0,1:3) - Aext1_old )  &
                                      + coef* ( ( fw%Ac_zt_m(lg%is(3),1,1:3) - Aext1 ) &
                                              - ( fw%Ac_zt_boundary_bottom_old(1:3) - Aext0_old ) )

      fw%Ac_zt_boundary_top = fw%Ac_zt_m(lg%ie(3),0,1:3)   &
                                  + coef* ( fw%Ac_zt_m(lg%ie(3),1,1:3) - fw%Ac_zt_boundary_top_old(1:3) )

      fw%Ac_zt_boundary_bottom_old(1:3) = fw%Ac_zt_boundary_bottom(1:3)
      fw%Ac_zt_boundary_top_old   (1:3) = fw%Ac_zt_boundary_top   (1:3)

    end do ! ii=1,mstep

  !$OMP parallel do collapse(2) private(ix,iy,iz)
    do iz=mg%is(3),mg%ie(3)
    do iy=mg%is(2),mg%ie(2)
    do ix=mg%is(1),mg%ie(1)

      fw%vec_Ac_m(1,ix,iy,iz,1:3) = fw%Ac_zt_m(iz,1,1:3) + fw%Ac_fourier(ix,iy,iz,1:3)

    end do
    end do
    end do

    if(mg%is(3)==lg%is(3))then
  !$OMP parallel do collapse(2) private(ix,iy,iz)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        fw%vec_Ac_boundary_bottom(ix,iy,1:3) = fw%Ac_zt_boundary_bottom
      end do
      end do
    end if

    if(mg%ie(3)==lg%ie(3))then
!$OMP parallel do collapse(2) private(ix,iy,iz)
      do iy=mg%is(2),mg%ie(2)
      do ix=mg%is(1),mg%ie(1)
        fw%vec_Ac_boundary_top(ix,iy,1:3) = fw%Ac_zt_boundary_top
      end do
      end do
    end if
    call timer_end(LOG_SS_FDTD_CALC)

    return
  end subroutine fdtd_gbp

end subroutine fdtd_singlescale

!===================================================================================================================================

subroutine fourier_singlescale(lg,mg,info,fg,rho,j_e,Vh,poisson,singlescale)
  use structures
  use math_constants,only : zi,pi
  use salmon_global,only: dt,method_singlescale
  use phys_constants, only: cspeed_au
  use communication, only: comm_summation,comm_bcast
  implicit none
  type(s_rgrid)          ,intent(in) :: lg
  type(s_rgrid)          ,intent(in) :: mg
  type(s_parallel_info)  ,intent(in) :: info
  type(s_reciprocal_grid),intent(in) :: fg
  type(s_scalar)         ,intent(in) :: rho
  type(s_vector)         ,intent(in) :: j_e ! electron number current density (without rho*A/c)
  type(s_scalar)                     :: Vh
  type(s_poisson)                    :: poisson
  type(s_singlescale)                :: singlescale
  !
  integer :: ix,iy,iz
  integer :: iiy,iiz,iix
  real(8) :: inv_lgnum3,vec_je(3),rho_t
  integer :: i
  real(8) :: sin,cos
  complex(8) :: j0,F0,c0,s0,F_old,c_old,s_old,Ac

  if(info%isize_x < 4) stop "isize(1) must be > 3"

  i = mod(info%id_x,4) ! i=0,1,2,3

  inv_lgnum3=1.d0/(lg%num(1)*lg%num(2)*lg%num(3))

  singlescale%b_ffte=0.d0
  !$OMP parallel do collapse(2) private(iiz,iiy,ix,iy,iz,vec_je,rho_t)
  do iz=1,mg%num(3)
  do iy=1,mg%num(2)
    iiz=iz+mg%is(3)-1
    iiy=iy+mg%is(2)-1
    do ix=mg%is(1),mg%ie(1)
      singlescale%b_ffte(ix,iy,iz,0) = dcmplx(rho%f(ix,iiy,iiz)) ! charge density rho
      vec_je = ( j_e%v(1:3,ix,iiy,iiz) + singlescale%vec_je_old(1:3,ix,iiy,iiz) )*0.5d0 ! j(t) = ( j(t+dt/2) + j(t-dt/2) )/2
      rho_t  = ( rho%f(ix,iiy,iiz) + singlescale%rho_old(ix,iiy,iiz) )*0.5d0 ! rho(t) = ( rho(t+dt/2) + rho(t-dt/2) )/2
      vec_je = vec_je + rho_t * singlescale%vec_Ac_m(1,ix,iiy,iiz,1:3) ! electron number current density
      singlescale%b_ffte(ix,iy,iz,1) = dcmplx(vec_je(1)) ! current density j_x
      singlescale%b_ffte(ix,iy,iz,2) = dcmplx(vec_je(2)) ! current density j_y
      singlescale%b_ffte(ix,iy,iz,3) = dcmplx(vec_je(3)) ! current density j_z
    end do
  end do
  end do

  call comm_summation(singlescale%b_ffte,singlescale%a_ffte,size(singlescale%a_ffte),info%icomm_x)

  CALL PZFFT3DV_MOD(singlescale%a_ffte(:,:,:,i),singlescale%b_ffte(:,:,:,i),lg%num(1),lg%num(2),lg%num(3),   &
                    info%isize_y,info%isize_z,-1, &
                    info%icomm_y,info%icomm_z)

  call comm_bcast(singlescale%b_ffte(:,:,:,0),info%icomm_x, 0)

! Poisson eq.: singlescale%b_ffte(ix,iy,iz,0)=rho(G) --> poisson%b_ffte(ix,iy,iz)=Vh(G)
  poisson%zrhoG_ele=0d0
  !$omp parallel do collapse(2) default(none) &
  !$omp             private(iz,iy,ix,iiy,iiz,iix) &
  !$omp             shared(mg,lg,poisson,singlescale,inv_lgnum3,fg)
  do iz=1,mg%num(3)
  do iy=1,mg%num(2)
    iiz=iz+mg%is(3)-1
    iiy=iy+mg%is(2)-1
    do ix=1,mg%num(1)
      iix=ix+mg%is(1)-1
      poisson%zrhoG_ele(iix,iiy,iiz) = singlescale%b_ffte(iix,iy,iz,0)*inv_lgnum3
    end do
    do ix=1,lg%num(1)
      poisson%b_ffte(ix,iy,iz) = singlescale%b_ffte(ix,iy,iz,0)*fg%coef(ix,iiy,iiz)
    end do
  end do
  end do
  !$omp end parallel do

! Maxwell eq.: singlescale%b_ffte(ix,iy,iz,i)=j(G,t) --> singlescale%b_ffte(ix,iy,iz,i)=Ac(G,t+dt)
  if(i/=0) then
    !$omp parallel do collapse(2) private(iz,iy,ix,iiz,iiy,j0,F0,c0,s0,F_old,c_old,s_old,Ac,sin,cos)
    do iz=1,mg%num(3)
    do iy=1,mg%num(2)
    do ix=1,lg%num(1)
      iiz=iz+mg%is(3)-1
      iiy=iy+mg%is(2)-1
    ! j(transverse) = j - (1/(4*pi))* d(grad(phi))/dt
      j0 = singlescale%b_ffte(ix,iy,iz,i) &
      & - (1d0/(4d0*pi))* fg%coef_nabla(ix,iiy,iiz,i) * ( poisson%b_ffte(ix,iy,iz) - singlescale%Vh_ffte_old(ix,iy,iz) )/dt
    ! F_{i} = - (4*pi/(c*G)**2)* j(t)
      F0 = - (1/cspeed_au**2)*fg%coef(ix,iiy,iiz) * j0
    ! old variables
      F_old = singlescale%zf_old(ix,iy,iz,i)
      c_old = singlescale%zc_old(ix,iy,iz,i)
      s_old = singlescale%zs_old(ix,iy,iz,i)
    ! cos(c*G*dt), sin(c*G*dt)
      cos = fg%cos_cGdt(ix,iiy,iiz)
      sin = fg%sin_cGdt(ix,iiy,iiz)
    ! c_{i} = c_{i-1}* cos(c*G*dt) + s_{i-1}* sin(c*G*dt) - (F_{i}-F_{i-1})
      c0 = c_old* cos + s_old* sin - (F0-F_old)
    ! s_{i} = - c_{i-1}* sin(c*G*dt) + s_{i-1}* cos(c*G*dt)
      s0 = - c_old* sin + s_old* cos
    ! Ac(t+dt) = c_{i}* cos(c*G*dt) + s_{i}* sin(c*G*dt) + F_{i}
      Ac = c0* cos + s0* sin + F0
    ! Ac(gx=gy=0) --> 0
      singlescale%b_ffte(ix,iy,iz,i) = Ac * fg%coef_gxgy0(ix,iiy,iiz)
    ! F_old,c_old,s_old for next step
      singlescale%zf_old(ix,iy,iz,i) = F0
      singlescale%zc_old(ix,iy,iz,i) = c0
      singlescale%zs_old(ix,iy,iz,i) = s0
    end do
    end do
    end do
    !$omp end parallel do
  end if

  !$omp parallel do collapse(2) private(iz,iy,ix,iiy,iiz,iix)
  do iz=1,mg%num(3)
  do iy=1,mg%num(2)
  do ix=1,lg%num(1)
    singlescale%b_ffte(ix,iy,iz,0)    = poisson%b_ffte(ix,iy,iz)
    singlescale%Vh_ffte_old(ix,iy,iz) = poisson%b_ffte(ix,iy,iz)
  end do
  end do
  end do
  !$omp end parallel do

  CALL PZFFT3DV_MOD(singlescale%b_ffte(:,:,:,i),singlescale%a_ffte(:,:,:,i),lg%num(1),lg%num(2),lg%num(3), &
                    info%isize_y,info%isize_z,1, &
                    info%icomm_y,info%icomm_z)

  call comm_bcast(singlescale%a_ffte(:,:,:,0),info%icomm_x, 0)
  call comm_bcast(singlescale%a_ffte(:,:,:,1),info%icomm_x, 1)
  call comm_bcast(singlescale%a_ffte(:,:,:,2),info%icomm_x, 2)
  call comm_bcast(singlescale%a_ffte(:,:,:,3),info%icomm_x, 3)

  !$OMP parallel do private(iiz,iiy,ix,iy,iz)
  do iz=1,mg%num(3)
  do iy=1,mg%num(2)
    iiz=iz+mg%is(3)-1
    iiy=iy+mg%is(2)-1
    Vh%f(mg%is(1):mg%ie(1),iiy,iiz) = singlescale%a_ffte(mg%is(1):mg%ie(1),iy,iz,0)
    singlescale%Ac_fourier(mg%is(1):mg%ie(1),iiy,iiz,1:3) = singlescale%a_ffte(mg%is(1):mg%ie(1),iy,iz,1:3)
  end do
  end do

  if(method_singlescale=='1d') then
    singlescale%Ac_fourier = 0d0
  end if

  return
end subroutine fourier_singlescale

!===================================================================================================================================

subroutine init_singlescale(mg,lg,info,hgs,matrix_B,rho,Vh,srg_scalar,fw,Ac,div_Ac)
  use structures
  use sendrecv_grid, only: update_overlap_real8
  use stencil_sub, only: calc_gradient_field
  use salmon_global, only: sysname,base_directory,yn_restart,yn_ffte,method_singlescale,yn_symmetrized_stencil
  use parallelization, only: nproc_id_global
  use communication, only: comm_is_root
  use initialization_sub, only: set_bn
  use filesystem, only: open_filehandle
  use inputoutput, only: t_unit_time
  use checkpoint_restart_sub, only: restart_singlescale
  use sendrecv_grid, only: init_sendrecv_grid
  implicit none
  type(s_rgrid)         ,intent(in) :: lg,mg
  type(s_parallel_info) ,intent(in) :: info
  real(8)               ,intent(in) :: hgs(3),matrix_B(3,3)
  type(s_scalar)        ,intent(in) :: rho,Vh ! electron number density & Hartree potential
  type(s_sendrecv_grid)             :: srg_scalar
  type(s_singlescale)               :: fw
  type(s_vector)                    :: Ac
  type(s_scalar)                    :: div_Ac
  !
  character(100) :: filename
  integer :: ii,jj,ix,iy,iz
  real(8) :: bnmat(4,4)

  call allocate_scalar(mg,div_Ac)

  if(yn_symmetrized_stencil=='y') then
    call allocate_vector_with_ovlp(mg,Ac)
  else
    call allocate_vector(mg,Ac)
  end if

  fw%Energy_poynting = 0d0
  fw%Energy_joule = 0d0
  fw%curr_ave = 0d0

  call set_bn(bnmat)
  do jj=1,3
    do ii=1,4
      fw%coef_nab(ii,jj) = bnmat(ii,4)/hgs(jj)
    end do
  end do

  allocate( fw%vec_Ac_old (3,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) )
  allocate( fw%curr         (mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3),3) )
  allocate( fw%vec_je_old (3,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) )
  allocate( fw%rho_old      (mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) )
  allocate( fw%current4pi   (mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3),3) )
  allocate( fw%grad_Vh    (3,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) )
  allocate( fw%grad_Vh_old(3,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) )

!1st element: time step (-1->m-1, 0->m, 1->m+1)
!5th element: components of A vector (1->Ax, 2->Ay, 3->Az)
  allocate( fw%vec_Ac_m(-1:1,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3),1:3) )

  allocate( fw%vec_Ac_boundary_bottom(mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),1:3) &
           ,fw%vec_Ac_boundary_bottom_old(mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),1:3) &
           ,fw%vec_Ac_boundary_top(mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),1:3) &
           ,fw%vec_Ac_boundary_top_old(mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),1:3) )

  allocate(fw%integral_poynting(lg%is(3):lg%ie(3)))

  allocate(fw%rot_Ac    (3,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) &
   & ,fw%poynting_vector(3,mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) &
        & ,fw%div_Ac      (mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) &
        & ,fw%div_Ac_old  (mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3)) &
        & ,fw%integral_poynting_tmp(lg%num(3)),fw%integral_poynting_tmp2(lg%num(3)) &
!        & ,fw%Ac_zt(3,lg%is(3):lg%ie(3)) &
        & ,fw%tmp_zt(lg%is(3):lg%ie(3),3))

  fw%vec_Ac_old = 0d0
  fw%vec_Ac_m = 0d0
  fw%vec_Ac_boundary_bottom = 0d0
  fw%vec_Ac_boundary_bottom_old = 0d0
  fw%vec_Ac_boundary_top = 0d0
  fw%vec_Ac_boundary_top_old = 0d0
  fw%div_Ac_old = 0d0
  fw%integral_poynting = 0d0
  fw%curr = 0d0
  fw%vec_je_old = 0d0
!  fw%Ac_zt = 0d0
  fw%tmp_zt = 0d0

! gbp
  if(method_singlescale/='3d') then
    allocate( fw%curr4pi_zt(lg%is(3):lg%ie(3),3) )
    allocate(fw%Ac_zt_m(lg%is(3)-1:lg%ie(3)+1,-1:1,1:3))
    allocate(fw%Ac_fourier(mg%is(1):mg%ie(1),mg%is(2):mg%ie(2),mg%is(3):mg%ie(3),3))
    allocate(fw%a_ffte(lg%num(1),mg%num(2),mg%num(3),0:3),fw%b_ffte(lg%num(1),mg%num(2),mg%num(3),0:3))
    allocate(fw%Vh_ffte_old(lg%num(1),mg%num(2),mg%num(3)))
    allocate(fw%zf_old(lg%num(1),mg%num(2),mg%num(3),0:3), &
           & fw%zc_old(lg%num(1),mg%num(2),mg%num(3),0:3), &
           & fw%zs_old(lg%num(1),mg%num(2),mg%num(3),0:3))
    fw%curr4pi_zt = 0d0
    fw%Ac_zt_m = 0d0
    fw%Ac_zt_boundary_bottom = 0d0
    fw%Ac_zt_boundary_top = 0d0
    fw%Ac_zt_boundary_bottom_old = 0d0
    fw%Ac_zt_boundary_top_old = 0d0
    fw%Ac_fourier = 0d0
    fw%zf_old = 0d0
    fw%zc_old = 0d0
    fw%zs_old = 0d0
    if(yn_ffte=='y') then
    ! Vh --> fw%Vh_ffte_old
      call calc_Vh_ffte
    end if
  end if

  if(comm_is_root(nproc_id_global)) then
    write(filename,"(2A,'_rt_micro.data')") trim(base_directory),trim(SYSname)
    fw%fh_rt_micro = open_filehandle(filename)
    write(fw%fh_rt_micro, '("#",99(1X,I0,":",A,"[",A,"]"))') &
      & 1, "time", trim(t_unit_time%name), &
      & 2, "Ac_tot_x(z=0)", "a.u.", &
      & 3, "Ac_tot_y(z=0)", "a.u.", &
      & 4, "Ac_tot_z(z=0)", "a.u.", &
      & 5, "Ac_tot_x(z=L)", "a.u.", &
      & 6, "Ac_tot_y(z=L)", "a.u.", &
      & 7, "Ac_tot_z(z=L)", "a.u.", &
      & 8, "Ac_ext_x",      "a.u.", &
      & 9, "Ac_ext_y",      "a.u.", &
     & 10, "Ac_ext_z",      "a.u.", &
     & 11, "J_x(w/o Ac)",   "a.u.", &
     & 12, "J_y(w/o Ac)",   "a.u.", &
     & 13, "J_z(w/o Ac)",   "a.u.", &
     & 14, "E_electron",    "a.u.", &
     & 15, "E_poynting(z=0)", "a.u.", &
     & 16, "E_poynting(z=L)", "a.u.", &
     & 17, "E_em",            "a.u.", &
     & 18, "E_joule",         "a.u."

!  ! for spatial distribution of excitation energy
!    write(filename,"(2A,'_excitation.data')") trim(base_directory),trim(SYSname)
!    fw%fh_excitation = open_filehandle(filename)

!  ! for the vector potential Ac(z,t)
!    write(filename,"(2A,'_Ac_zt.data')") trim(base_directory),trim(SYSname)
!    fw%fh_Ac_zt = open_filehandle(filename)
  end if

  allocate(fw%box1(mg%is_array(1):mg%ie_array(1), &
  & mg%is_array(2):mg%ie_array(2), &
  & mg%is_array(3):mg%ie_array(3)))
  fw%box1 = 0d0

!$OMP parallel do collapse(2) private(ix,iy,iz)
  do iz=mg%is(3),mg%ie(3)
  do iy=mg%is(2),mg%ie(2)
  do ix=mg%is(1),mg%ie(1)
    fw%box1(ix,iy,iz) = Vh%f(ix,iy,iz)
    fw%rho_old(ix,iy,iz) = rho%f(ix,iy,iz)
  end do
  end do
  end do
  if(info%if_divide_rspace) call update_overlap_real8(srg_scalar, mg, fw%box1)
  call calc_gradient_field(mg,fw%coef_nab,matrix_B,fw%box1,fw%grad_Vh_old)

  if(yn_restart=='y') then
    call restart_singlescale(info%icomm_rko,lg,mg,fw,Ac,div_Ac)
  end if

! specialized in FDTD timestep
  fw%eg%nd = 1
  fw%eg%is = mg%is
  fw%eg%ie = mg%ie
  fw%eg%is_array = mg%is - 1
  fw%eg%ie_array = mg%ie + 1
  call init_sendrecv_grid(fw%srg_eg, fw%eg, 1, srg_scalar%icomm, srg_scalar%neig)
  allocate(fw%box(fw%eg%is_array(1):fw%eg%ie_array(1), &
                & fw%eg%is_array(2):fw%eg%ie_array(2), &
                & fw%eg%is_array(3):fw%eg%ie_array(3)))

  return

contains

  subroutine calc_Vh_ffte
    use communication, only: comm_summation
    implicit none
    integer :: ix,iy,iz
    integer :: iiy,iiz,iix

    fw%b_ffte = 0d0
  !$OMP parallel do private(iiz,iiy,ix)
    do iz=1,mg%num(3)
    do iy=1,mg%num(2)
      iiz=iz+mg%is(3)-1
      iiy=iy+mg%is(2)-1
      fw%b_ffte(mg%is(1):mg%ie(1),iy,iz,0) = Vh%f(mg%is(1):mg%ie(1),iiy,iiz)
    end do
    end do
    call comm_summation(fw%b_ffte,fw%a_ffte,size(fw%a_ffte),info%icomm_x)

    CALL PZFFT3DV_MOD(fw%a_ffte(:,:,:,0),fw%Vh_ffte_old,lg%num(1),lg%num(2),lg%num(3),   &
                      info%isize_y,info%isize_z,-1, &
                      info%icomm_y,info%icomm_z)

    return
  end subroutine calc_Vh_ffte

end subroutine init_singlescale

!===================================================================================================================================

! line integral \int_{{\bf r}_0 \rightarrow {\bf r}_1} {\bf A}({\bf r}) \cdot d{\bf x}
! path: r0 --> r1 = (ix1*Hx,iy1*Hy,iz1*Hz)
subroutine line_integral(integral,r0,A,nx,ny,nz,ix1,iy1,iz1,Hx,Hy,Hz &
                        ,A_lerp,line,wrk,n_max,index)
  implicit none
  integer,intent(in)  :: nx,ny,nz,ix1,iy1,iz1,n_max
  real(8),intent(in)  :: r0(3),A(3,0:nx-1,0:ny-1,0:nz-1),Hx,Hy,Hz
  real(8),intent(out) :: integral
  integer             :: index(n_max)
  real(8)             :: A_lerp(3,n_max),line(3,n_max),wrk(n_max)
  !
  integer :: ix,iy,iz,ix0,iy0,iz0,i,j,k,n,ixp,iyp,izp
  real(8) :: r1(3),p(3),q(3),r(3),h(3),l(3),t

# define IMOD(ix,iy,iz) modulo(ix,nx),modulo(iy,ny),modulo(iz,nz)

  r1 = (/ ix1*Hx, iy1*Hy, iz1*Hz /)

  ix0 = floor( r0(1)/Hx )
  iy0 = floor( r0(2)/Hy )
  iz0 = floor( r0(3)/Hz )

  h = (/ Hx, Hy, Hz /)

! trilinear interpolation @ r0
  ix = ix0
  iy = iy0
  iz = iz0
  ixp = ix + 1
  iyp = iy + 1
  izp = iz + 1
  p = (/ ix*Hx, iy*Hy, iz*Hz /)
  r = ( r0 - p )/h
  A_lerp(:,1) = A(:,IMOD(ix,iy,iz)) * ( 1d0 - r(1) ) * ( 1d0 - r(2) ) * ( 1d0 - r(3) ) &
              + A(:,IMOD(ixp,iy,iz)) * r(1) * ( 1d0 - r(2) ) * ( 1d0 - r(3) ) &
              + A(:,IMOD(ix,iyp,iz)) * ( 1d0 - r(1) ) * r(2) * ( 1d0 - r(3) ) &
              + A(:,IMOD(ix,iy,izp)) * ( 1d0 - r(1) ) * ( 1d0 - r(2) ) * r(3) &
              + A(:,IMOD(ix,iyp,izp)) * ( 1d0 - r(1) ) * r(2) * r(3) &
              + A(:,IMOD(ixp,iy,izp)) * r(1) * ( 1d0 - r(2) ) * r(3) &
              + A(:,IMOD(ixp,iyp,iz)) * r(1) * r(2) * ( 1d0 - r(3) ) &
              + A(:,IMOD(ixp,iyp,izp)) * r(1) * r(2) * r(3)
  line(:,1) = 0d0
  wrk(1) = 0d0

! bilinear interpolation @ intersections
  i = 1
  do ix=min(ix0,ix1),max(ix0,ix1)
    if(ix0==ix1) cycle
    t = ( ix*Hx - r0(1) ) / ( r1(1) - r0(1) )
    if(t<0d0 .or. 1d0<t) cycle
    l = r0 + t * ( r1 - r0 )
    iy = floor( l(2)/Hy )
    iz = floor( l(3)/Hz )
    iyp = iy + 1
    izp = iz + 1
    q = (/ ix*Hx, iy*Hy, iz*Hz /)
    r = ( l - q )/h
    i = i + 1
    if(i>n_max) stop "n_max is too small"
    A_lerp(:,i) = A(:,IMOD(ix,iy,iz)) * ( 1d0 - r(2) ) * ( 1d0 - r(3) ) &
                + A(:,IMOD(ix,iyp,iz)) * r(2) * ( 1d0 - r(3) ) &
                + A(:,IMOD(ix,iy,izp)) * ( 1d0 - r(2) ) * r(3) &
                + A(:,IMOD(ix,iyp,izp)) * r(2) * r(3)
    line(:,i) = l - r0
    wrk(i) = line(1,i)**2 + line(2,i)**2 + line(3,i)**2
  end do
  do iy=min(iy0,iy1),max(iy0,iy1)
    if(iy0==iy1) cycle
    t = ( iy*Hy - r0(2) ) / ( r1(2) - r0(2) )
    if(t<0d0 .or. 1d0<t) cycle
    l = r0 + t * ( r1 - r0 )
    ix = floor( l(1)/Hx )
    iz = floor( l(3)/Hz )
    ixp = ix + 1
    izp = iz + 1
    q = (/ ix*Hx, iy*Hy, iz*Hz /)
    r = ( l - q )/h
    i = i + 1
    if(i>n_max) stop "n_max is too small"
    A_lerp(:,i) = A(:,IMOD(ix,iy,iz)) * ( 1d0 - r(1) ) * ( 1d0 - r(3) ) &
                + A(:,IMOD(ixp,iy,iz)) * r(1) * ( 1d0 - r(3) ) &
                + A(:,IMOD(ix,iy,izp)) * ( 1d0 - r(1) ) * r(3) &
                + A(:,IMOD(ixp,iy,izp)) * r(1) * r(3)
    line(:,i) = l - r0
    wrk(i) = line(1,i)**2 + line(2,i)**2 + line(3,i)**2
  end do
  do iz=min(iz0,iz1),max(iz0,iz1)
    if(iz0==iz1) cycle
    t = ( iz*Hz - r0(3) ) / ( r1(3) - r0(3) )
    if(t<0d0 .or. 1d0<t) cycle
    l = r0 + t * ( r1 - r0 )
    ix = floor( l(1)/Hx )
    iy = floor( l(2)/Hy )
    ixp = ix + 1
    iyp = iy + 1
    q = (/ ix*Hx, iy*Hy, iz*Hz /)
    r = ( l - q )/h
    i = i + 1
    if(i>n_max) stop "n_max is too small"
    A_lerp(:,i) = A(:,IMOD(ix,iy,iz)) * ( 1d0 - r(1) ) * ( 1d0 - r(2) ) &
                + A(:,IMOD(ixp,iy,iz)) * r(1) * ( 1d0 - r(2) ) &
                + A(:,IMOD(ix,iyp,iz)) * ( 1d0 - r(1) ) * r(2) &
                + A(:,IMOD(ixp,iyp,iz)) * r(1) * r(2)
    line(:,i) = l - r0
    wrk(i) = line(1,i)**2 + line(2,i)**2 + line(3,i)**2
  end do
  if(i==1) then
    i = i + 1
    if(i>n_max) stop "n_max is too small"
    A_lerp(:,i) = A(:,IMOD(ix0,iy0,iz0))
    line(:,i) = p - r0
    wrk(i) = line(1,i)**2 + line(2,i)**2 + line(3,i)**2
  end if
  n = i

  call heapsort(n,wrk,index)

  integral = 0d0
  do i=2,n
    j = index(i)
    k = index(i-1)
    integral = integral + 0.5d0* ( A_lerp(1,j) + A_lerp(1,k) ) &
                        * ( line(1,j) - line(1,k) ) &
                        + 0.5d0* ( A_lerp(2,j) + A_lerp(2,k) ) &
                        * ( line(2,j) - line(2,k) ) &
                        + 0.5d0* ( A_lerp(3,j) + A_lerp(3,k) ) &
                        * ( line(3,j) - line(3,k) )
  end do

  return
contains

  subroutine heapsort(n,array,index)
    implicit none
    integer,intent(in)    :: n
    real(8),intent(inout) :: array(n)
    integer,intent(out)   :: index(n)
    !
    integer :: i,k,j,l,m
    real(8) :: t

    if(n < 2) return

    do i=1,n
      index(i) = i
    end do

    l = n/2 + 1
    k = n
    do while(k /= 1)
      if(l > 1)then
        l = l - 1
        t = array(l)
        m = index(l)
      else
        t = array(k)
        m = index(k)
        array(k) = array(1)
        index(k) = index(1)
        k = k - 1
        if(k == 1) then
          array(1) = t
          index(1) = m
          exit
        end if
      end if
      i = l
      j = l + l
      do while(j <= k)
        if(j < k) then
          if(array(j) < array(j+1)) j = j + 1
        endif
        if (t < array(j))then
          array(i) = array(j)
          index(i) = index(j)
          i = j
          j = j + j
        else
          j = k + 1
        end if
      end do
      array(i) = t
      index(i) = m
    end do

    return
  end subroutine heapsort
end subroutine line_integral

end module fdtd_coulomb_gauge


text:
      MODULE OBSCHDMODULE
         INTEGER, SAVE, POINTER    ::NQCH,NQCCH,NQTCH,IUCHOBSV,IPRT
         INTEGER, SAVE, DIMENSION(:),   POINTER ::NQOBCH
         INTEGER, SAVE, DIMENSION(:),   POINTER ::NQCLCH
         INTEGER, SAVE, DIMENSION(:),   POINTER ::IOBTS
         REAL,    SAVE, DIMENSION(:),   POINTER ::FLWSIM
         REAL,    SAVE, DIMENSION(:),   POINTER ::FLWOBS
         REAL,    SAVE, DIMENSION(:),   POINTER ::TOFF
         REAL,    SAVE, DIMENSION(:),   POINTER ::OTIME
         REAL,    SAVE, DIMENSION(:,:), POINTER ::QCELL
         CHARACTER*12,SAVE,DIMENSION(:),POINTER ::OBSNAM
      TYPE OBSCHDTYPE
         INTEGER,            POINTER ::NQCH,NQCCH,NQTCH,IUCHOBSV,IPRT
         INTEGER,     DIMENSION(:),  POINTER ::NQOBCH
         INTEGER,     DIMENSION(:),  POINTER ::NQCLCH
         INTEGER,     DIMENSION(:),  POINTER ::IOBTS
         REAL,        DIMENSION(:),  POINTER ::FLWSIM
         REAL,        DIMENSION(:),  POINTER ::FLWOBS
         REAL,        DIMENSION(:),  POINTER ::TOFF
         REAL,        DIMENSION(:),  POINTER ::OTIME
         REAL,        DIMENSION(:,:),POINTER ::QCELL
         CHARACTER*12,DIMENSION(:),  POINTER ::OBSNAM
      END TYPE
      TYPE(OBSCHDTYPE), SAVE  ::OBSCHDDAT(10)
      END MODULE OBSCHDMODULE

      SUBROUTINE OBS2CHD7AR(IUCHOB,IGRID)
C     ******************************************************************
C     ALLOCATE AND READ DATA FOR FLOW OBSERVATIONS AT CONSTANT-HEAD
C     BOUNDARY CELLS
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL, ONLY: NCOL,NROW,NLAY,NPER,NSTP,PERLEN,TSMULT,ISSFLG,
     1                  IOUT,ITRSS
      USE OBSCHDMODULE
C
      CHARACTER*200 LINE
C     ------------------------------------------------------------------
      ALLOCATE(NQCH,NQTCH,NQCCH,IUCHOBSV,IPRT)
C
C1------INITIALIZE VARIABLEA.
      ZERO=0.0
      IERR=0
      NT=0
      NC=0
C
C2------IDENTIFY PROCESS
      WRITE(IOUT,14) IUCHOB
   14 FORMAT(/,' OBS2CHD7 -- CONSTANT-HEAD BOUNDARY FLOW OBSERVATIONS',
     &    /,' VERSION 2.0, 02/28/2006       INPUT READ FROM UNIT ',I3)
C
C3------ITEM 1
      CALL URDCOM(IUCHOB,IOUT,LINE)
      LLOC = 1
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NQCH,DUM,IOUT,IUCHOB)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NQCCH,DUM,IOUT,IUCHOB)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,NQTCH,DUM,IOUT,IUCHOB)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,2,IUCHOBSV,DUM,IOUT,IUCHOB)
      CALL URWORD(LINE,LLOC,ISTART,ISTOP,1,IDUM,DUM,IOUT,IUCHOB)
      IPRT=1
      IF(LINE(ISTART:ISTOP).EQ.'NOPRINT') THEN
        IPRT=0
        WRITE(IOUT,*) 'NOPRINT option for CONSTANT-HEAD OBSERVATIONS'
      END IF
      WRITE (IOUT,17) NQCH, NQCCH, NQTCH
   17 FORMAT (/,
     &    ' NUMBER OF FLOW-OBSERVATION CONSTANT-HEAD-CELL GROUPS:',I5,/,
     &    '   NUMBER OF CELLS IN CONSTANT-HEAD-CELL GROUPS......:',I5,/,
     &    '   NUMBER OF CONSTANT-HEAD-CELL FLOWS................:',I5)
      IF(NQTCH.LE.0) THEN
         WRITE(IOUT,*) ' NQTCH LESS THAN OR EQUAL TO 0'
         CALL USTOP(' ')
      END IF
      IF(IUCHOBSV.GT.0) THEN
         WRITE(IOUT,21) IUCHOBSV
   21    FORMAT(1X,
     1      'CH OBSERVATIONS WILL BE SAVED ON UNIT...............:',I5)
      ELSE
         WRITE(IOUT,22)
   22    FORMAT(1X,'CH OBSERVATIONS WILL NOT BE SAVED IN A FILE')
      END IF
C
C4------ALLOCATE ARRAYS
      ALLOCATE (NQOBCH(NQCH))
      ALLOCATE (NQCLCH(NQCH))
      ALLOCATE (IOBTS(NQTCH))
      ALLOCATE (FLWSIM(NQTCH))
      ALLOCATE (FLWOBS(NQTCH))
      ALLOCATE (TOFF(NQTCH))
      ALLOCATE (OTIME(NQTCH))
      ALLOCATE (QCELL(4,NQCCH))
      ALLOCATE (OBSNAM(NQTCH))
      DO 19 N=1,NQTCH
      OTIME(N)=ZERO
      FLWSIM(N)=ZERO
   19 CONTINUE
C
C5------READ AND WRITE TIME-OFFSET MULTIPLIER FOR FLOW-OBSERVATION TIMES
      READ(IUCHOB,*) TOMULTCH
      IF(IPRT.NE.0) WRITE (IOUT,520) TOMULTCH
  520 FORMAT (/,' OBSERVED CONSTANT-HEAD-CELL FLOW DATA',/,
     &' -- TIME OFFSETS ARE MULTIPLIED BY: ',G12.5)
C
C6------LOOP THROUGH CELL GROUPS.
      DO 120 IQ = 1,NQCH
C
C7------READ ITEM 3
        READ (IUCHOB,*) NQOBCH(IQ), NQCLCH(IQ)
        IF(IPRT.NE.0) WRITE (IOUT,525) IQ, 'CHD', NQCLCH(IQ), NQOBCH(IQ)
  525   FORMAT (/,'   GROUP NUMBER: ',I3,'   BOUNDARY TYPE: ',A,
     &         '   NUMBER OF CELLS IN GROUP: ',I5,/,
     &         '   NUMBER OF FLOW OBSERVATIONS: ',I5,//,
     &         40X,'OBSERVED',/,
     &         20X,'REFER.',12X,'BOUNDARY FLOW',/,
     &      7X,'OBSERVATION',2X,'STRESS',4X,'TIME',5X,'GAIN (-) OR',/,
     &         2X,'OBS#    NAME',6X,'PERIOD   OFFSET',5X,'LOSS (+)')
C
C8------SET FLAG FOR SETTING ALL FACTORS TO 1
        IFCTFLG = 0
        IF (NQCLCH(IQ).LT.0) THEN
          IFCTFLG = 1
          NQCLCH(IQ) = -NQCLCH(IQ)
        ENDIF
C
C9------READ TIME STEPS, MEASURED FLOWS, AND WEIGHTS.
        NT1 = 1 + NT
        NT2 = NQOBCH(IQ) + NT
        DO 30 N = NT1, NT2
C
C10-----READ ITEM 4
          READ (IUCHOB,*) OBSNAM(N), IREFSP, TOFFSET, FLWOBS(N)
          IF(IPRT.NE.0) WRITE (IOUT,535) N, OBSNAM(N), IREFSP, TOFFSET,
     1                                   FLWOBS(N)
  535     FORMAT(1X,I5,1X,A12,2X,I4,2X,G11.4,1X,G11.4)
          CALL UOBSTI(OBSNAM(N),IOUT,ISSFLG,ITRSS,NPER,NSTP,IREFSP,
     &                IOBTS(N),PERLEN,TOFF(N),TOFFSET,TOMULTCH,TSMULT,1,
     &                OTIME(N))
   30   CONTINUE
C
C11-----READ LAYER, ROW, COLUMN, AND FACTOR (ITEM 5)
        NC1 = NC + 1
        NC2 = NC + NQCLCH(IQ)
        IF(IPRT.NE.0) WRITE (IOUT,540)
  540   FORMAT (/,'       LAYER  ROW  COLUMN    FACTOR')
        DO 40 L = NC1, NC2
          READ (IUCHOB,*) (QCELL(I,L),I=1,4)
          IF(QCELL(4,L).EQ.0. .OR. IFCTFLG.EQ.1) QCELL(4,L) = 1.
          IF(IPRT.NE.0) WRITE (IOUT,550) (QCELL(I,L),I=1,4)
  550     FORMAT (4X,F8.0,F6.0,F7.0,F9.2)
          K = QCELL(1,L)
          I = QCELL(2,L)
          J = QCELL(3,L)
          IF (K.LE.0 .OR. K.GT.NLAY .OR .J.LE.0 .OR. J.GT.NCOL .OR.
     &        I.LE.0 .OR. I.GT.NROW) THEN
            WRITE (IOUT,590)
  590       FORMAT (/,' ROW OR COLUMN NUMBER INVALID',
     &        ' -- STOP EXECUTION (OBS2CHD7AR)',/)
            IERR = 1
          ENDIF
   40   CONTINUE
C
C12-----UPDATE COUNTERS.
        NC = NC2
        NT = NT2
  120 CONTINUE
C
C13-----STOP IF THERE WERE ANY ERRORS WHILE READING.
      IF (IERR.GT.0) THEN
        WRITE(IOUT,620)
  620   FORMAT (/,' ERROR:  SEE ABOVE FOR ERROR MESSAGE AND "STOP',
     &        ' EXECUTION" (OBS2CHD7AR)')
        CALL USTOP(' ')
      ENDIF
C
C14-----RETURN.
      CALL SOBS2CHD7PSV(IGRID)
      RETURN
      END
      SUBROUTINE OBS2CHD7SE(KKPER,IGRID)
C     ******************************************************************
C     CALCULATE SIMULATED EQUIVALENTS TO OBSERVED CONSTANT-HEAD FLOWS
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,    ONLY:IBOUND,IOUT
      USE OBSBASMODULE,ONLY:ITS
      USE OBSCHDMODULE
C
      DOUBLE PRECISION RATE
C     ------------------------------------------------------------------
      CALL SOBS2CHD7PNT(IGRID)
C
C1------INITIALIZE VARIABLES
      ZERO = 0.0
      NC = 0
      NT1 = 1
C
C2------LOOP THROUGH BOUNDARY FLOW CELL GROUPS
      DO 60 IQ = 1, NQCH
        NT2 = NT1 + NQOBCH(IQ) - 1
C
C3--------LOOP THROUGH THE OBSERVATION TIMES FOR THIS CELL GROUP.
        DO 40 NT = NT1, NT2
C
C4--------WAS THERE A MEASUREMENT AT THIS BOUNDARY THIS TIME STEP?
          IF (IOBTS(NT).EQ.ITS .OR.
     &        (IOBTS(NT).EQ.ITS-1.AND.TOFF(NT).GT.ZERO)) THEN
C
C5------YES -- LOOP THROUGH CELLS.
            NC1 = NC + 1
            NC2 = NC + NQCLCH(IQ)
            DO 30 N = NC1, NC2
              K = QCELL(1,N)
              I = QCELL(2,N)
              J = QCELL(3,N)
              IF (IBOUND(J,I,K).GE.0) THEN
                WRITE(IOUT,500) K,I,J,KKPER
  500           FORMAT(/,
     &' *** ERROR: CONSTANT-HEAD FLOW OBSERVATION SPECIFIED FOR CELL (',
     &I3,',',I5,',',I5,'),',/,
     &12X,'BUT THIS CELL IS NOT CONSTANT-HEAD IN STRESS PERIOD ',I4,/
     &12X,'-- STOP EXECUTION (OBS2CHD7SE)')
                CALL USTOP(' ')
              ENDIF
C
C6------CALL SUBROUTINE TO CALCULATE CONSTANT-HEAD FLOW FOR CELL
              CALL SOBS2CHD7FFLW(J,I,K,RATE)
C
C7------SUM VALUES FROM INDIVIDUAL CELLS.
C7------CALCULATE FACTOR FOR TEMPORAL INTERPOLATION
   20         FACT = 1.0
              IF (TOFF(NT).GT.ZERO) THEN
                IF (IOBTS(NT).EQ.ITS) FACT = 1. - TOFF(NT)
                IF (IOBTS(NT).EQ.ITS-1) FACT = TOFF(NT)
              ENDIF
C
C8------ACCUMULATE FLOWS FOR THE SIMULATED OBSERVATION.
              FLWSIM(NT) = FLWSIM(NT) + RATE*FACT*QCELL(4,N)
C
Cx------END OF LOOP FOR CELLS IN ONE GROUP.
   30       CONTINUE
C
          ENDIF
C
Cx------END OF LOOP FOR OBSERVATION TIMES FOR ONE GROUP
   40   CONTINUE
C
C10-----UPDATE CELL AND TIME COUNTERS
        NC = NC + NQCLCH(IQ)
        NT1 = NT2 + 1
C
Cx------END OF LOOP FOR CELL GROUPS.
   60 CONTINUE
C
C11-----RETURN.
      RETURN
      END
      SUBROUTINE SOBS2CHD7FFLW(J,I,K,RATE)
C     ******************************************************************
C     CALCULATE CONSTANT-HEAD BOUNDARY FLOW FOR A GIVEN CELL
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL,       ONLY:IBOUND,HNEW,CR,CC,CV,BOTM,NBOTM,
     1                       NCOL,NROW,NLAY,LAYHDT,LBOTM
      USE GWFBASMODULE, ONLY:ICHFLG
C
      DOUBLE PRECISION HD,X1,X2,X3,X4,X5,X6,RATE
C     ------------------------------------------------------------------
C
C6------CLEAR VALUES FOR FLOW RATE THROUGH EACH FACE OF CELL.
      ZERO=0.
      X1=ZERO
      X2=ZERO
      X3=ZERO
      X4=ZERO
      X5=ZERO
      X6=ZERO
C
C7------CALCULATE FLOW THROUGH THE LEFT FACE.
C7------COMMENTS A-C APPEAR ONLY IN THE SECTION HEADED BY COMMENT 7,
C7------BUT THEY APPLY IN A SIMILAR MANNER TO SECTIONS 8-12.
C
C7A-----IF THERE IS NO FLOW TO CALCULATE THROUGH THIS FACE, THEN GO ON
C7A-----TO NEXT FACE.  NO FLOW OCCURS AT THE EDGE OF THE GRID, TO AN
C7A-----ADJACENT NO-FLOW CELL, OR TO AN ADJACENT CONSTANT-HEAD CELL
C7A-----WHEN ICHFLG IS 0.
      IF(J.EQ.1) GO TO 30
      IF(IBOUND(J-1,I,K).EQ.0) GO TO 30
      IF(ICHFLG.EQ.0 .AND. IBOUND(J-1,I,K).LT.0) GO TO 30
C
C7B-----CALCULATE FLOW THROUGH THIS FACE INTO THE ADJACENT CELL.
      HDIFF=HNEW(J,I,K)-HNEW(J-1,I,K)
      X1=HDIFF*CR(J-1,I,K)
C
C8------CALCULATE FLOW THROUGH THE RIGHT FACE.
   30 IF(J.EQ.NCOL) GO TO 60
      IF(IBOUND(J+1,I,K).EQ.0) GO TO 60
      IF(ICHFLG.EQ.0 .AND. IBOUND(J+1,I,K).LT.0) GO TO 60
      HDIFF=HNEW(J,I,K)-HNEW(J+1,I,K)
      X2=HDIFF*CR(J,I,K)
C
C9------CALCULATE FLOW THROUGH THE BACK FACE.
   60 IF(I.EQ.1) GO TO 90
      IF (IBOUND(J,I-1,K).EQ.0) GO TO 90
      IF(ICHFLG.EQ.0 .AND. IBOUND(J,I-1,K).LT.0) GO TO 90
      HDIFF=HNEW(J,I,K)-HNEW(J,I-1,K)
      X3=HDIFF*CC(J,I-1,K)
C
C10-----CALCULATE FLOW THROUGH THE FRONT FACE.
   90 IF(I.EQ.NROW) GO TO 120
      IF(IBOUND(J,I+1,K).EQ.0) GO TO 120
      IF(ICHFLG.EQ.0 .AND. IBOUND(J,I+1,K).LT.0) GO TO 120
      HDIFF=HNEW(J,I,K)-HNEW(J,I+1,K)
      X4=HDIFF*CC(J,I,K)
C
C11-----CALCULATE FLOW THROUGH THE UPPER FACE.
  120 IF(K.EQ.1) GO TO 150
      IF (IBOUND(J,I,K-1).EQ.0) GO TO 150
      IF(ICHFLG.EQ.0 .AND. IBOUND(J,I,K-1).LT.0) GO TO 150
      HD=HNEW(J,I,K)
      IF(LAYHDT(K).EQ.0) GO TO 122
      TMP=HD
      TOP=BOTM(J,I,LBOTM(K)-1)
      IF(TMP.LT.TOP) HD=TOP
  122 HDIFF=HD-HNEW(J,I,K-1)
      X5=HDIFF*CV(J,I,K-1)
C
C12-----CALCULATE FLOW THROUGH THE LOWER FACE.
  150 IF(K.EQ.NLAY) GO TO 180
      IF(IBOUND(J,I,K+1).EQ.0) GO TO 180
      IF(ICHFLG.EQ.0 .AND. IBOUND(J,I,K+1).LT.0) GO TO 180
      HD=HNEW(J,I,K+1)
      IF(LAYHDT(K+1).EQ.0) GO TO 152
      TMP=HD
      TOP=BOTM(J,I,LBOTM(K+1)-1)
      IF(TMP.LT.TOP) HD=TOP
  152 HDIFF=HNEW(J,I,K)-HD
      X6=HDIFF*CV(J,I,K)
C
C13-----SUM THE FLOWS THROUGH SIX FACES OF CONSTANT HEAD CELL
 180  RATE=X1+X2+X3+X4+X5+X6
C
C-----RETURN
      RETURN
      END
      SUBROUTINE OBS2CHD7OT(IGRID)
C     ******************************************************************
C     WRITE ALL OBSERVATIONS TO LISTING FILE.
C     ******************************************************************
C        SPECIFICATIONS:
C     ------------------------------------------------------------------
      USE GLOBAL, ONLY: IOUT
      USE OBSCHDMODULE
      DOUBLE PRECISION SQ,SUMSQ
C     ------------------------------------------------------------------
      CALL SOBS2CHD7PNT(IGRID)
C
C1------WRITE OBSERVATIONS TO LISTING FILE.
      IF(IPRT.NE.0) WRITE(IOUT,17)
   17 FORMAT(1X,/,1X,'CONSTANT HEAD FLOW OBSERVATIONS',/,
     1  1X,'OBSERVATION       OBSERVED           SIMULATED',/
     2  1X,'  NAME              VALUE              VALUE',
     3     '             DIFFERENCE',/
     4  1X,'----------------------------------------------',
     5     '----------------------')
      SUMSQ=0.
      DO 100 N=1,NQTCH
      DIFF=FLWOBS(N)-FLWSIM(N)
      SQ=DIFF*DIFF
      SUMSQ=SUMSQ+SQ
      IF(IPRT.NE.0) WRITE(IOUT,27) OBSNAM(N),FLWOBS(N),FLWSIM(N),DIFF
   27 FORMAT(1X,A,1P,3G20.11)
  100 CONTINUE
      WRITE(IOUT,28) SUMSQ
   28 FORMAT(1X,/,1X,'CHD FLOW SUM OF SQUARED DIFFERENCE:',1P,E15.5)
C
C2------WRITE OBSERVATIONS TO SEPARATE FILE.
      IF(IUCHOBSV.GT.0) CALL UOBSSV(IUCHOBSV,NQTCH,FLWSIM,FLWOBS,
     1                              OBSNAM,0)
C
C3------RETURN.
      RETURN
      END
      SUBROUTINE OBS2CHD7DA(IGRID)
C  Deallocate OBSCHD memory
      USE OBSCHDMODULE
C
      DEALLOCATE(NQCH)
      DEALLOCATE(NQTCH)
      DEALLOCATE(NQCCH)
      DEALLOCATE(IUCHOBSV)
      DEALLOCATE(IPRT)
      DEALLOCATE(NQOBCH)
      DEALLOCATE(NQCLCH)
      DEALLOCATE(IOBTS)
      DEALLOCATE(FLWSIM)
      DEALLOCATE(FLWOBS)
      DEALLOCATE(TOFF)
      DEALLOCATE(OTIME)
      DEALLOCATE(QCELL)
      DEALLOCATE(OBSNAM)
C
      RETURN
      END
      SUBROUTINE SOBS2CHD7PNT(IGRID)
C  Change OBSCHD data to a different grid.
      USE OBSCHDMODULE
C
      NQCH=>OBSCHDDAT(IGRID)%NQCH
      NQTCH=>OBSCHDDAT(IGRID)%NQTCH
      NQCCH=>OBSCHDDAT(IGRID)%NQCCH
      IUCHOBSV=>OBSCHDDAT(IGRID)%IUCHOBSV
      IPRT=>OBSCHDDAT(IGRID)%IPRT
      NQOBCH=>OBSCHDDAT(IGRID)%NQOBCH
      NQCLCH=>OBSCHDDAT(IGRID)%NQCLCH
      IOBTS=>OBSCHDDAT(IGRID)%IOBTS
      FLWSIM=>OBSCHDDAT(IGRID)%FLWSIM
      FLWOBS=>OBSCHDDAT(IGRID)%FLWOBS
      TOFF=>OBSCHDDAT(IGRID)%TOFF
      OTIME=>OBSCHDDAT(IGRID)%OTIME
      QCELL=>OBSCHDDAT(IGRID)%QCELL
      OBSNAM=>OBSCHDDAT(IGRID)%OBSNAM
C
      RETURN
      END
      SUBROUTINE SOBS2CHD7PSV(IGRID)
C  Save OBSCHD data for a grid.
      USE OBSCHDMODULE
C
      OBSCHDDAT(IGRID)%NQCH=>NQCH
      OBSCHDDAT(IGRID)%NQTCH=>NQTCH
      OBSCHDDAT(IGRID)%NQCCH=>NQCCH
      OBSCHDDAT(IGRID)%IUCHOBSV=>IUCHOBSV
      OBSCHDDAT(IGRID)%IPRT=>IPRT
      OBSCHDDAT(IGRID)%NQOBCH=>NQOBCH
      OBSCHDDAT(IGRID)%NQCLCH=>NQCLCH
      OBSCHDDAT(IGRID)%IOBTS=>IOBTS
      OBSCHDDAT(IGRID)%FLWSIM=>FLWSIM
      OBSCHDDAT(IGRID)%FLWOBS=>FLWOBS
      OBSCHDDAT(IGRID)%TOFF=>TOFF
      OBSCHDDAT(IGRID)%OTIME=>OTIME
      OBSCHDDAT(IGRID)%QCELL=>QCELL
      OBSCHDDAT(IGRID)%OBSNAM=>OBSNAM
C
      RETURN
      END

ext: f

max_stars_repo_name: cneyens/MODFLOW-MAXDRN

max_issues_repo_name: cneyens/MODFLOW-MAXDRN

max_forks_repo_name: cneyens/MODFLOW-MAXDRN


text:
      subroutine cal_displace
      include 'moldy.h'
c
      integer*4 i, j, ii, ihe, ijk
      real*8 xij, yij, zij, rxij, ryij, rzij, r2
      real*8 xhea, yhea, zhea, rij
      integer*4 idya1, iikey

      do i = 1, nboost
         displace(i) = 0.0
      end do

      do i = 1, idya
         idi = idya_n(i)
         xi = x0(idi)
         yi = y0(idi)
         zi = z0(idi)
         xj = xjp(i)
         yj = yjp(i)
         zj = zjp(i)
         xij = (xi - xj)*b0(1,1)
         yij = (yi - yj)*b0(2,2)
         zij = (zi - zj)*b0(3,3)

         if(xij > b0(1,1)/2.0d0) xij = xij - b0(1,1)
         if(xij <-b0(1,1)/2.0d0) xij = xij + b0(1,1)
         if(yij > b0(2,2)/2.0d0) yij = yij - b0(2,2)
         if(yij <-b0(2,2)/2.0d0) yij = yij + b0(2,2)
         if(zij > b0(3,3)/2.0d0) zij = zij - b0(3,3)
         if(zij <-b0(3,3)/2.0d0) zij = zij + b0(3,3)

         dist = sqrt(xij*xij + yij*yij + zij*zij)
         displace(i) = dist
      end do

      return
      end subroutine


text:
subroutine initur(j         ,nmmaxj    ,nmmax     ,kmax      ,icx       , &
                & icy       ,ltur      ,lturi     ,rtur0     , &
                & s0        ,dps       ,hu        ,hv        ,u0        , &
                & v0        ,thick     ,windsu    ,windsv    ,z0urou    , &
                & z0vrou    ,kfu       ,kfv       ,kfs       ,kcs       , &
                & dudz      ,dvdz      ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program 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 General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: initur.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/inichk/initur.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Initialises turbulent energy and/or
!              turbulent dissipation
! Method used: Reference: R.E. Uittenbogaard, J.A.Th.M. van
!              Kester, G.S. Stelling, Implementation of three
!              turbulence models in 3D-FLOW for rectangular
!              grids, Deltares report Z81, april 1992)
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    use dfparall
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp)               , pointer :: rhow
    real(fp)               , pointer :: z0
    real(fp)               , pointer :: z0v
    real(fp)               , pointer :: vonkar
    real(fp)               , pointer :: cde
    real(fp)               , pointer :: cmukep
    real(fp)               , pointer :: cewall
    real(fp)               , pointer :: zwi
    real(fp)               , pointer :: ck
    integer                , pointer :: inpzw
!
! Global variables
!
    integer                                                   , intent(in) :: icx    !!  Increment in the X-dir., if ICX= NMAX then
                                                                                !!  computation proceeds in the X-dir.
                                                                                !!  If icx=1 then computation proceeds in the Y-dir.
    integer                                                   , intent(in) :: icy    !!  Increment in the Y-dir. (see ICX)
    integer                                                                :: j      !!  Begin pointer for arrays which have been
                                                                                !!  transformed into 1D arrays. Due to the shift
                                                                                !!  in the 2nd (M-) index, J = -2*NMAX + 1
    integer                                                   , intent(in) :: kmax   !  Description and declaration in esm_alloc_int.f90
    integer                                                   , intent(in) :: ltur   !  Description and declaration in esm_alloc_int.f90
    integer                                                   , intent(in) :: lturi  !  Description and declaration in tricom.igs
    integer                                                   , intent(in) :: nmmax  !  Description and declaration in dimens.igs
    integer                                                                :: nmmaxj !  Description and declaration in dimens.igs
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: kcs    !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: kfs    !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: kfu    !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: kfv    !  Description and declaration in esm_alloc_int.f90
    real(prec), dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: dps    !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: hu     !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: hv     !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: s0     !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: windsu !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: windsv !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: z0urou !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)              , intent(in) :: z0vrou !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)                   :: dudz   !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)                   :: dvdz   !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax, ltur)             :: rtur0  !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)        , intent(in) :: u0     !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)        , intent(in) :: v0     !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(kmax)                               , intent(in) :: thick  !  Description and declaration in esm_alloc_real.f90
!
! Local variables
!
    integer                               :: k
    integer                               :: kup
    integer                               :: maskval
    integer                               :: ndm
    integer                               :: nm
    integer                               :: nmd
    integer                               :: nm_pos    ! indicating the array to be exchanged has nm index at the 2nd place, e.g., dbodsd(lsedtot,nm)
    real(fp)                              :: dz
    real(fp)                              :: ene
    real(fp)                              :: h0
    real(fp)                              :: rl
    real(fp)                              :: shear
    real(fp)                              :: sigw
    real(fp)                              :: tke
    real(fp)                              :: tsg
    real(fp)                              :: ust
    real(fp)                              :: utot
    real(fp)                              :: uuu
    real(fp)                              :: vvv
    real(fp), dimension(:,:), allocatable :: z0scratch
    real(fp)                              :: z00
    real(fp)                              :: zw
    real(fp)                              :: zwc
!
!! executable statements -------------------------------------------------------
!
    cde         => gdp%gdturcoe%cde
    cmukep      => gdp%gdturcoe%cmukep
    cewall      => gdp%gdturcoe%cewall
    zwi         => gdp%gdturcoe%zwi
    ck          => gdp%gdturcoe%ck
    inpzw       => gdp%gdturcoe%inpzw
    rhow        => gdp%gdphysco%rhow
    z0          => gdp%gdphysco%z0
    z0v         => gdp%gdphysco%z0v
    vonkar      => gdp%gdphysco%vonkar
    !
    nm_pos      =  1
    !
    allocate(z0scratch(gdp%d%nmlb:gdp%d%nmub, kmax))
    !
    nmd = -icx
    ndm = -icy
    do nm = 1, nmmax
       nmd = nmd + 1
       ndm = ndm + 1
       if (kfs(nm)==1) then
          z0scratch(nm, kmax) = (  kfu(nm)*z0urou(nm) + kfu(nmd)*z0urou(nmd)   &
                              &  + kfv(nm)*z0vrou(nm) + kfv(ndm)*z0vrou(ndm) ) &
                              &  /(kfu(nm) + kfu(nmd) + kfv(nm) + kfv(ndm))
       endif
    enddo
    !
    !***PRODUCTION AND BUOYANCY TERM (ONLY VERTICAL GRADIENTS)
    !
    if (lturi>0) then
       do k = 1, kmax-1
          kup = k + 1
          tsg = 0.5*(thick(k) + thick(kup))
          do nm = j, nmmax
             dudz(nm, k) = 0.0
             if (kfu(nm)==1) then
                dudz(nm, k) = (u0(nm, k) - u0(nm, kup))/(hu(nm)*tsg)
             endif
             dvdz(nm, k) = 0.0
             if (kfv(nm)==1) then
                dvdz(nm, k) = (v0(nm, k) - v0(nm, kup))/(hv(nm)*tsg)
             endif
          enddo
       enddo
       !
       !***KOLMOGOROV-PRANDTL MIXING LENGTH MODEL
       !
       sigw = 0.0
       do k = 1, kmax-1
          sigw = sigw - thick(k)
          nmd = -icx
          ndm = -icy
          do nm = 1, nmmax
             nmd = nmd + 1
             ndm = ndm + 1
             if (kfs(nm)==1) then
                h0 = max(0.01_fp, real(dps(nm),fp) + s0(nm))
                if (kcs(nm)==3) then
                   maskval = kcs(nm) - 2
                else
                   maskval = kcs(nm)
                endif
                !
                !***MIXING LENGTH
                !
                rl = vonkar*h0*(1. + sigw)*sqrt( - sigw)
                !
                !***INITIALISATION TKE and EPSILON from PRANDTL MIXING LENGTH MODEL
                !
                shear = maskval*0.5 *(  dudz(nm, k)**2 + dudz(nmd, k)**2   &
                      &               + dvdz(nm, k)**2 + dvdz(ndm, k)**2 )
                ene   = maskval*0.25*(  u0(nm, k)**2 + u0(nmd, k)**2       &
                      &               + v0(nm, k)**2 + v0(ndm, k)**2 )
                tke   = min(ene, rl*rl*shear/sqrt(cmukep))
                rtur0(nm, k, 1) = tke
             endif
          enddo
       enddo
       !
       !***TKE AT FREE SURFACE AND BOTTOM
       !
       nmd = -icx
       ndm = -icy
       do nm = 1, nmmax
          nmd = nmd + 1
          ndm = ndm + 1
          if (kfs(nm)==1) then
             h0 = max(0.01_fp, real(dps(nm),fp) + s0(nm))
             if (kcs(nm)==3) then
                maskval = kcs(nm) - 2
             else
                maskval = kcs(nm)
             endif
             !
             !***AT FREE SURFACE
             !
             tke = sqrt((windsu(nm)**2 + windsv(nm)**2)/cmukep)/rhow
             rtur0(nm, 0, 1) = tke
             !
             !***AT BOTTOM
             !
             z00  = z0scratch(nm, kmax)
             uuu  = 0.5*(u0(nmd, kmax) + u0(nm, kmax))*maskval
             vvv  = 0.5*(v0(ndm, kmax) + v0(nm, kmax))*maskval
             utot = sqrt(uuu*uuu + vvv*vvv)
             dz   = 0.5*thick(kmax)*h0
             ust  = utot*vonkar/log(1. + dz/z00)
             tke  = ust*ust/sqrt(cmukep)
             rtur0(nm, kmax, 1) = tke
          endif
       enddo
    endif
    !
    !***KOLMOGOROV-PRANDTL MIXING LENGTH MODEL FOR EPS
    !
    if (ltur>1) then
       sigw = 0.0
       do k = 1, kmax - 1
          sigw = sigw - thick(k)
          do nm = 1, nmmax
             if (kfs(nm)==1) then
                h0 = max(0.01_fp, real(dps(nm),fp) + s0(nm))
                !
                !***MIXING LENGTH
                !
                rl = vonkar*h0*(1. + sigw)*sqrt( - sigw)
                !
                !***INITIALISATION EPSILON from PRANDTL MIXING LENGTH MODEL
                !
                rtur0(nm, k, 2) = max(1.E-7_fp, cde*(rtur0(nm, k, 1))**1.5/rl)
             endif
          enddo
       enddo
       !
       !***EPS AT FREE SURFACE AND BOTTOM
       !
       do nm = 1, nmmax
          if (kfs(nm)==1) then
             h0 = max(0.01_fp, real(dps(nm),fp) + s0(nm))
             !
             !***AT FREE SURFACE
             !
             zwc = 0.5*thick(1)*h0
             zw  = inpzw*zwi + (1 - inpzw)*zwc
             rtur0(nm, 0, 2) = max(1.E-7_fp, cewall*(rtur0(nm, 0, 1))**1.5/zw)
             !
             !***AT BOTTOM
             !
             z00 = z0scratch(nm, kmax)
             rtur0(nm, kmax, 2) = max(1.E-7_fp, cewall*(rtur0(nm, kmax, 1))**1.5/z00)
          endif
       enddo
    endif
    !
    if (lturi>0) then
       do k = 1, kmax
          do nm = 1, nmmax
             rtur0(nm, k, 1) = max(1.E-7_fp, rtur0(nm, k, 1))
          enddo
       enddo
    endif
    !
    ! For domain decomposition (in couple points)
    do k = 1, kmax
       do nm = 1, nmmax
          if (kcs(nm)==3) then
             rtur0(nm, k, 1) = 1.E-7
             rtur0(nm, k, 2) = 1.E-7
          endif
       enddo
    enddo
    !
    deallocate(z0scratch)
    !
    ! Delft3D-16494: NOT NECESSARY? No implicit loop was performed
    ! Exchange turbulence values with neighbours for parallel runs
    !
    call dfexchg ( rtur0(:,:,1), 0, kmax, dfloat, nm_pos, gdp )
    call dfexchg ( rtur0(:,:,2), 0, kmax, dfloat, nm_pos, gdp )
end subroutine initur

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


module functions
    contains

subroutine spaceevenly(start, finish, number, arrayout)
    implicit none

    integer, intent(in) :: start, finish
    real(kind=8), intent(in) :: number
    real(kind=8) :: div
    real(kind=8), dimension(0:finish), intent(out) :: arrayout
    integer :: i

    div = number / ((finish-1) - (start-1))

    do i=start, finish
        arrayout(i) = div * i
    end do

end subroutine spaceevenly

subroutine diffuse(timesteps, u, nu, dt, dx, dy, nx, ny)
    implicit none
    integer, intent(in) :: nx, ny, timesteps
    real(kind=8), intent(in) :: nu, dt, dx, dy
    real(kind=8), intent(out), dimension(nx,ny) :: u
    real(kind=8), dimension(nx,ny) :: un
    integer :: i, j, time
    real(kind=8) :: Ydisip, Xdisip, Yderiv, Xderiv
    time = timesteps

    do time = 1, time
            un = u

            do j=2, nx-1
                do i=2, ny-1
                    Xderiv = (nu * dt) / (dx**2)
                    Yderiv = (nu * dt) / (dy**2)
                    Xdisip = un(j,i+1) - (2 * un(j,i)) + un(j,i-1)
                    Ydisip = un(j+1,i) - (2 * un(j,i)) + un(j-1,i)

                    u(j,i) = un(j,i) + (Xderiv * Xdisip) + (Yderiv * Ydisip)

                    u(1,j) = 1
                    u(nx,j) = 1
                    u(i,1) = 1
                    u(i,ny) = 1

                    u(1,1) = 1
                    u(1,nx) = 1
                    u(ny,1) = 1
                    u(ny,nx) = 1
                end do
            end do
        end do
end subroutine

end module functions


text:
Function rlu(idum)
  Common /ludatr/mrlu(6), rrlu(100)
  Save /ludatr/
  Equivalence (mrlu1, mrlu(1)), (mrlu2, mrlu(2)), (mrlu3, mrlu(3)), (mrlu4, mrlu(4)), (mrlu5, mrlu(5)), (mrlu6, mrlu(6)), (rrlu98, rrlu(98)), (rrlu99, rrlu(99)), (rrlu00, rrlu(100))
  If (mrlu2==0) Then
    ij = mod(mrlu1/30082, 31329)
    kl = mod(mrlu1, 30082)
    i = mod(ij/177, 177) + 2
    j = mod(ij, 177) + 2
    k = mod(kl/169, 178) + 1
    l = mod(kl, 169)
    Do ii = 1, 97
      s = 0.
      t = 0.5
      Do jj = 1, 24
        m = mod(mod(i*j,179)*k, 179)
        i = j
        j = k
        k = m
        l = mod(53*l+1, 169)
        If (mod(l*m,64)>=32) s = s + t
        t = 0.5*t
      End Do
      rrlu(ii) = s
    End Do
    twom24 = 1.
    Do i24 = 1, 24
      twom24 = 0.5*twom24
    End Do
    rrlu98 = 362436.*twom24
    rrlu99 = 7654321.*twom24
    rrlu00 = 16777213.*twom24
    mrlu2 = 1
    mrlu3 = 0
    mrlu4 = 97
    mrlu5 = 33
  End If
  130 runi = rrlu(mrlu4) - rrlu(mrlu5)
  If (runi<0.) runi = runi + 1.
  rrlu(mrlu4) = runi
  mrlu4 = mrlu4 - 1
  If (mrlu4==0) mrlu4 = 97
  mrlu5 = mrlu5 - 1
  If (mrlu5==0) mrlu5 = 97
  rrlu98 = rrlu98 - rrlu99
  If (rrlu98<0.) rrlu98 = rrlu98 + rrlu00
  runi = runi - rrlu98
  If (runi<0.) runi = runi + 1.
  If (runi<=0 .Or. runi>=1.) Goto 130
  mrlu3 = mrlu3 + 1
  If (mrlu3==1000000000) Then
    mrlu2 = mrlu2 + 1
    mrlu3 = 0
  End If
  rlu = runi
  Return
End Function rlu

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
      program t
      integer::i
      character(5)::s
      open(unit=77,file='tmpfile',status='new',form='unformatted')
      write(77) 77
      rewind(77)
      read(77) i
      print '(i0)',i
      close (77,status='delete')
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
sharonstahl can describe herself here.

Hello and welcome to the wiki/business owner welcome to the Wiki. Please check out the link to the business owner welcome page. Thank you.Users/hankim

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! { dg-do compile }

      SUBROUTINE A21_GOOD(N)
      INTEGER N
!$OMP DO ORDERED
        DO I = 1,N
          IF (I <= 10) THEN
!$OMP ORDERED
              CALL WORK(I)
!$OMP END ORDERED
          ENDIF
          IF (I > 10) THEN
!$OMP ORDERED
              CALL WORK(I+1)
!$OMP END ORDERED
          ENDIF
        ENDDO
      END SUBROUTINE A21_GOOD

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      Subroutine gravity(lat,g)
c       calculates g as a funciton of latitude using the 1980 IUGG formula
c
c       Bulletin Geodesique, Vol 62, No 3, 1988 (Geodesist's Handbook)
c       p 356, 1980 Gravity Formula (IUGG, H. Moritz)
c       units are in m/sec^2 and have a relative precision of 1 part
c       in 10^10 (0.1 microGal)
c       code by M. Zumberge.
c
c       check values are:
c
c        g = 9.780326772 at latitude  0.0
c        g = 9.806199203 at latitude 45.0
c        g = 9.832186368 at latitude 90.0
c
      real*8 gamma, c1, c2, c3, c4, phi, lat, g

      gamma = 9.7803267715
      c1 = 0.0052790414
      c2 = 0.0000232718
      c3 = 0.0000001262
      c4 = 0.0000000007
      phi = lat * 3.14159265358979 / 180.0
      g = gamma * (1.0
     $ + c1 * ((sin(phi))**2)
     $ + c2 * ((sin(phi))**4)
     $ + c3 * ((sin(phi))**6)
     $ + c4 * ((sin(phi))**8))
c
      return
      end


program special_characters
!! This program shows a few special ASCII characters in Fortran.
!! https://en.wikipedia.org/wiki/ASCII#Character_groups

use iso_c_binding, only: c_null_char

implicit none (type, external)

character, parameter :: &
nul = c_null_char, &
etx = char(3), &
tab = char(9), &
backslash = char(92)   ! necessary for strict compilers like PGI and Flang in strings

print *, 'null'
print '(A1)', nul

print *, 'etx'
print '(A1)', etx

print *, 'tab'
print '(A1)', tab

print *, 'backslash'
print '(A1)', backslash

end program


text:
! { dg-do compile }
! Tests a number of error messages relating to derived type objects
! in common blocks.  Originally due to PR 33198

subroutine one
type a
   sequence
   integer :: i = 1
end type a
type(a) :: t ! { dg-error "Derived type variable .t. in COMMON at ... may not have default initializer" }
common /c/ t
end

subroutine first
type a
   integer :: i
   integer :: j
end type a
type(a) :: t  ! { dg-error "Derived type variable .t. in COMMON at ... has neither the SEQUENCE nor the BIND.C. attribute" }
common /c/ t
end

subroutine prime
type a
   sequence
   integer, allocatable :: i(:)
   integer :: j
end type a
type(a) :: t  ! { dg-error "Derived type variable .t. in COMMON at ... has an ultimate component that is allocatable" }
common /c/ t
end

subroutine source
parameter(x=0.) ! { dg-error "COMMON block .x. at ... is used as PARAMETER at ..." }
common /x/ i  ! { dg-error "COMMON block .x. at ... is used as PARAMETER at ..." }
intrinsic sin
common /sin/ j ! { dg-error "COMMON block .sin. at ... is also an intrinsic procedure" }
end subroutine source


text:
	SUBROUTINE TI_C2I  ( dattim, iyyyy, immdd, ihhmm, iret )
C************************************************************************
C* TI_C2I								*
C*									*
C* This subroutine converts a standard GEMPAK time into three		*
C* integers representing the year, month and day, and hour and		*
C* minute.								*
C*									*
C* Example:  DATTIM = 840317/1200 					*
C*           IYYYY = 1984  IMMDD = 317  IHMM = 1200			*
C*									*
C* TI_C2I  ( DATTIM, IYYYY, IMMDD, IHHMM, IRET )			*
C*									*
C* Input parameters:							*
C*	DATTIM		CHAR*		GEMPAK time 			*
C*									*
C* Output parameters:							*
C*	IYYYY		INTEGER 	Year  (YYYY)			*
C*	IMMDD		INTEGER 	Month and day  (MMDD)		*
C*	IHHMM		INTEGER 	Hour and minute  (HHMM)		*
C*	IRET		INTEGER 	Return code			*
C*					  0 = normal return		*
C*					 -1 = invalid time		*
C*					 -7 = invalid year		*
C*					 -8 = invalid month		*
C*					 -9 = invalid day		*
C*					-10 = invalid hour		*
C*					-11 = invalid minute		*
C**									*
C* Log:									*
C* M. Goodman/RDS	 4/84	Original source code			*
C* M. desJardins/GSFC	 6/84						*
C************************************************************************
	CHARACTER*(*) 	dattim
	INTEGER		idtarr (5)
C-----------------------------------------------------------------------
	iret= 0
C
C*	Break into 5 integers.
C
	CALL TI_CTOI  ( dattim, idtarr, iret )
	IF  ( iret .ne. 0 )  THEN
	    iyyyy = 0
	    immdd = 0
	    ihhmm = 0
	    RETURN
	END IF
C
C*	Form output variables
C
	iyyyy = idtarr (1)
	immdd = idtarr (2) * 100 + idtarr (3)
	ihhmm = idtarr (4) * 100 + idtarr (5)
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak