text:
module GMS_mod_spectrum

    !
    ! wavy - A spectral ocean wave modeling and development framework
    ! Copyright (c) 2017, Wavebit Scientific LLC
    ! All rights reserved.
    !
    ! Licensed under the BSD-3 clause license. See LICENSE for details.
    ! Modified by Bernard Gingold (beniekg@gmail.com) on 16/03/2019
    !

    use ,intrinsic :: ISO_C_BINDING, Only :  c_int, c_double
    use GMS_mod_kinds,     only : int32_t, dp
    use GMS_mod_utility,   only : diff, diff_periodic
    use GMS_mod_datetime,  only : datetime
    use GMS_mod_timedelta, only : timedelta
    implicit none

    private

    type, public :: spectrum_type

          public

          type(datetime)  :: start_time
          type(datetime)  :: end_time
          type(timedelta) :: time_step
!DIR$   ATTRIBUTES ALIGN : 64 :: spec
          real(kind=dp), allocatable, dimension(:,:) :: spec ! 2-d spectrum
!DIR$   ATTRIBUTES ALIGN : 64 :: f
          real(kind=dp), allocatable, dimension(:)   :: f    ! Frequency [Hz]
!DIR$   ATTRIBUTES ALIGN : 64 :: df
          real(kind=dp), allocatable, dimension(:)   :: df   ! Frequency spacing [Hz]
!DIR$   ATTRIBUTES ALIGN : 64 :: k
          real(kind=dp), allocatable, dimension(:)   :: k    ! Wavenumber [rad/m]
!DIR$   ATTRIBUTES ALIGN : 64 :: dk
          real(kind=dp), allocatable, dimension(:)   :: dk   ! Wavenumber spacing [rad/m]
!DIR$   ATTRIBUTES ALIGN : 64 ::  th
          real(kind=dp), allocatable, dimension(:)   :: th   ! Direction [rad]
!DIR$   ATTRIBUTES ALIGN : 64 :: dth
          real(kind=dp), allocatable, dimension(:)   :: dth  ! Directional spacing [rad]
!DIR$   ATTRIBUTES ALIGN : 64 :: cp
          real(kind=dp), allocatable, dimension(:)   :: cp   ! Phase speed [m/s]
!DIR$   ATTRIBUTES ALIGN : 64 :: cg
          real(kind=dp), allocatable, dimension(:)   :: cg   ! Group speed [m/s]
!DIR$   ATTRIBUTES ALIGN : 64 :: u
          real(kind=dp), allocatable, dimension(:)   :: u    ! Mean current velocity in x-direction [m/s]
!DIR$   ATTRIBUTES ALIGN : 64 :: v
          real(kind=dp), allocatable, dimension(:)   :: v    !  Mean current velocity in y-direction [m/s]
!DIR$   ATTRIBUTES ALIGN : 64 :: z
          real(kind=dp), allocatable, dimension(:)   :: z    ! Depth levels for current array [m]
          real(kind=dp) :: air_density     !! Air density [kg/m^3]
          real(kind=dp) :: depth           !! Mean water depth [m]
          real(kind=dp) :: elevation       !! Mean surface elevation [m]
          real(kind=dp) :: grav            !! Gravitational acceleration [m/s^2]
          real(kind=dp) :: surface_tension !! Surface tension [N/m]
          real(kind=dp) :: water_density   !! Water density [kg/m^3]

    end type spectrum_type

    contains

    subroutine constructor(spectrum,fmin,fmax,df,ndirs,depth,grav,    &
                 air_density,water_density,surface_tension,           &
                 errstate,iounit,logging,verbose,append,fname)
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: constructor
          use GMS_mod_print_error,  only : handle_fatal_memory_error,  &
                                       print_non_fatal_error
          type(spectrum_type),      intent(inout) :: spectrum
          real(kind=dp),            intent(in)    :: fmin
          real(kind=dp),            intent(in)    :: fmax
          real(kind=dp),            intent(in)    :: df
          integer(kind=int32_t),    intent(in)    :: ndirs
          real(kind=dp),            intent(in)    :: depth
          real(kind=dp),            intent(in), optional :: grav
          real(kind=dp),            intent(in), optional :: air_density
          real(kind=dp),            intent(in), optional :: water_density
          real(kind=dp),            intent(in), optional :: surface_tension
          logical(kind=int32_t),    intent(inout) :: errstate
          integer(kind=int32_t),    intent(in)    :: iounit
          logical(kind=int32_t),    intent(in)    :: logging
          logical(kind=int32_t),    intent(in)    :: verbose
          logical(kind=int32_t),    intent(in)    :: append
          character(len=*),         intent(in)    :: fname
          ! Locals
!DIR$   ATTRIBUTES ALIGN : 64 :: tmp_f,tmp_k,th_tmp
          real(kind=dp), allocatable, dimension(:) :: tmp_f,tmp_k,th_tmp
          character(len=256),    automatic :: emsg

          integer(kind=int32_t), automatic :: aerr
          integer(kind=int32_t), automatic :: n
          integer(kind=int32_t), automatic :: nfreqs
          ! EXec code.....
          if(present(grav)) then
                spectrum.grav = grav
          else
                spectrum.grav = 9.8_dp
          end if
          if(present(air_density)) then
                spectrum.air_density = air_density
          else
                spectrum.air_density = 1.2_dp
          end if
          if(present(water_density)) then
                spectrum.water_density = water_density
          else
                spectrum.water_density = 1.0e3_dp
          end if
          if(present(surface_tension)) then
                spectrum.surface_tension = surface_tension
          else
                spectrum.surface_tension = 0.07_dp
          end if
          spectrum.depth = depth
          if(fmin == fmax) then
                nfreqs = 1
          else
                nfreqs = int((log(fmax)-log(fmin))/log(df))
          end if
          if(allocated(spectrum.spec)) then
                deallocate(spectrum.spec)
                allocate(spectrum.spec(nfreqs,ndirs),  &
                            STAT=aerr,ERRMSG=emsg)
                if(aerr /= 0) goto 9999
          else
                allocate(spectrum.spec(nfreqs,ndirs),  &
                            STAT=aerr,ERRMSG=emsg)
                if(aerr /= 0) goto 9999
          end if
          spectrum.spec = 0.0_dp
          if(allocated(spectrum.f)) then
                 deallocate(spectrum.f)
                 allocate(spectrum.f(nfreqs),  &
                             STAT=aerr,ERRMSG=emsg)
                  if(aerr /= 0) goto 9999
          else
                  allocate(spectrum.f(nfreqs),  &
                             STAT=aerr,ERRMSG=emsg)
                  if(aerr /= 0) goto 9999
           end if
           if(allocated(spectrum.df)) then
                   deallocate(spectrum.df)
                   allocate(spectrum.df(nfreqs), &
                               STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           else
                   allocate(spectrum.df(nfreqs), &
                               STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           end if
           if(allocated(spectrum.k)) then
                   deallocate(spectrum.k))
                   allocate(spectrum.k(nfreqs), &
                               STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           else
                   allocate(spectrum.k(nfreqs), &
                               STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           end if
           if(allocated(spectrum.dk)) then
                   deallocate(spectrum.dk)
                   allocate(spectrum.dk(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           else
                   allocate(spectrum.dk(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           end if
           if(allocated(spectrum.cp)) then
                   deallocate(spectrum.cp)
                   allocate(spectrum.cp(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           else
                   allocate(spectrum.cp(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           end if
           if(allocated(spectrum.cg)) then
                   deallocate(spectrum.cg)
                   allocate(spectrum.cg(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           else
                   allocate(spectrum.cg(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           end if
           if(allocated(spectrum.th)) then
                   deallocate(spectrum.th)
                   allocate(spectrum.th(ndirs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           else
                   allocate(spectrum.th(nfreqs), &
                                STAT=aerr,ERRMSG=emsg)
                   if(aerr /= 0) goto 9999
           end if
           allocate(tmp_f(nfreqs))
           allocate(tmp_k(nfreqs))
           tmp_f = 0.0_dp; tmp_k = 0.0_dp
           if(nfreqs == 1) then
                spectrum.f(1) = fmin
           else
                do concurrent(n=1:nfreqs)
                    spectrum.f(n) = exp(log(fmin)+(n-1)*log(df))
                end do
           end if
           do concurrent(n=1:nfreqs)
                spectrum.k(n) = wavenumber(spectrum.f(n),          &
                                           spectrum.depth,         &
                                           spectrum.water_density, &
                                           spectrum.grav,          &
                                           spectrum.surface_tension)
           end do
           spectrum.cp = pi2_const*spectrum.f/spectrum.k
           call diff(spectrum.f,tmp_f,nfreqs)
           call diff(spectrum.k,tmp_k,nfreqs)
           spectrum.cg = pi2_const*tmp_f/tmp_k

           do concurrent(n=1:ndirs)
                spectrum.th(n) = (n-0.5_dp*(ndirs+1))*2pi_const/ndirs
           end do
           spectrum.df = tmp_f
           spectrum.dk = tmp_k
           if(ndirs > 1) then
                allocate(th_tmp(ndirs))
                call diff_periodic(spectrum.th.th_tmp,ndirs)
                spectrum.dth = th_tmp
                spectrum.dth(1) = spectrum.dth(2)
                spectrum.dth(ndirs) = spectrum.dth(ndirs-1)
           else
                spectrum.dth = [1]
           end if
           call setCurrent2d(spectrum,[0.0_dp],[0.0_dp],[0.0_dp])
           return
9999
           call handle_fatal_memory_error(iounit, logging,verbose,append,fname, &
                     "logger: "// __FILE__ // "module: mod_spectrum, subroutine: constructor -- Memory Allocation Failure !!", &                                                        &
                              "module: mod_spectrum, subroutine: constructor -- Memory Allocation Failure !!", &
                              emsg,243)
    end subroutine constructor

!DIR$   ATTRIBUTES INLINE :: isAllocated
    pure elemental logical function isAllocated(spectrum)
!! Returns the allocation status of the spectrum array.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: isAllocated
          type(spectrum_type),      intent(in) :: spectrum
          ! Exec code ....
          isAllocated = allocated(spectrum.spec)
    end function isAllocated

!DIR$   ATTRIBUTES INLINE :: isMonochromatic
    pure elemental function isMonochromatic(spectrum) result(monochromatic)
  !! Returns `.true.` if only one frequency bin is allocated,
  !! and `.false.` otherwise.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: isMonochromatic
          type(spectrum_type),      intent(in) :: spectrum
          ! Locals
          logical(kind=int32_t) :: monochromatic
          ! Exec code ....
          if(size(spectrum.f) == 1) then
                monochromatic = .true.
          else
                monochromatic = .false.
          end if
    end function isMonochromatic

!DIR$   ATTRIBUTES INLINE :: isOmnidirectional
    pure elemental function isOmnidirectional(spectrum) result(omnidirect)
 !! Returns `.true.` if only one direction bin is allocated,
  !! and `.false.` otherwise.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: isOmnidirectional
          type(spectrum_type),      intent(in) :: spectrum
          ! LOcals
          logical(kind=int32_t) :: omnidirect
          ! Exec code ...
          if(size(spectrum.th) == 1) then
                omnidirect = .true.
          else
                omnidirect = .false.
          end if
    end function isOmnidirectional

    subroutine getFrequency2d(spectrum,f_out,nfreqs,ndirs)
  !!Extract frequency [Hz] array of the spectrum instance, reshaped to
  !! match the spectrum array shape. This method is most useful for conforming
  !! shape array in 2-d spectrum computations.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getFRequency2d
          type(spectrum_type),      intent(in) :: spectrum
 !DIR$  ASSUME_ALIGNED f_out:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: f_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t), automatic :: ndir
          do concurrent(ndir=1:ndirs)
                f_out(:,ndir) = spectrum.f
          end do
    end subroutine getFrequency2d

    subroutine getWavenumber2d(spectrum,k_out,nfreqs,ndirs)
!! Returns the wavenumber [rad/m] array of the spectrum instance, reshaped to
  !! match the spectrum array shape. This method is most useful for conforming
  !! shape array in 2-d spectrum computations.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getWavenumber2d
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED k_out:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: k_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t), automatic :: ndir
          ! EXec code ....
          do concurrent(ndir=1:ndirs)
                k_out(:,ndir) = spectrum.k
          end do
    end subroutine getWavenumber2d

!DIR$   ATTRIBUTES INLINE :: getWavelength
    subroutine getWavelength(spectrum,lambda,nfreqs)
          use mod_constants, only : pi2_const
!! Returns the wavelength [m] array of the spectrum instance.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getWavelength
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED lambda:64
          real(kind=dp),    dimension(nfreqs), intent(out) :: lambda
          integer(kind=int32_t),    intent(in) ::nfreqs
          ! Exec code ....
          lambda = pi2_const/spectrum.k
    end subroutine getWavelength

    subroutine getDirection2d(spectrum,th_out,nfreqs,ndirs)
 !! Returns the directions [rad] array of the spectrum instance, reshaped to
  !! match the spectrum array shape. This method is most useful for conforming
  !! shape array in 2-d spectrum computations.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getDirection2d
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED th_out:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: th_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t), automatic :: ndir
          ! EXec code ....
          do concurrent(ndir=1:ndirs)
                th_out(nfreq,:) = spectrum.th
          end do
    end subroutine getDirection2d

    subroutine getPhaseSpeed2d(spectrum,cp_out,nfreqs,ndirs)
  !! Returns the phase speed [m/s] array of the spectrum instance.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getPhaseSpeed2d
          type(spectrum_type),      intent(in) :: spectrum
!DIR$ ASSUME_ALIGNED cp_out:64
          real(kind=dp),    dimension(nfreqs,ndirs),  intent(out) :: cp_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t), automatic :: ndir
          ! EXec code ....
          do concurrent(ndir=1:ndirs)
                cp_out(:,ndir) = spectrum.cp
          end do
    end subroutine getPhaseSpeed2d

    subroutine getGroupSpeed2d(spectrum,cg_out,nfreqs,ndirs)
 !! Returns the group speed [m/s] array of the spectrum instance.
 !DIR$  ATTRIBUTES CODE_ALIGN:32 :: getGroupSpeed2d
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED cg:64
          real(kind=dp), dimension(nfreqs,ndirs),   intent(out) :: cg_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t), automatic :: ndir
          ! EXec code ....
          do concurrent(ndir=1:ndirs)
                cg_out(:,ndir) = spectrum.cg
          end do
    end subroutine getGroupSpeed2d

!DIR$ ATTRIBUTES INLINE :: getWaveAction
    subroutine getWaveAction(spectrum,wave_action,freq2d,nfreqs,ndirs)
!! Returns the wave action spectrum, which corresponds to the the wave
  !! variance spectrum normalized by the intrinsic frequency.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getWaveAction
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED wave_action:64
          real(kind=dp), dimension(nfreqs,ndirs), intent(in) :: wave_action
!DIR$   ASSUME_ALIGNED freq2d:64
          real(kind=dp), dimension(nfreqs,ndirs), intent(in) :: freq2d
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Exec code
          wave_action = spectrum.spec/freq2d
    end subroutine getWaveAction

!DIR$   ATTRIBUTES INLINE :: getAmplitude
    subroutine getAmplitude(spectrum,a,nfreqs,ndirs)
  !! Returns the amplitude array.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: getAmplitude
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED a:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: a
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t) :: ndir,ndirs2
          ! Exec code ...
          if(isMonochromatic(spectrum)) then
                a = sqrt(spectrum.spec+spectrum.spec)
          else
                ndirs2 = size(spectrum.th)
                do concurrent(ndir=1:ndirs2)
                         a(:,ndir) = sqrt((spectrum.spec(:,ndir)+spectrum.spec(:,ndir)) * &
                                           spectrum.df)
                end do
          end if

    end subroutine getAmplitude

!DIR$   ATTRIBUTES INLINE :: omnidirectionalSpectrum
    subroutine omnidirectionalSpectrum(spectrum,spec_out,nfreqs,ndirs)
 !! Returns the omnidirectional spectrum that corresponds to the input
  !! directional spectrum, integrated over all directions.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: omnidirectionalSpectrum
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_out:64
          real(kind=dp),    dimension(nfreqs),  intent(out) :: spec_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t) :: ndir
!DIR$   VECTOR ALIGNED
!DIR$   SIMD
          do ndir = 1, ndirs
                spec_out(:) = spec_out(:)+spectrum.spec(:,ndir)*spectrum.dth(ndir)
          end do
    end subroutine omnidirectionalSpectrum

    subroutine setCurrent1d(spectrum,u,z)
!! Sets the 1-d current velocity field.
!DIR$   ATTRIBUTES CODE_ALIGN:32 ::  setCurrent1d
          type(spectrum_type),      intent(inout) :: spectrum
!DIR$   ASSUME_ALIGNED u:64
          real(kind=dp),    dimension(:), intent(in) :: u
!DIR$   ASSUME_ALIGNED z:64
          real(kind=dp),    dimension(:), intent(in) :: z
          ! Exec code ....
          spectrum.u = u
          spectrum.z = z
          allocate(spectrum.v(size(u))
          spectrum.v = 0.0_dp
    end subroutine setCurrent1d

    subroutine setCurrent2d(spectrum,u,v,z)
!! Sets the 2-d current velocity field.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: setCurrent2d
          type(spectrum_type),      intent(inout) :: spectrum
!DIR$   ASSUME_ALIGNED u:64
          real(kind=dp),    dimension(:), intent(in) :: u
!DIR$   ASSUME_ALIGNED v:64
          real(kind=dp),    dimension(:), intent(in) :: v
!DIR$   ASSUME_ALIGNED z:64
          real(kind=dp),    dimension(:), intent(in) :: z
          ! EXec code ....
          spectrum.u = u
          spectrum.v = v
          spectrum.z = z
    end subroutine setCurrent2d

!DIR$   ATTRIBUTES INLINE :: setDepth
    subroutine setDepth(spectrum,depth)
 !! Sets the mean surface elevation value.
 !DIR$  ATTRIBUTES CODE_ALIGN:32 :: setDepth
          type(spectrum_type),  intent(inout) :: spectrum
          real(kind=dp),        intent(in)    :: depth
          ! Exec code ...
          spectrum.depth = depth
    end subroutine setDepth

!DIR$   ATTRIBUTES INLINE :: setElevation
    subroutine setElevation(spectrum,elevation)
 !! Sets the mean surface elevation value.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: setElevation
          type(spectrum_type),   intent(inout) :: spectrum
          real(kind=dp),         intent(in)    :: elevation
          ! Exec code ...
          spectrum.elevation = elevation
    end subroutine setElevation

!DIR$   ATTRIBUTES INLINE :: setGravity
    subroutine setGravity(spectrum,gravity)
!! Sets the gravitational acceleration [m/s^2].
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: setGravity
          type(spectrum_type),    intent(inout) :: spectrum
          real(kind=dp),          intent(in)    :: gravity
          ! Exec code ....
          spectrum.gravity = gravity
    end subroutine setGravity

!DIR$   ATTRIBUTES INLINE :: setSurfaceTension
    subroutine setSurfaceTension(spectrum,surface_tension)
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: setSurfaceTension
          type(spectrum_type),  intent(inout) :: spectrum
          real(kind=dp),        intent(in)    :: surface_tension
          ! Exec code ...
          spectrum.surface_tension = surface_tension
    end subroutine setSurfaceTension

!DIR$   ATTRIBUTES INLINE :: setAirDensity
    subroutine setAirDensity(spectrum,air_density)
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: setAirDensity
          type(spectrum_type),  intent(inout) :: spectrum
          real(kind=dp),        intent(in)    :: air_density
          ! Exec code ....
          spectrum.air_density = air_density
    end subroutine setAirDensity

!DIR$   ATTRIBUTES INLINE :: setWaterDensity
    subroutine setWaterDensity(spectrum,water_density)
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: setWaterDensity
          type(spectrum_type),      intent(inout) :: spectrum
          real(kind=dp),            intent(in)    :: water_density
          ! Exec code ...
          spectrum.water_density = water_density
    end subroutine setWaterDensity

    subroutine meanSquareSlopeDirectional(spectrum,mss,wavenum_spec,dir_proj,  &
                                          nfreqs,ndirs)
!! For each directional frequency bin, computes the mean square slope of all
  !! all waves longer than that bin, projected to the direction of that bin.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: meanSquareSlopeDirectional
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED mss:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: mss
!DIR$   ASSUME_ALIGNED wavenum_spec:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(in)  :: wavenum_spec
!DIR$   ASSUME_ALIGNED dir_proj:64
          real(kind=dp),    dimension(ndirs,ndirs),  intent(in)  :: dir_proj
          integer(kind=int32_t),                     intent(in)  :: nfreqs
          integer(kind=int32_t),                     intent(in)  :: ndirs
          ! Locals
          integer(kind=int32_t), automatic :: ndir
          integer(kind=int32_t), automatic :: nfreq
          ! EXec code ....
          do concurrent(ndir=1:ndirs)
                dir_proj(:,ndir) = abs(cos(spectrum.th(ndir)-spectrum.th(ndir)))
          end do
          do ndir=1, ndirs
!DIR$   SIMD VECTORLENGTHFOR(REAL(KIND=dp))
             do nfreq=2,nfreqs
                 mss(nfreq,ndir) = mss(nfreq-1,ndir)                &
                 + sum(wavenum_spec(nfreq-1,:)*dir_proj(:,ndir)     &
                 * spectrum.k(nfreq-1)**2*spectrum.dk(nfreq-1)
             end do
          end do
    end subroutine meanSquareSlopeDirectional

!DIR$   ATTRIBUTES INLINE :: momentum_x
    pure elemental real(kind=dp) function momentum_x(spectrum)
!! Returns total wave momentum [kg/m/s] in x-direction.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: momentum_x
          type(spectrum_type),  intent(in) :: spectrum
          ! LOcals
          integer(kind=int32_t), automatic :: n,nfreqs
          ! Exec code ....
          nfreqs = size(spectrum.f)
          momentum_x = 0.0_dp
!DIR$ SIMD REDUCTION(+:momentum_x)
          do n=1, nfreqs
              momentum_x = momentum_x+sum(spectrum.spec(n,:)*spectrum.dth*cos(spectrum.th)) &
                * spectrum.df(n)/spectrum.cp(n)
          end do
          momentum_x = mometum_x*spectrum.water_density*spectrum.grav
    end function momentum_x

!DIR$   ATTRIBUTES INLINE :: momentum_y
    pure elemental real(kind=dp) function momentum_y(spectrum)
 !! Returns total wave momentum [kg/m/s] in y-direction.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: momentum_y
          type(spectrum_type),  intent(in) :: spectrum
          ! Locals
          integer(kind=int32_t), automatic :: n,nfreqs
          ! Exec code ....
          nfreqs = size(spectrum.f)
          momentum_y = 0.0_dp
!DIR$   SIMD REDUCTION(+:momentum_y)
          do n=1, nfreqs
               momentum_y = momentum_y+sum(spectrum.spec(n,:)*spectrum.dth*sin(spectrum.th)) &
                            * spectrum.df(n)/spectrum.cp(n)
          end do
          momentum_y = momentum_y*spectrum.water_density*spectrum.grav
    end function  momentum_y

!DIR$   ATTRIBUTES INLINE :: momentumFlux_xx
    pure elemental real(kind=dp) function momentumFlux_xx(spectrum)
!! Returns total advective flux [kg/m^2/s^2] in y-direction of momentum in
  !! y-direction.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: momentumFlux_xx
          type(spectrum_type),    intent(in) :: spectrum
          ! Locals
          integer(kind=int32_t), automatic :: n.nfreqs
          ! Exec code ....
          nfreqs = size(spectrum.f)
          momentumFlux_xx = 0.0_dp
!DIR$   SIMD REDUCTION(+:momentumFlux_xx)
          do n=1, nfreqs
               momentumFlux_xx = momentumFlux_xx            &
                + sum(spectrum.spec(n,:)*spectrum.dth*cos(spectrum.th)**2) &
                * spectrum.df(n)*spectrum.cg(n)/spectrum.cp(n)
          end do
          momentumFlux_xx = momentumFlux_xx*spectrum.water_density*spectrum.grav
    end function momentumFlux_xx

!DIR$   ATTRIBUTES INLINE :: momentumFlux_xy
    pure elemental real(kind=dp) function momentumFlux_xy(spectrum)
 !! Returns total advective flux [kg/m^2/s^2] in x-direction of momentum in
  !! y-direction and vice versa (flux in y-direction of momentum in
  !! y-direction), because \int{Cgx*My} == \int{Cgy*Mx}.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: momentumFlux_xy
          type(spectrum_type),      intent(in) :: spectrum
          ! Locals
          integer(kind=int32_t), automatic :: n,nfreqs
          ! Exec code ....
          nfreqs = size(spectrum.f)
          momentumFlux_xy = 0.0_dp
!DIR$   SIMD  REDUCTION(+:momentumFlux_xy)
          do n=1, nfreqs
              momentumFlux_xy = momentumFlux_xy             &
              + sum(spectrum.spec(n,:)*spectrum.dth*cos(spectrum.th)*sin(spectrum.th) &
              * spectrum.df(n)*spectrum.cg(n)/spectrum.cp(n)
          end do
          momentumFlux_xy = momentumFlux_xy*spectrum.water_density*spectrum.grav
    end function momentumFlux_xy

!DIR$   ATTRIBUTES INLINE :: frequencyMoment
    pure elemental real(kind=dp) function frequencyMoment(spectrum,n)
 !! Returns the spectral frequency moment of order n.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: frequencyMoment
          type(spectrum_type),   intent(in) :: spectrum
          integer(kind=int32_t), intent(in) :: n
          ! Exec code ....
          frequencyMoment = sum(spectrum.f**n*sum(spectrum.spec,dim=2)*spectrum.df)
    end function frequencyMoment

!DIR$   ATTRIBUTES INLINE :: peakedness
    pure elemental real(kind=dp) function peakedness(spectrum,spec_in,ndirs)
  !! Returns the peakedness parameter that quantifies the sharpness of the
  !! spectral peak, following Goda (1970).
  !!
  !! References:
  !!
  !! Goda, Y., 1970. Numerical experiments on waves statistics with spectral
  !! simulation. *Report. Port and Harbour Research Institute*, Japan, **9**,
  !! 3-57.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: peakedness
          use GMS_mod_constants, only : eps
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_in:64
          real(kind=dp), dimension(ndirs), intent(in) :: spec_in
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Exec code ....
          peakedness = 2*sum(spectrum.f*spec_in**2*spectrum.df) &
                       / (frequencyMoment(spectrum,0)**2+eps)
    end function peakedness

!DIR$   ATTRIBUTES INLINE :: peakFrequency
    pure elemental real(kind=dp) function peakFrequency(spectrum,spec_in,ndirs)
 !! Returns the peak frequency based on Young (1995).
  !!
  !! References:
  !!
  !! Young, I, 1995. The determination of confidence limits associated with
  !! estimates of the spectral peak frequency. *Ocean Engng.*, **22**, 669-686.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: peakFrequency
          use GMS_mod_constants, only : eps
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_in:64
          real(kind=dp),    dimension(ndirs), intent(in) :: spec_in
          integer(kind=int32_t),    intent(in) :: ndirs
          ! EXec code
          peakFrequency = sum(spectrum.f*spec_in**4*spectrum.df) &
                          /(sum(spec_in**4*spectrum.df)+eps)
    end function peakFrequency

!DIR$   ATTRIBUTES INLINE :: peakFrequencyDiscrete
    pure elemental real(kind=dp) function peakFrequencyDiscrete(spectrum,spec_in,ndirs)
!! Returns the peak frequency based on simple discrete maximum location of
  !! the spectrum array.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: peakFrequencyDiscrete
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_in:64
          real(kind=dp),    dimension(ndirs), intent(in) :: spec_in
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Exec code
          peakFrequencyDiscrete = &
                spectrum.f(maxloc(spec_in,dim=1))
    end function peakFrequencyDiscrete

!DIR$   ATTRIBUTES INLINE :: momentumFlux_yy
    pure elemental real(kind=dp) function momentumFlux_yy(spectrum)
 !! Returns total advective flux [kg/m^2/s^2] in y-direction of momentum in
  !! y-direction.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: momentumFlux_yy
          type(spectrum_type),  intent(in) :: spectrum

          ! Locals
          integer(kind=int32_t), automatic :: n,nfreqs
          ! Exec code ....
          nfreqs = size(spectrum.f)
          momentumFlux_yy = 0.0_dp
!DIR$   SIMD  REDUCTION(+:momentumFlux_yy)
          do n=1, nfreqs
               momentumFlux_yy = momentumFlux_yy            &
               + sum(spectrum.spec(n,:)*spectrum.dth*sin(spectrum.th)**2) &
               * spectrum.df(n)*spectrum.cg(n)/spectrum.cp(n)
          end do
          momentumFlux_yy = momentumFlux_yy*spectrum.water_density*spectrum.grav
    end function momentumFlux_yy

    pure elemental function wavenumber(f,depth,water_density,grav, &
                                        surface_tension)
          use mod_constants, only : pi2_const,eps
!! Solves the linear water wave dispersion relationship using a
 !! Newton-Raphson iteration loop.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: wavenumber
          real(kind=dp),    intent(in) :: f
          real(kind=dp),    intent(in), optional :: depth
          real(kind=dp),    intent(in), optional :: water_density
          real(kind=dp),    intent(in), optional :: grav
          real(kind=dp),    intent(in), optional :: surface_tension
          ! LOcals
          real(kind=dp), automatic :: dk,b,fnd,t
          real(kind=dp) :: wavenumber
          integer(kind=int32_t), automatic :: counter
          ! EXec code ....
          associate(k=>wavenumber)
                fnd = pi2_const*f*sqrt(depth/grav)
                k   = fnd*fnd
                b = surface_tension/(water_density*grav*depth**2)
                counter = 1
                dk = 2.0e-3_dp
                newton_raphson: do
                        t = tanh(k)
                        dk = -(fnd**2-k*t*(1.0_dp+b*k**2)) &
                           /(3.0_dp*b*k**2*t+t+k*(1.0_dp+b*k**2)*(1.0_dp-t**2))
                        k = k-dk
                        if(abs(dk) < eps .or. counter > 100) then
                                exit newton_raphson
                        end if
                        counter = counter + 1
               end do newton_raphson
               k = k/depth
          end associate
    end function wavenumber

!DIR$ ATTRIBUTES INLINE :: wavenumberMoment
    pure elemental real(kind=dp) function wavenumberMoment(spectrum,spec_in,nfreqs,ndirs,n)
!! Returns the spectral wavenumber moment of order n.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: wavenumberMoment
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_in:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(in) :: spec_in
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          integer(kind=int32_t),    intent(in) :: n
          ! Exec code ...
          wavenumberMoment = sum(spectrum.k**n*sum(spec_in,dim=2) &
                             * spectrum.dk
    end function wavenumberMoment

!DIR$   ATTRIBUTES INLINE :: significantWaveHeight
    pure elemental real(kind=dp) function significantWaveHeight(spectrum)
 !! Returns the significant wave height [m].
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: significantWaveHeight
          type(spectrum_type),      intent(in) :: spectrum
          ! Exec code ....
          significantWaveHeight =  4.0_dp*sqrt(frequencyMoment(spectrum,0))
    end function significantWaveHeight

!DIR$   ATTRIBUTES INLINE :: significantSurfaceOrbitalVelocity
    pure elemental real(kind=dp) function significantSurfaceOrbitalVelocity(spectrum) &
            result(uorb)
 !! Returns the significant surface orbital velocity [m/s].
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: significantSurfaceOrbitalVelocity
          GMS_use mod_constants, only : pi2_const
          type(spectrum_type),      intent(in) :: spectrum
          ! Exec code
          uorb = 2*sqrt(sum((pi2_const*spectrum.f)**2*sum(spectrum.spec,dim=2)*spectrum.df))
    end function significantSurfaceOrbitalVelocity

!DIR$   ATTRIBUTES INLINE :: meanPeriod
    pure elemental real(kind=p) function meanPeriod(spectrum)
 !! Returns the mean wave period [s].
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: meanPeriod
          use mod_constants, only : eps
          type(spectrum_type),      intent(in) :: spectrum
          ! Exec code ...
          meanPeriod = frequencyMoment(spectrum,0)/(frequencyMoment(spectrum,1)+eps)
    end function meanPeriod

!DIR$   ATTRIBUTES INLINE :: meanPeriodZeroCrossing
    pure elemental real(kind=dp) function meanPeriodZeroCrossing(spectrum)
 !! Returns the zero-crossing mean wave period [s]:
  !!
  !! Tm02 = \sqrt(m_0 / m_2)
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: meanPeriodZeroCrossing
          use GMS_mod_constants, only : eps
          type(spectrum_type),  intent(in) :: spectrum
          ! Exec code ...
          meanPeriodZeroCrossing = sqrt(frequencyMoment(spectrum,0) &
                                    /(frequencyMoment(2)+eps)
    end function meanPeriodZeroCrossing

    subroutine wavenumberSpectrum(spectrum,spec_out,nfreqs,ndirs)
          use GMS_mod_constants, only : v1_over_pi
!! Returns the wavenumber spectrum array of the spectrum instance.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: wavenumberSpectrum
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_out:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: spec_out
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t) :: ndir
          ! Exec code
          do concurrent(ndir=1:ndirs)
                spec_out(:,ndir) =  spectrum.spec(nfreq,:)*spectrum.cg(nfreq)*v1_over_pi
          end do
    end subroutine wavenumberSpectrum

    subroutine saturationSpectrum(spectrum,spec_out,wavenum_spec,nfreqs,ndirs)
!! Returns the saturation spectrum B(k) = F(k)k^4.
!DIR$   ATTRIBUTES CODE_ALIGN:32 :: saturationSpectrum
          type(spectrum_type),      intent(in) :: spectrum
!DIR$   ASSUME_ALIGNED spec_out:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(out) :: spec_out
!DIR$   ASSUME_ALIGNED wavenum_spec:64
          real(kind=dp),    dimension(nfreqs,ndirs), intent(in)  :: wavenum_spec
          integer(kind=int32_t),    intent(in) :: nfreqs
          integer(kind=int32_t),    intent(in) :: ndirs
          ! Locals
          integer(kind=int32_t) :: n
          ! Exec code...
          do concurrent(n=1:ndirs)
                spec_out(:,n) = wavenum_spec(:,n)*spectrum.k**4
          end do
    end subroutine saturationSpectrum

end module GMS_mod_spectrum

max_stars_repo_path: Wave model/GMS_mod_spectrum.f90

max_issues_repo_path: Wave model/GMS_mod_spectrum.f90

max_forks_repo_path: Wave model/GMS_mod_spectrum.f90


text:
*
*     -------------------------------------------------------------
*      C L E 0 S M
*     -------------------------------------------------------------
*
*                                                                  *
*     THIS PACKAGE DETERMINES THE VALUES OF FOLLOWING              *
*                                                 ---         ---  *
*                                                 I  Q   S  C   I  *
*     CLEBSCH - GORDAN COEFFICIENT:               I             I  *
*                                                 I  QM  0  CM  I  *
*                                                 ---         ---  *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vilnius,  Lithuania                          December 1993   *
*
      SUBROUTINE CLE0SM(Q,QM,S,C,CM,A)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      A=ZERO
      IIQ=TWO*Q+TENTH
      IIC=TWO*C+TENTH
      IIS=TWO*S+TENTH
      IF(ITTK(IIQ,IIC,IIS).EQ.0)RETURN
      IF(S.LT.EPS)GO TO 1
      CALL C1E0SM(Q,QM,C,CM,A)
      RETURN
    1 IF((Q+TENTH).LT.ABS(QM))RETURN
      IF((C+TENTH).LT.ABS(CM))RETURN
      IF(ABS(Q-C).GT.EPS)RETURN
      IF(ABS(QM-CM).GT.EPS)RETURN
      A=ONE
      RETURN
      END

ext: f


text:
module potential_mumps

!NOTES:
! - IF ONE REALLY WANTED TO CLEAN THIS UP IT MIGHT BE MORE EFFICIENT TO USE HARWELL-BOEING FORMAT
!   FOR MATRICES...

!SOME SUPERFLUOUS ARGUMENTS THAT ARE LEFT IN TO MAINTAIN UNIFORMITY ACROSS CALLS
!/home/zettergm/zettergmdata/GEMINI/numerical/potential/potential_mumps.f90:1291:0:
!warning: unused parameter ‘vminx1’ [-Wunused-parameter]
!   function
!elliptic2D_nonint_curv(srcterm,sig0,sigP,Vminx1,Vmaxx1,Vminx3,Vmaxx3,x,flagdirich,perflag,it)
! ^
!/home/zettergm/zettergmdata/GEMINI/numerical/potential/potential_mumps.f90:764:0:
!warning: unused parameter ‘vminx3’ [-Wunused-parameter]
!   function
!elliptic2D_pol_conv_curv_periodic2(srcterm,SigP,SigH,Cm,v2,v3,Vminx2,Vmaxx2,Vminx3,Vmaxx3,dt,x,Phi0,perflag,it)
! ^
!/home/zettergm/zettergmdata/GEMINI/numerical/potential/potential_mumps.f90:764:0:
!warning: unused parameter ‘vmaxx3’ [-Wunused-parameter]
!/home/zettergm/zettergmdata/GEMINI/numerical/potential/potential_mumps.f90:22:0:
!warning: unused parameter ‘vminx1’ [-Wunused-parameter]
!   function
!elliptic3D_curv(srcterm,sig0,sigP,sigH,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3,
!&
! ^

use, intrinsic:: iso_fortran_env, only: stderr=>error_unit, stdout=>output_unit

use phys_consts, only: wp, debug
use calculus, only: grad3D1, grad3D2, grad3D3
use meshobj, only: curvmesh
use interpolation, only: interp1
use PDEelliptic, only: elliptic3D_cart, elliptic3D_cart_periodic

implicit none (type, external)
private
public :: potential3D_fieldresolved_decimate, potential2D_polarization, potential2D_polarization_periodic, &
            potential2D_fieldresolved, potential3D_fieldresolved, potential3D_fieldresolved_truncate

integer, dimension(:), pointer, protected, save :: mumps_perm   !cached permutation, unclear whether save is necessary...

interface ! potential2d.f90
  module function potential2D_polarization(srcterm,SigP2,SigP3,SigH,Cm,v2,v3,Vminx2,Vmaxx2,Vminx3,Vmaxx3,dt,x,Phi0,perflag,it)
    real(wp), dimension(:,:), intent(in) :: srcterm,SigP2,SigP3,SigH,Cm,v2,v3
    !! ZZZ - THESE WILL NEED TO BE MODIFIED CONDUCTIVITIES, AND WE'LL NEED THREE OF THEM
    real(wp), dimension(:), intent(in) :: Vminx2,Vmaxx2
    real(wp), dimension(:), intent(in) :: Vminx3,Vmaxx3
    real(wp), intent(in) :: dt
    class(curvmesh), intent(in) :: x
    real(wp), dimension(:,:), intent(in) :: Phi0
    logical, intent(in) :: perflag
    integer, intent(in) :: it
    real(wp), dimension(size(SigP2,1),size(SigP2,2)) :: potential2D_polarization
  end function potential2D_polarization

  module function potential2D_polarization_periodic(srcterm,SigP,SigH,Cm,v2,v3,Vminx2,Vmaxx2,Vminx3,Vmaxx3,dt,x,Phi0,perflag,it)
    real(wp), dimension(:,:), intent(in) :: srcterm,SigP,SigH,Cm,v2,v3
    real(wp), dimension(:), intent(in) :: Vminx2,Vmaxx2
    real(wp), dimension(:), intent(in) :: Vminx3,Vmaxx3
    real(wp), intent(in) :: dt
    class(curvmesh), intent(in) :: x
    real(wp), dimension(:,:), intent(in) :: Phi0
    logical, intent(in) :: perflag
    integer, intent(in) :: it
    real(wp), dimension(size(SigP,1),size(SigP,2)) :: potential2D_polarization_periodic
  end function potential2D_polarization_periodic

  module function potential2D_fieldresolved(srcterm,sig0,sigP,Vminx1,Vmaxx1,Vminx3,Vmaxx3,x,flagdirich,perflag,it)
    real(wp), dimension(:,:,:), intent(in) :: srcterm,sig0,sigP   !arrays passed in will still have full rank 3
    real(wp), dimension(:,:), intent(in) :: Vminx1,Vmaxx1
    real(wp), dimension(:,:), intent(in) :: Vminx3,Vmaxx3
    class(curvmesh), intent(in) :: x
    integer, intent(in) :: flagdirich
    logical, intent(in) :: perflag
    integer, intent(in) :: it
    real(wp), dimension(size(sig0,1),size(sig0,2),size(sig0,3)) :: potential2D_fieldresolved
  end function potential2D_fieldresolved
end interface

contains
  function potential3D_fieldresolved_decimate(srcterm,sig0,sigP,sigH,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
                    x,flagdirich,perflag,it)
    !! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 3D USING MUMPS
    !! ASSUME THAT WE ARE RESOLVING THE POTENTIAL ALONG THE FIELD
    !! LINE.  THIS IS MOSTLY INEFFICIENT/UNWORKABLE FOR MORE THAN 1M
    !! GRID POINTS.

    real(wp), dimension(:,:,:), intent(in) :: srcterm,sig0,sigP,sigH
    real(wp), dimension(:,:), intent(in) :: Vminx1,Vmaxx1
    real(wp), dimension(:,:), intent(in) :: Vminx2,Vmaxx2
    real(wp), dimension(:,:), intent(in) :: Vminx3,Vmaxx3
    class(curvmesh), intent(in) :: x
    integer, intent(in) :: flagdirich
    logical, intent(in) :: perflag
    integer, intent(in) :: it

    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsigP2,gradsigP3
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsigH2,gradsigH3
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsig01
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: Ac,Bc,Cc,Dc,Ec,Fc

    integer :: lx1,lx2,lx3,ix1,ix2,ix3
    integer, parameter :: ldec=11
    real(wp), dimension(:), allocatable :: x1dec
    real(wp), dimension(:), allocatable :: dx1dec
    real(wp), dimension(:), allocatable :: x1idec
    real(wp), dimension(:), allocatable :: dx1idec
    real(wp), dimension(:,:,:), allocatable :: Acdec,Bcdec,Ccdec,Dcdec,Ecdec,Fcdec,srctermdec
    real(wp), dimension(:,:), allocatable :: Vminx2dec,Vmaxx2dec
    real(wp), dimension(:,:), allocatable :: Vminx3dec, Vmaxx3dec
    real(wp), dimension(:,:,:), allocatable :: Phidec

    real(wp), dimension(1:size(Vminx1,1),1:size(Vminx1,2)) :: Vminx1pot,Vmaxx1pot

    real(wp), dimension(size(srcterm,1),size(srcterm,2),size(srcterm,3)) :: potential3D_fieldresolved_decimate

    !SYSTEM SIZES
    lx1=x%lx1    !These will be full grid sizes if called from root (only acceptable thing)
    lx2=x%lx2all
    lx3=x%lx3all

    !COMPUTE AUXILIARY COEFFICIENTS TO PASS TO CART SOLVER
    if (debug) print *, 'Prepping coefficients for elliptic equation...'
    gradsig01=grad3D1(sig0,x,1,lx1,1,lx2,1,lx3)
    gradsigP2=grad3D2(sigP,x,1,lx1,1,lx2,1,lx3)
    gradsigP3=grad3D3(sigP,x,1,lx1,1,lx2,1,lx3)
    gradsigH2=grad3D2(sigH,x,1,lx1,1,lx2,1,lx3)
    gradsigH3=grad3D3(sigH,x,1,lx1,1,lx2,1,lx3)

    Ac=sigP
    Bc=sigP
    Cc=sig0
    Dc=gradsigP2+gradsigH3
    Ec=gradsigP3-gradsigH2
    Fc=gradsig01

    !DEFINE A DECIMATED MESH (THIS IS HARDCODED FOR NOW)
    if (debug) print*, 'Decimating parallel grid...'
    allocate(x1dec(-1:ldec+2),dx1dec(0:ldec+2),x1idec(1:ldec+1),dx1idec(1:ldec))
    x1dec(-1:lx1+2)=[x%x1(-1),x%x1(0),x%x1(1),81.8e3_wp,84.2e3_wp,87.5e3_wp,93.3e3_wp,106.0e3_wp,124.0e3_wp, &
                144.6e3_wp,206.7e3_wp,882.2e3_wp,x%x1(lx1),x%x1(lx1+1),x%x1(lx1+2)]
    !x1dec(-1:lx1+2)=[x%x1(-1),x%x1(0),x%x1(1),81.8e3_wp,84.2e3_wp,87.5e3_wp,93.3e3_wp,106.0e3_wp,124.0e3_wp, &
    !            144.6e3_wp,175e3_wp,206.7e3_wp,250e3_wp,400e3_wp,600e3_wp,882.2e3_wp,x%x1(lx1),x%x1(lx1+1),x%x1(lx1+2)]

    dx1dec(0:ldec+2)=x1dec(0:ldec+2)-x1dec(-1:ldec+1)
    x1idec(1:ldec+1)=0.5_wp*(x1dec(0:ldec)+x1dec(1:ldec+1))
    dx1idec(1:ldec)=x1idec(2:ldec+1)-x1idec(1:ldec)

    !INTERPOLATE COEFFICIENTS AND SOURCE TERM ONTO DECIMATED GRID
    if (debug) print*, 'Interpolating coefficients...'
    allocate(Acdec(1:ldec,1:lx2,1:lx3),Bcdec(1:ldec,1:lx2,1:lx3),Ccdec(1:ldec,1:lx2,1:lx3), &
             Dcdec(1:ldec,1:lx2,1:lx3),Ecdec(1:ldec,1:lx2,1:lx3),Fcdec(1:ldec,1:lx2,1:lx3), &
             srctermdec(1:ldec,1:lx2,1:lx3))
    do ix2=1,lx2
      do ix3=1,lx3
        Acdec(:,ix2,ix3)=interp1(x%x1(1:lx1),Ac(:,ix2,ix3),x1dec(1:ldec))
        Bcdec(:,ix2,ix3)=interp1(x%x1(1:lx1),Bc(:,ix2,ix3),x1dec(1:ldec))
        Ccdec(:,ix2,ix3)=interp1(x%x1(1:lx1),Cc(:,ix2,ix3),x1dec(1:ldec))
        Dcdec(:,ix2,ix3)=interp1(x%x1(1:lx1),Dc(:,ix2,ix3),x1dec(1:ldec))
        Ecdec(:,ix2,ix3)=interp1(x%x1(1:lx1),Ec(:,ix2,ix3),x1dec(1:ldec))
        Fcdec(:,ix2,ix3)=interp1(x%x1(1:lx1),Fc(:,ix2,ix3),x1dec(1:ldec))
        srctermdec(:,ix2,ix3)=interp1(x%x1(1:lx1),srcterm(:,ix2,ix3),x1dec(1:ldec))
      end do
    end do

    !INTERPOLATE BOUNDARY CONDITIONS ONTO DECIMATED GRID
    allocate(Vminx2dec(1:ldec,1:lx3),Vmaxx2dec(1:ldec,1:lx3))
    do ix3=1,lx3
      Vminx2dec(:,ix3)=interp1(x%x1(1:lx1),Vminx2(:,ix3),x1dec(1:ldec))
      Vmaxx2dec(:,ix3)=interp1(x%x1(1:lx1),Vmaxx2(:,ix3),x1dec(1:ldec))
    end do
    allocate(Vminx3dec(1:ldec,1:lx2),Vmaxx3dec(1:ldec,1:lx2))
    do ix2=1,lx2
      Vminx3dec(:,ix2)=interp1(x%x1(1:lx1),Vminx3(:,ix2),x1dec(1:ldec))
      Vmaxx3dec(:,ix2)=interp1(x%x1(1:lx1),Vmaxx3(:,ix2),x1dec(1:ldec))
    end do

    !FOR WHATEVER REASON THE EDGE VALUES GET MESSED UP BY INTERP1
    Acdec(1,:,:)=Ac(1,:,:)
    Acdec(ldec,:,:)=Ac(lx1,:,:)
    Bcdec(1,:,:)=Bc(1,:,:)
    Bcdec(ldec,:,:)=Bc(lx1,:,:)
    Ccdec(1,:,:)=Cc(1,:,:)
    Ccdec(ldec,:,:)=Cc(lx1,:,:)
    Dcdec(1,:,:)=Dc(1,:,:)
    Dcdec(ldec,:,:)=Dc(lx1,:,:)
    Ecdec(1,:,:)=Ec(1,:,:)
    Ecdec(ldec,:,:)=Ec(lx1,:,:)
    Fcdec(1,:,:)=Fc(1,:,:)
    Fcdec(ldec,:,:)=Fc(lx1,:,:)
    Vminx2dec(1,:)=Vminx2(1,:)
    Vminx2dec(ldec,:)=Vminx2(lx1,:)
    Vmaxx2dec(1,:)=Vmaxx2(1,:)
    Vmaxx2dec(ldec,:)=Vmaxx2(lx1,:)
    Vminx3dec(1,:)=Vminx3(1,:)
    Vminx3dec(ldec,:)=Vminx3(lx1,:)
    Vmaxx3dec(1,:)=Vmaxx3(1,:)
    Vmaxx3dec(ldec,:)=Vmaxx3(lx1,:)
    srctermdec(1,:,:)=srcterm(1,:,:)
    srctermdec(ldec,:,:)=srcterm(lx1,:,:)

    !
    !print*, minval(Acdec),maxval(Acdec)
    !print*, minval(Ac),maxval(Ac)
    !print*, minval(Bcdec),maxval(Bcdec)
    !print*, minval(Bc),maxval(Bc)
    !print*, minval(Ccdec),maxval(Ccdec)
    !print*, minval(Cc),maxval(Cc)
    !print*, minval(Dcdec),maxval(Dcdec)
    !print*, minval(Dc),maxval(Dc)
    !print*, minval(Ecdec),maxval(Ecdec)
    !print*, minval(Ec),maxval(Ec)
    !print*, minval(Fcdec),maxval(Fcdec)
    !print*, minval(Fc),maxval(Fc)
    !print*, minval(srctermdec),maxval(srctermdec)
    !print*, minval(srcterm),maxval(srcterm)
    !print*, minval(Vminx2dec),maxval(Vminx2dec)
    !print*, minval(Vmaxx2dec),maxval(Vmaxx2dec)
    !print*, minval(Vminx3dec),maxval(Vminx3dec)
    !print*, minval(Vmaxx3dec),maxval(Vmaxx3dec)
    !print*, minval(dx1dec),maxval(dx1dec)
    !print*, minval(dx1idec),maxval(dx1idec)
    !print*, x1dec(-1:ldec+2)
    !print*, dx1dec(0:ldec+2)
    !

    !ADJUST THE BOUNDARY CONDITION TO POTENTIAL DERIVATIVE INSTEAD OF CURRENT DENSITY
    Vminx1pot=-1*Vminx1/sig0(1,:,:)
    Vmaxx1pot=-1*Vmaxx1/sig0(lx1,:,:)

    !CALL CARTESIAN SOLVER ON THE DECIMATED GRID
    if (debug) print*, 'Calling solve on decimated grid...'
    allocate(Phidec(1:ldec,1:lx2,1:lx3))
    Phidec=elliptic3D_cart(srctermdec,Acdec,Bcdec,Ccdec,Dcdec,Ecdec,Fcdec,Vminx1pot,Vmaxx1pot, &
                    Vminx2dec,Vmaxx2dec,Vminx3dec,Vmaxx3dec, &
                    dx1dec,dx1idec,x%dx2all,x%dx2iall,x%dx3all,x%dx3iall,flagdirich,perflag,it)

    !INTERPOLATE BACK UP TO MAIN GRID
    if (debug) print*, 'Upsampling potential...'
    do ix2=1,lx2
      do ix3=1,lx3
        potential3D_fieldresolved_decimate(:,ix2,ix3)=interp1(x1dec(1:ldec),Phidec(:,ix2,ix3),x%x1(1:lx1))
      end do
    end do

    !AGAIN NEED TO FIX THE EDGES...
    potential3D_fieldresolved_decimate(1,:,:)=Phidec(1,:,:)
    potential3D_fieldresolved_decimate(lx1,:,:)=Phidec(ldec,:,:)

    ! open(newunit=u, form='unformatted', access='stream',file='Phidec.raw8',status='replace', action='write')
    ! write(u) potential3D_fieldresolved_decimate,Phidec,Ac,Acdec,Bc,Bcdec,Cc,Ccdec,Dc,Dcdec,Ec,Ecdec,Fc,Fcdec,srcterm,srctermdec
    ! write(u) Vminx1pot,Vmaxx1pot,Vminx2dec,Vmaxx2dec,Vminx3dec,Vmaxx3dec
    ! close(u)

    !CLEAN UP THE ALLOCATED ARRAYS
    deallocate(srctermdec,Acdec,Bcdec,Ccdec,Dcdec,Ecdec,Fcdec,dx1dec,x1dec,x1idec,dx1idec,Vminx2dec,Vmaxx2dec,Vminx3dec,Vmaxx3dec)
    deallocate(Phidec)
  end function potential3D_fieldresolved_decimate

  function potential3D_fieldresolved(srcterm,sig0,sigP,sigH,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
                    x,flagdirich,perflag,it)
    !! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 3D USING MUMPS
    !! ASSUME THAT WE ARE RESOLVING THE POTENTIAL ALONG THE FIELD
    !! LINE.  THIS IS MOSTLY INEFFICIENT/UNWORKABLE FOR MORE THAN 1M
    !! GRID POINTS.
    real(wp), dimension(:,:,:), intent(in) :: srcterm,sig0,sigP,sigH
    real(wp), dimension(:,:), intent(in) :: Vminx1,Vmaxx1
    real(wp), dimension(:,:), intent(in) :: Vminx2,Vmaxx2
    real(wp), dimension(:,:), intent(in) :: Vminx3,Vmaxx3
    class(curvmesh), intent(in) :: x
    integer, intent(in) :: flagdirich
    logical, intent(in) :: perflag
    integer, intent(in) :: it
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsigP2,gradsigP3
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsigH2,gradsigH3
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsig01
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: Ac,Bc,Cc,Dc,Ec,Fc
    integer :: lx1,lx2,lx3,ix1,ix2,ix3
    integer, parameter :: ldec=11
    real(wp), dimension(:), allocatable :: x1dec
    real(wp), dimension(:), allocatable :: dx1dec
    real(wp), dimension(:), allocatable :: x1idec
    real(wp), dimension(:), allocatable :: dx1idec
    real(wp), dimension(:,:,:), allocatable :: Acdec,Bcdec,Ccdec,Dcdec,Ecdec,Fcdec,srctermdec
    real(wp), dimension(:,:), allocatable :: Vminx2dec,Vmaxx2dec
    real(wp), dimension(:,:), allocatable :: Vminx3dec, Vmaxx3dec
    real(wp), dimension(:,:,:), allocatable :: Phidec
    real(wp), dimension(1:size(Vminx1,1),1:size(Vminx1,2)) :: Vminx1pot,Vmaxx1pot
    real(wp), dimension(size(srcterm,1),size(srcterm,2),size(srcterm,3)) :: potential3D_fieldresolved

    !SYSTEM SIZES
    lx1=x%lx1    !These will be full grid sizes if called from root (only acceptable thing)
    lx2=x%lx2all
    lx3=x%lx3all

    !COMPUTE AUXILIARY COEFFICIENTS TO PASS TO CART SOLVER
    if (debug) print *, 'Prepping coefficients for elliptic equation...'
    gradsig01=grad3D1(sig0,x,1,lx1,1,lx2,1,lx3)
    gradsigP2=grad3D2(sigP,x,1,lx1,1,lx2,1,lx3)
    gradsigP3=grad3D3(sigP,x,1,lx1,1,lx2,1,lx3)
    gradsigH2=grad3D2(sigH,x,1,lx1,1,lx2,1,lx3)
    gradsigH3=grad3D3(sigH,x,1,lx1,1,lx2,1,lx3)

    ! Need to correct the derivatives if periodic in x3 chosen; FIXME: assume cartesian for now
    !!  FIXME: ideally this needs to be a routine from within the gradent submodule of calculus
    if (x%flagper) then
      if (debug) print*, 'Adjusting conductivity derivatives to account for periodic grid...'
      gradsigP3(:,:,1)=(sigP(:,:,2)-sigP(:,:,lx3))/(x%dx3all(2)+x%dx3all(1))
      gradsigH3(:,:,1)=(sigH(:,:,2)-sigH(:,:,lx3))/(x%dx3all(2)+x%dx3all(1))
      gradsigP3(:,:,lx3)=(sigP(:,:,1)-sigP(:,:,lx3-1))/(x%dx3all(1)+x%dx3all(lx3))
      gradsigH3(:,:,lx3)=(sigH(:,:,1)-sigH(:,:,lx3-1))/(x%dx3all(1)+x%dx3all(lx3))
    end if

    ! coefficients for 3D solve
    Ac=sigP
    Bc=sigP
    Cc=sig0
    Dc=gradsigP2+gradsigH3
    Ec=gradsigP3-gradsigH2
    Fc=gradsig01

    !ADJUST THE BOUNDARY CONDITION TO POTENTIAL DERIVATIVE INSTEAD OF CURRENT DENSITY
    Vminx1pot=-1*Vminx1/sig0(1,:,:)
    Vmaxx1pot=-1*Vmaxx1/sig0(lx1,:,:)

    ! call solver on the full grid
    if (x%flagper) then
      if (debug) print*, 'Calling 3D solve on full grid, periodic'
      potential3D_fieldresolved=elliptic3D_cart_periodic(srcterm,Ac,Bc,Cc,Dc,Ec,Fc,Vminx1pot,Vmaxx1pot, &
                      Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
                      x%dx1,x%dx1i,x%dx2all,x%dx2iall,x%dx3all,x%dx3iall,flagdirich,perflag,it)
    else
      if (debug) print*, 'Calling 3D solve on full grid, non-periodic'
      potential3D_fieldresolved=elliptic3D_cart(srcterm,Ac,Bc,Cc,Dc,Ec,Fc,Vminx1pot,Vmaxx1pot, &
                      Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
                      x%dx1,x%dx1i,x%dx2all,x%dx2iall,x%dx3all,x%dx3iall,flagdirich,perflag,it)
    end if
  end function potential3D_fieldresolved

  function potential3D_fieldresolved_truncate(srcterm,sig0,sigP,sigH,Vminx1,Vmaxx1,Vminx2,Vmaxx2,Vminx3,Vmaxx3, &
                    x,flagdirich,perflag,it)
    !! SOLVE IONOSPHERIC POTENTIAL EQUATION IN 3D USING MUMPS
    !! ASSUME THAT WE ARE RESOLVING THE POTENTIAL ALONG THE FIELD
    !! LINE.  THIS IS MOSTLY INEFFICIENT/UNWORKABLE FOR MORE THAN 1M
    !! GRID POINTS.
    real(wp), dimension(:,:,:), intent(in) :: srcterm,sig0,sigP,sigH
    real(wp), dimension(:,:), intent(in) :: Vminx1,Vmaxx1
    real(wp), dimension(:,:), intent(in) :: Vminx2,Vmaxx2
    real(wp), dimension(:,:), intent(in) :: Vminx3,Vmaxx3
    class(curvmesh), intent(in) :: x
    integer, intent(in) :: flagdirich
    logical, intent(in) :: perflag
    integer, intent(in) :: it
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsigP2,gradsigP3
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsigH2,gradsigH3
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: gradsig01
    real(wp), dimension(1:size(srcterm,1),1:size(srcterm,2),1:size(srcterm,3)) :: Ac,Bc,Cc,Dc,Ec,Fc
    integer :: lx1,lx2,lx3,ix1,ix2,ix3
    real(wp), dimension(:), allocatable :: x1trunc
    real(wp), dimension(:), allocatable :: dx1trunc
    real(wp), dimension(:), allocatable :: x1itrunc
    real(wp), dimension(:), allocatable :: dx1itrunc
    real(wp), dimension(:,:,:), allocatable :: Actrunc,Bctrunc,Cctrunc,Dctrunc,Ectrunc,Fctrunc,srctermtrunc
    real(wp), dimension(:,:), allocatable :: Vminx2trunc,Vmaxx2trunc
    real(wp), dimension(:,:), allocatable :: Vminx3trunc, Vmaxx3trunc
    real(wp), dimension(:,:,:), allocatable :: Phitrunc
    real(wp), dimension(1:size(Vminx1,1),1:size(Vminx1,2)) :: Vminx1pot,Vmaxx1pot
    real(wp), dimension(size(srcterm,1),size(srcterm,2),size(srcterm,3)) :: potential3D_fieldresolved_truncate
    integer :: lx1trunc
    real(wp) :: alttrunc,alt

    !SYSTEM SIZES
    lx1=x%lx1    !These will be full grid sizes if called from root (only acceptable thing)
    lx2=x%lx2all
    lx3=x%lx3all

    !COMPUTE AUXILIARY COEFFICIENTS TO PASS TO CART SOLVER
    if (debug) print *, 'Prepping coefficients for elliptic equation...'
    gradsig01=grad3D1(sig0,x,1,lx1,1,lx2,1,lx3)
    gradsigP2=grad3D2(sigP,x,1,lx1,1,lx2,1,lx3)
    gradsigP3=grad3D3(sigP,x,1,lx1,1,lx2,1,lx3)
    gradsigH2=grad3D2(sigH,x,1,lx1,1,lx2,1,lx3)
    gradsigH3=grad3D3(sigH,x,1,lx1,1,lx2,1,lx3)

    ! Need to correct the derivatives if periodic in x3 chosen; FIXME: assume cartesian for now
    !!  FIXME: ideally this needs to be a routine from within the gradent submodule of calculus
    if (x%flagper) then
      if (debug) print*, 'Adjusting conductivity derivatives to account for periodic grid...'
      gradsigP3(:,:,1)=(sigP(:,:,2)-sigP(:,:,lx3))/(x%dx3all(2)+x%dx3all(1))
      gradsigH3(:,:,1)=(sigH(:,:,2)-sigH(:,:,lx3))/(x%dx3all(2)+x%dx3all(1))
      gradsigP3(:,:,lx3)=(sigP(:,:,1)-sigP(:,:,lx3-1))/(x%dx3all(1)+x%dx3all(lx3))
      gradsigH3(:,:,lx3)=(sigH(:,:,1)-sigH(:,:,lx3-1))/(x%dx3all(1)+x%dx3all(lx3))
    end if

    ! coefficients for 3D solve
    Ac=sigP
    Bc=sigP
    Cc=sig0
    Dc=gradsigP2+gradsigH3
    Ec=gradsigP3-gradsigH2
    Fc=gradsig01

    ! Now truncate the problem at some appropriate altitude
    ix1=1
    alttrunc=300e3
    do while(x%alt(ix1,1,1)<alttrunc)
      ix1=ix1+1
    end do
    if (ix1>=1) then
      lx1trunc=ix1-1
      print*, 'Truncating field aligned solve at index:  ',lx1trunc
    else
      print*, 'Unable to truncate field-resolved solve at sensible altitude...',ix1,alttrunc
      error stop
    end if

    !ADJUST THE BOUNDARY CONDITION TO POTENTIAL DERIVATIVE INSTEAD OF CURRENT DENSITY
    Vminx1pot=-1*Vminx1/sig0(1,:,:)
    Vmaxx1pot=-1*Vmaxx1/sig0(lx1trunc,:,:)

    ! allocate and assign truncated variables
    allocate(Actrunc(1:lx1trunc,1:lx2,1:lx3))
    allocate(Bctrunc,Cctrunc,Dctrunc,Ectrunc,Fctrunc,srctermtrunc, mold=Actrunc)
    allocate(Vminx2trunc(1:lx1trunc,1:lx3))
    allocate(Vmaxx2trunc, mold=Vminx2trunc)
    allocate(Vminx3trunc(1:lx1trunc,1:lx2))
    allocate(Vmaxx3trunc, mold=Vminx3trunc)
    allocate(dx1trunc(0:lx1trunc+2))
    allocate(dx1itrunc(1:lx1trunc+1))
    Actrunc=Ac(1:lx1trunc,1:lx2,1:lx3)
    Bctrunc=Bc(1:lx1trunc,1:lx2,1:lx3)
    Cctrunc=Cc(1:lx1trunc,1:lx2,1:lx3)
    Dctrunc=Dc(1:lx1trunc,1:lx2,1:lx3)
    Ectrunc=Ec(1:lx1trunc,1:lx2,1:lx3)
    Fctrunc=Fc(1:lx1trunc,1:lx2,1:lx3)
    srctermtrunc=srcterm(1:lx1trunc,1:lx2,1:lx3)
    Vminx2trunc=Vminx2(1:lx1trunc,1:lx3)
    Vmaxx2trunc=Vmaxx2(1:lx1trunc,1:lx3)
    Vminx3trunc=Vminx3(1:lx1trunc,1:lx2)
    Vmaxx3trunc=Vmaxx3(1:lx1trunc,1:lx2)
    dx1trunc=x%dx1(0:lx1trunc+2)
    dx1itrunc=x%dx1i(1:lx1trunc+1)

    ! call the solver for the truncated grid
    allocate(Phitrunc, mold=Actrunc)
    if (x%flagper) then
      if (debug) print*, 'Calling 3D solve on full grid, periodic'
      Phitrunc=elliptic3D_cart_periodic(srctermtrunc,Actrunc,Bctrunc,Cctrunc,Dctrunc,Ectrunc,Fctrunc,Vminx1pot,Vmaxx1pot, &
                      Vminx2trunc,Vmaxx2trunc,Vminx3trunc,Vmaxx3trunc, &
                      dx1trunc,dx1itrunc,x%dx2all,x%dx2iall,x%dx3all,x%dx3iall,flagdirich,perflag,it)
    else
      if (debug) print*, 'Calling 3D solve on full grid, non-periodic'
      Phitrunc=elliptic3D_cart(srctermtrunc,Actrunc,Bctrunc,Cctrunc,Dctrunc,Ectrunc,Fctrunc,Vminx1pot,Vmaxx1pot, &
                      Vminx2trunc,Vmaxx2trunc,Vminx3trunc,Vmaxx3trunc, &
                      dx1trunc,dx1itrunc,x%dx2all,x%dx2iall,x%dx3all,x%dx3iall,flagdirich,perflag,it)
    end if
    potential3D_fieldresolved_truncate(1:lx1trunc,1:lx2,1:lx3)=Phitrunc(1:lx1trunc,1:lx2,1:lx3)
    do ix1=lx1trunc+1,lx1
      potential3D_fieldresolved_truncate(ix1,1:lx2,1:lx3)=Phitrunc(lx1trunc,1:lx2,1:lx3)
    end do

    ! deallocate temporary arrays
    deallocate(Actrunc,Bctrunc,Cctrunc,Dctrunc,Ectrunc,Fctrunc,srctermtrunc,Vminx2trunc,Vmaxx2trunc, &
               Vminx3trunc,Vmaxx3trunc,dx1trunc,dx1itrunc,Phitrunc)
  end function potential3D_fieldresolved_truncate
end module potential_mumps

max_forks_repo_name: mattzett/GEMINI


text:
! { dg-do compile }
! Tests the fix for PR25951, a regression caused by the assumed
! size patch.
! Test case provided by Mark Hesselink  <mhesseli@caltech.edu>
PROGRAM loc_1
  integer i(10)
  call f (i)
CONTAINS
   SUBROUTINE f (x)
      INTEGER, DIMENSION(*)   :: x
      INTEGER                 :: address
! The next line would cause:
! Error: The upper bound in the last dimension must appear in the
! reference to the assumed size array 'x' at (1)
      address=LOC(x)
   END SUBROUTINE f
END PROGRAM loc_1

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program t
! test for nested do-loops
  integer::i,j,k,l
  do 77 i=1,2
    do 88 j=1,2
      do 88 k=1,2
        do 88 l=1,2
          print *,1
88        continue
77        continue
        end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


text:
! ur8_def.for - Ulysses REAL*8 time format constants for WIND/WAVES processing

	real*8		ur8_day, ur8_hr, ur8_min, ur8_sec, ur8_msec
	parameter	(ur8_day=1.0,
	1		ur8_hr=ur8_day/24.0,
	1		ur8_min=ur8_hr/60.0,
	1		ur8_sec=ur8_min/60.0,
	1		ur8_msec=ur8_sec/1000.0)

	real*8		ur8_1x_mf_diff
	parameter	(ur8_1x_mf_diff=ur8_sec*92.0/250.0)
	real*8		ur8_2x_mf_diff
	parameter	(ur8_2x_mf_diff=ur8_sec*46.0/250.0)

ext: for


text:
! { dg-do run }
! { dg-options "-g" }
  integer :: a(10), b(12)
  call sub (a, 10)
  call sub (b, 12)
  write (*,*) a, b
end

subroutine sub (a, n)
  integer :: a(n), n
  do i = 1, n
    a(i) = i
  end do
  write (*,*) a	! { dg-final { gdb-test 14 "a(10)" "10" } }
end subroutine


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  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 registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine actrep( noalg   , noprot   , namprot, nampact, nopralg,
     +                   nampralg, constants)

!     Deltares Software Centre

!>/File
!>      replace active proto processes with actual processes

      use timers         !< performance timers
      use dlwq_data      !< data definitions
      use processet      !< use processet definitions
      implicit none

      ! decalaration of arguments

      integer  ( 4)    , intent(in   ) :: noalg             !< number of algae types
      integer  ( 4)    , intent(in   ) :: noprot            !< number of proto processes ono to one
      integer  ( 4)    , intent(in   ) :: nopralg           !< number of proto processes per algae type
      character(len=*) , intent(in   ) :: namprot(noprot)   !< name proto processes
      character(len=*) , intent(in   ) :: nampact(noprot)   !< name actual processes
      character(len=*) , intent(in   ) :: nampralg(noprot)  !< name proto processes per algae type
      type(t_dlwq_item), intent(inout) :: constants         !< delwaq constants list

      ! local decalaration

      character*10                     :: name1
      character*20                     :: name20
      character*80                     :: line
      integer                          :: nocons
      integer                          :: nocon2
      integer                          :: ico
      integer                          :: ipro
      integer                          :: ilen
      integer                          :: ialg
      integer                          :: ierr2
      integer(4)                       :: ithndl = 0

      if (timon) call timstrt( "actrep", ithndl )

      ! replace active proto processes with actual processes

      nocons = constants%no_item
      nocon2 = constants%no_item
      do ico = 1 , nocons
         call dhucas(constants%name(ico),name20,20)
         if ( name20(1:6) .eq. 'ACTIVE' ) then

            ! the one to one processes

            do ipro = 1 , noprot
               if ( namprot(ipro) .eq. name20(8:17) ) then
                  constants%name(ico)(8:) = nampact(ipro)
               endif
            enddo

            ! the processes that need to be expanded for every algal type

            do ipro = 1 , nopralg
               if ( nampralg(ipro) .eq. name20(8:17) ) then
                  name1 = nampralg(ipro)
                  call dhslen(name1,ilen)
                  do ialg = 1 , noalg
                     write(name1(ilen+1:),'(i2.2)') ialg
                     if ( ialg .eq. 1 ) then
                        constants%name(ico)(8:) = name1
                     else
                        nocon2 = nocon2 + 1
                        ierr2 = dlwq_resize(constants,nocon2)
                        if ( ierr2 .gt. 0 ) then
                           write(line,'(a,i10)') ' ERROR: actrep resize error constants size:',nocon2
                           call monsys(line,1)
                           call srstop(1)
                        endif
                        constants%no_item          = nocon2
                        constants%name(nocon2)     = constants%name(ico)
                        constants%name(nocon2)(8:) = name1
                     endif
                  enddo
               endif
            enddo
         endif
      enddo

      if (timon) call timstop( ithndl )
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd