text:
C$Procedure    EXTMSI  ( print message with integer and do error exit )

      SUBROUTINE EXTMSI ( MESSGE, MARKER, ICODE )

C$ Abstract
C
C     This subroutine prints an error message with an integer
C     and then exits.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     ERROR
C
C$ Declarations

      IMPLICIT NONE

      INTEGER               ICODE

      CHARACTER*(*)         MESSGE
      CHARACTER*(*)         MARKER

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     MESSGE     I   Character string with message and an imbedded
C                    marker.
C     MARKER     I   Character string marker to be replaced.
C     ICODE      I   Integer that will be inserted into the message.
C
C$ Detailed_Input
C
C     MESSGE     is a character string with a message, and a marker
C                where the integer ICODE will be inserted.
C
C     MARKER     is a character string to be replaced by an integer.
C
C     ICODE      is an integer that will be inserted into the message
C                at the location of the marker.
C
C$ Detailed_Output
C
C     None.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     Error free.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This is a common task in MKPLAT.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     J.A. Bytof      (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.0, 08-OCT-2009 (NJB)
C
C        Re-ordered header sections.
C
C-    SPICELIB Version 1.0.0, 10-APR-1997 (JAB)
C
C-&

C$ Index_Entries
C
C     Error exit message with integer.
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN

C
C     Local variables.
C

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'EXTMSI' )
      END IF

      CALL SETMSG ( MESSGE              )
      CALL ERRINT ( MARKER, ICODE       )
      CALL SIGERR ( 'SPICE(ERROREXIT)'  )

      CALL CHKOUT ( 'EXTMSI' )

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!===============================================================================
!> Optimized restrict operations, intended for guard zones (not conservative)
!>
!>    1   2   3   4   5   6   7   8   9   0  11    17  18  19  20  21  22
!>    +---+---+-|-o---o---o---o---o---o---0---0-| --0---o---o---o-|-+---+---+
!>  -o-o-o-o-o-o|+-+-+                          |            +-+-+|o-o-o-o-o-o
!>   4 5 6 7 8 9 0 1 2                                       1 2 3 4 5 6 7 8 9
!>         o-o-o|o-o-o-o-o-o-o-o-o-o-o-o-o-o-o-o|o-o-o
!>         1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2
!>
!> When boundaries are lined up for 16x16x16 patches, restrict operations should
!> involve interpolations over index 14-15, 16-17, and 18-19 on the high source
!> side, and interpolations over index 4-5, 6-7, and 8-9 on the low source side,
!> while of coure 1-3 and 20-22 on the low and high target side.  When the
!> source overlaps with the target it covers target index 4-11 left and 12-19
!> right, with its own range 4-19.
!===============================================================================
MODULE guard_zones_mod
  USE io_mod
  USE io_unit_mod
  USE trace_mod
  USE patch_mod
  USE mesh_mod
  USE omp_timer_mod
  implicit none
  private
  type, public:: guard_zones_t
    procedure(restrict_simple), pointer:: method => null()
  contains
    procedure:: init
    procedure:: restrict_simple
    procedure:: different_opt
    procedure:: restrict_vect
    procedure, nopass:: info
  end type
  real(8):: update_t(3)=0d0
  integer:: verbose=0
  logical:: detailed_timer=.false.
  type(guard_zones_t), public:: guard_zones
CONTAINS

!===============================================================================
!> Choose method and verbosity
!===============================================================================
SUBROUTINE init (self)
  class(guard_zones_t):: self
  character(len=6):: method='opt'
  namelist /guard_zone_params/ method, detailed_timer, verbose
  integer:: iostat
  logical, save:: first_time=.true.
  !-----------------------------------------------------------------------------
  call trace%begin ('guard_zones_t%init')
  if (first_time) then
    !$omp critical (init_cr)
    if (first_time) then
      rewind (io%input)
      read (io%input, guard_zone_params, iostat=iostat)
      write (io%output, guard_zone_params)
      first_time = .false.
    end if
    !$omp end critical (init_cr)
  end if
  select case (trim(method))
  case ('simple')
    self%method => restrict_simple
  case ('opt')
    self%method => different_opt
  case ('vect')
    self%method => restrict_vect
  case default
    call io%abort ('guard zone method unknown')
  end select
  call trace%end ()
END SUBROUTINE init

!===============================================================================
!> Compute interpolated guard zone values by setting up interface to the
!> slow but correct patch_t%interpolate() point-by-point interpolation.
!===============================================================================
SUBROUTINE restrict_simple (self, target, source, jt, pt, same)
  class(guard_zones_t):: self
  class(patch_t), pointer:: target, source
  integer:: jt(2)
  real(4):: pt(2)
  logical, optional:: same
  !.............................................................................
  real, dimension(:,:,:), pointer:: tmem, smem, td, sd
  class(mesh_t), pointer:: mesh
  integer:: i1, i2, i3, iv, jv, tl(3), tu(3), sl(3), su(3)
  real(8):: tlc(3), tuc(3), slc(3), suc(3), dist(3)
  real(4):: pl(3), pu(3)
  real(8):: used
  integer, save:: itimer=0
  !-----------------------------------------------------------------------------
  ! Interpolate in time and space, point-by-point
  !-----------------------------------------------------------------------------
  call trace%begin ('guard_zones_t%restrict_simple', itimer=itimer)
  if (detailed_timer) &
    used = wallclock()
  call source%time_interval (target%time, jt, pt)
  do iv=1,target%nv
    jv = 0
    !---------------------------------------------------------------------------
    ! Target lower and upper corner in real space
    !---------------------------------------------------------------------------
    tlc = target%myposition (target%mesh%lb, jv)
    tuc = target%myposition (target%mesh%ub, jv)
    !---------------------------------------------------------------------------
    ! Map to source index space and limit to source interior
    !---------------------------------------------------------------------------
    sl = source%index_only (tlc, jv, nowrap=.true.)
    su = source%index_only (tuc, jv, nowrap=.true.)
    sl = min(max(sl,source%mesh%li),source%mesh%ui-1)
    su = min(max(su,source%mesh%li),source%mesh%ui-1)
    !---------------------------------------------------------------------------
    ! Map source points to target legal index space
    !---------------------------------------------------------------------------
    slc = source%myposition(sl, jv)
    suc = source%myposition(su, jv)
    tl = target%index_only (slc, jv, nowrap=.true., roundup=.true.)
    tu = target%index_only (suc, jv, nowrap=.true., roundup=.true.)
    tl = min(max(tl,target%mesh%lb),target%mesh%ub)
    tu = min(max(tu,target%mesh%lb),target%mesh%ub)
    if (verbose > 0) &
      write(io_unit%log,'(a,i3,2(2x,3i4),2(2x,3f7.3))') &
        'guard_zones_t%restrict_simple1: iv, tl, tu =', iv, tl, tu
    !call dnload_range (target, source, iv, tl, tu)
    !if (verbose > 0) &
    !  write(io_unit%log,'(a,i3,2(2x,3i4),2(2x,3f7.3))') &
    !    'guard_zones_t%restrict_simple2: iv, tl, tu =', iv, tl, tu
    do i3=tl(3),tu(3); do i2=tl(2),tu(2); do i1=tl(1),tu(1)
      call interpolate (source, target, [i1,i2,i3], iv, jt, pt)
    end do; end do; end do
  end do
  if (detailed_timer) then
    used = wallclock()-used
    write (stdout,'(a,6p,f6.3," mus/cell")') &
      "guard_zone_t%restrict_simple:", used/product(tu-tl+1)
  end if
  call trace%end (itimer)
END SUBROUTINE restrict_simple

!===============================================================================
!> Compute interpolated guard zone values by setting up interface to the
!> slow but correct patch_t%interpolate() point-by-point interpolation.
!===============================================================================
SUBROUTINE different_opt (self, target, source, jt, pt, same)
  class(guard_zones_t):: self
  class(patch_t), pointer:: target, source
  integer:: jt(2)
  real(4):: pt(2)
  logical, optional:: same
  !.............................................................................
  real, parameter:: eps=1e-4
  real, dimension(:,:,:), pointer:: tmem, smem, td, sd
  class(mesh_t), pointer:: mesh
  integer:: i, i1, i2, i3, iv, jv, tl(3), tu(3), sl(3), su(3), cells
  real:: pi(16)
  real(8):: tlc, tuc, slc, suc, dist, box, tc
  class(mesh_t), pointer:: tm, sm
  real(8):: start
  integer, save:: itimer=0
  !-----------------------------------------------------------------------------
  ! Interpolate in time and space, point-by-point
  !-----------------------------------------------------------------------------
  call trace%begin ('guard_zones_t%different_opt', itimer=itimer)
  if (detailed_timer) then
    start = wallclock()
  end if
  call source%time_interval (target%time, jt, pt)
  do iv=1,target%nv
    jv = 0
    do i=1,3
      !-------------------------------------------------------------------------
      ! Source distance from target center
      !-------------------------------------------------------------------------
      tm => target%mesh(i)
      sm => source%mesh(i)
      box = tm%b
      dist = sm%p - tm%p
      if (target%periodic(i)) &
        dist = modulo(dist+0.5d0*box,box) - 0.5d0*box
      !-------------------------------------------------------------------------
      ! Target corner distances from center
      !-------------------------------------------------------------------------
      tlc = tm%r(tm%lb)
      tuc = tm%r(tm%ub)
      !-------------------------------------------------------------------------
      ! Source lower and upper li and ui points, relative to target center
      !-------------------------------------------------------------------------
      slc = sm%r(sm%li) + dist
      suc = sm%r(sm%ui) + dist
      !-------------------------------------------------------------------------
      ! Lower and upper index distance from target center
      !-------------------------------------------------------------------------
      tlc = max(tlc,slc)
      tuc = min(tuc,suc)
      !-------------------------------------------------------------------------
      ! Target indices to lower and upper guard zone limits
      !-------------------------------------------------------------------------
      tc = 1d0/tm%d
      tl(i) = ceiling(tlc*tc + tm%o - eps)
      tu(i) = ceiling(tuc*tc + tm%o - eps)
      if (verbose > 0) &
        write(io_unit%log,'(a,i3,2x,5f8.3,2x,2i4)') &
        'guard_zones_t%different_opt: i,dist,tlc,tuc,slc,suc,tl,tu =', &
        i, dist, tlc, tuc, slc, suc, tl(i), tu(i)
    end do
    !call dnload_range (target, source, iv, tl, tu)
    !if (verbose > 0) &
    !  write(io_unit%log,'(a,i3,2(2x,3i4),2(2x,3f7.3))') &
    !    'guard_zones_t%restrict_simple2: iv, tl, tu =', iv, tl, tu
    if (present(same)) then
      call interpolate (source, target, [i1,i2,i3], iv, jt, pt, pi=pi, set=.true.)
      do i3=tl(3),tu(3); do i2=tl(2),tu(2); do i1=tl(1),tu(1)
        call interpolate (source, target, [i1,i2,i3], iv, jt, pt, pi=pi)
      end do; end do; end do
    else
      do i3=tl(3),tu(3); do i2=tl(2),tu(2); do i1=tl(1),tu(1)
        call interpolate (source, target, [i1,i2,i3], iv, jt, pt)
      end do; end do; end do
    end if
  end do
  if (detailed_timer) then
    cells = product(tu-tl+1)
    write (stdout,'(a,6p,2(2x,i6,i3),i6," cells,",f9.3," mus/cell")') &
      "guard_zones_t%different_opt: target, level, source, level,", &
      target%id, target%level, source%id, source%level, cells, &
      (wallclock()-start)/max(1,cells)
  end if
  call trace%end (itimer)
END SUBROUTINE different_opt

!===============================================================================
!> Compute interpolated guard zone values
!===============================================================================
SUBROUTINE restrict_vect (self, target, source, jt, pt, same)
  class(guard_zones_t):: self
  class(patch_t), pointer:: target, source
  integer:: jt(2)
  real(4):: pt(2)
  logical, optional:: same
  !.............................................................................
  real, dimension(:,:,:), pointer:: tmem, smem, td, sd
  class(mesh_t), pointer:: mesh
  integer:: i, i1, i2, i3, j, k, iv, tl(3), tu(3), sl(3), su(3), j0(3)
  real(8):: tlc(3), tuc(3), slc(3), suc(3)
  real(4):: p(3), p1, p2, p3, w(8)
  real(8):: used
  integer, save:: itimer=0
  !-----------------------------------------------------------------------------
  ! Map the source corners into target index space
  !-----------------------------------------------------------------------------
  call trace%begin ('guard_zones_t%restrict_vect', itimer=itimer)
  if (detailed_timer) &
    used = wallclock()
  !-----------------------------------------------------------------------------
  ! Interpolate in time and space
  !-----------------------------------------------------------------------------
  call source%time_interval (target%time, jt, pt)
  if (verbose > 0) &
    write(io_unit%log,'(a,2(2x,i6,i3,2x,3f7.3),2x,3f7.3)') &
      'target, level, lc, source, level, lc, relative =', &
      target%id, target%level, target%llc_cart, &
      source%id, source%level, source%llc_cart, &
      source%distance(target)/(0.5d0*(target%size+source%size))
  do iv=1,target%nv
    !---------------------------------------------------------------------------
    ! Target lower and upper corner in real space; these are the positions of
    ! the first and last target points.
    !---------------------------------------------------------------------------
    tlc = target%myposition (target%mesh%lb, iv)
    tuc = target%myposition (target%mesh%ub, iv)
    !---------------------------------------------------------------------------
    ! Map to source index space and limit to source interior. These are the last
    ! source points to the left of a target point, with the next source point to
    ! the right of the target point, and still inside the source interior.
    !---------------------------------------------------------------------------
    sl = source%index_only (tlc, 1, nowrap=.true.)
    su = source%index_only (tuc, 1, nowrap=.true.)
    sl = min(max(sl,source%mesh%li),source%mesh%ui-1)
    su = min(max(su,source%mesh%li),source%mesh%ui-1)
    !---------------------------------------------------------------------------
    ! Map those source points back to the target index space, rounding up to
    ! find the target points between the source points.
    !---------------------------------------------------------------------------
    slc = source%myposition(sl, iv)
    suc = source%myposition(su, iv)
    !---------------------------------------------------------------------------
    ! Compensate for periodic wrapping of corner points
    !---------------------------------------------------------------------------
    !where (target%position < source%position .and. slc > source%position)
      !slc = slc - target%mesh%b
    !end where
    !where (target%position > source%position .and. suc < source%position)
      !suc = suc + target%mesh%b
    !end where
    tl = target%index_only (slc, 1, roundup=.true., nowrap=.true.)
    tu = target%index_only (suc, 1, roundup=.true., nowrap=.true.)
    !---------------------------------------------------------------------------
    ! The brackets below should not really be needed, since the points should
    ! all be inside the target legal range
    !---------------------------------------------------------------------------
    tl = min(max(tl,target%mesh%lb),target%mesh%ub)
    tu = min(max(tu,target%mesh%lb),target%mesh%ub)
    if (verbose > 0) then
      write(io_unit%log,1) &
        'guard_zones_t%restrict_vect: iv, tlc, sl, slc, tl =', &
        iv, tlc, sl, slc, tl
      1 format(a,i3,2(2(2x,3f8.4,2x,3i4)))
      write(io_unit%log,1) &
        'guard_zones_t%restrict_vect: iv, tuc, su, suc, tu =', &
        iv, tuc, su, suc, tu
    end if
    !---------------------------------------------------------------------------
    ! Compute indices and interpolation weights for the first target point, and
    ! save the weights for the 8 nearest neighbor points in source space
    !---------------------------------------------------------------------------
    call source%index_space (target%myposition(tl), 1, j0, p)
    j = 0
    do i3=0,1
      p3 = merge (p(3), 1.0-p(3), i3==1)
      do i2=0,1
        p2 = merge (p(2), 1.0-p(2), i2==1)
        do i1=0,1
          p1 = merge (p(1), 1.0-p(1), i1==1)
          j = j+1
          w(j) = p1*p2*p3
        end do
      end do
    end do
    !---------------------------------------------------------------------------
    if (verbose > 0) &
      write(io_unit%log,'(a,i3,4(2x,3i4))') &
        'guard_zones_t%restrict_vect: iv, tl, tu, j0 =', &
        iv, tl, tu, j0
    !---------------------------------------------------------------------------
    ! Map target and source memory to contiguous 3-D arrays
    !---------------------------------------------------------------------------
    tmem => target%mem(:,:,:,iv,target%it,1)
    td => target%mem(:,:,:,target%idx%d,target%it,1)
    do k=1,2
      smem => source%mem(:,:,:,iv,jt(k),1)
      sd => source%mem(:,:,:,source%idx%d,jt(k),1)
      call update (target, tmem, td, tl, tu, iv, source, smem, sd, j0, k, &
                   pt(k), w, verbose)
    end do
  end do
  if (detailed_timer) then
    used = wallclock()-used
    write (stdout,'(a,6p,f6.3," mus/cell")') &
      "guard_zone_t%restrict_vect:", used/product(tu-tl+1)
  end if
  call trace%end (itimer)
END SUBROUTINE restrict_vect

!===============================================================================
!> Update interpolated guard zone values by looping over target indices,
!> incrementing source indices by two (FIXME: this assumes factor two AMR)
!> for each step in target indices.
!>
!> For a 16x16x16 test case, the expected source index ranges are:
!>  case  source   target   j0    l   u
!>  -----------------------------------
!>  left:   4-9      1-3     4    1   3
!> right:  14-19    19-22   14   19  22
!>  full:   4-19     4-11    4    4  11
!>  full:   4-19    12-19    4   12  19
!===============================================================================
SUBROUTINE update (target, tmem, td, l, u, iv, source, smem, sd, j0, k, &
                   pt, w, verbose)
  class(patch_t), pointer:: target, source
  real, dimension(:,:,:), contiguous, intent(out):: tmem
  real, dimension(:,:,:), contiguous, intent(in):: smem
  real, dimension(:,:,:), contiguous, intent(inout):: td, sd
  integer:: l(3), u(3), iv, j0(3), k, verbose
  real:: pt, w(8)
  !.............................................................................
  real(8):: wc(4)
  real, allocatable:: f(:,:), tot(:)
  integer:: i, j, n, i1 , i2, i3, j1, j2, j3, jj(3)
  !-----------------------------------------------------------------------------
  call trace%begin ('guard_zones_t%update')
  if (detailed_timer) then
    wc(1) = wallclock()
  end if
  n = product(u-l+1)
  allocate (f(n,8), tot(n))
  !-----------------------------------------------------------------------------
  ! Gather
  !-----------------------------------------------------------------------------
  i = 0
  j3 = j0(3)
  do i3=l(3),u(3)
    j2 = j0(2)
    do i2=l(2),u(2)
      j1 = j0(1)
      if (target%pervolume(iv)) then
        do i1=l(1),u(1)
          i = i+1
!define VALIDATE
#ifdef VALIDATE
          jj = [j1,j2,j3]
          if (any(jj < source%li .or. jj > source%ui)) then
            write(stderr     ,1) iv, jj, j0, l, u
            write(io_unit%log,1) iv, jj, j0, l, u
            1 format("guard_zones_t%update: ERROR iv,jj,j0,l,u =",i3,4(2x,3i4))
          end if
#endif
          f(i,1) = smem(j1  ,j2  ,j3  )/sd(j1  ,j2  ,j3  )
          f(i,2) = smem(j1+1,j2  ,j3  )/sd(j1+1,j2  ,j3  )
          f(i,3) = smem(j1  ,j2+1,j3  )/sd(j1  ,j2+1,j3  )
          f(i,4) = smem(j1+1,j2+1,j3  )/sd(j1+1,j2+1,j3  )
          f(i,5) = smem(j1  ,j2  ,j3+1)/sd(j1  ,j2  ,j3+1)
          f(i,6) = smem(j1+1,j2  ,j3+1)/sd(j1+1,j2  ,j3+1)
          f(i,7) = smem(j1  ,j2+1,j3+1)/sd(j1  ,j2+1,j3+1)
          f(i,8) = smem(j1+1,j2+1,j3+1)/sd(j1+1,j2+1,j3+1)
          j1 = j1+2
        end do
      else
        do i1=l(1),u(1)
          i = i+1
#ifdef VALIDATE
          jj = [j1,j2,j3]
          if (any(jj < source%li .or. jj > source%ui)) then
            write(stderr     ,1) iv, jj, j0, l, u
            write(io_unit%log,1) iv, jj, j0, l, u
          end if
#endif
          f(i,1) = smem(j1  ,j2  ,j3  )
          f(i,2) = smem(j1+1,j2  ,j3  )
          f(i,3) = smem(j1  ,j2+1,j3  )
          f(i,4) = smem(j1+1,j2+1,j3  )
          f(i,5) = smem(j1  ,j2  ,j3+1)
          f(i,6) = smem(j1+1,j2  ,j3+1)
          f(i,7) = smem(j1  ,j2+1,j3+1)
          f(i,8) = smem(j1+1,j2+1,j3+1)
          j1 = j1+2
        end do
      end if
      j2 = j2+2
    end do
    j3 = j3+2
  end do
  if (detailed_timer) then
    wc(2) = wallclock()
  end if
  !-----------------------------------------------------------------------------
  ! Compute weighted results, vectorized
  !-----------------------------------------------------------------------------
  tot = 0.0
  if (target%unsigned(iv)) then
    do j=1,8
    do i=1,n
      tot(i) = tot(i) + pt*w(j)*log(f(i,j))
    end do
    end do
    tot = exp(tot)
  else
    do j=1,8
    do i=1,n
      tot(i) = tot(i) + pt*w(j)*f(i,j)
    end do
    end do
  end if
  if (detailed_timer) then
    wc(3) = wallclock()
  end if
  !-----------------------------------------------------------------------------
  ! Scatter, vectorized
  !-----------------------------------------------------------------------------
  i = 0
  do i3=l(3),u(3)
    do i2=l(2),u(2)
      if (k==1) then
        do i1=l(1),u(1)
          i = i+1
          tmem(i1,i2,i3) = tot(i)
        end do
      else
        do i1=l(1),u(1)
          i = i+1
          tmem(i1,i2,i3) = tmem(i1,i2,i3) + tot(i)
        end do
      end if
    end do
  end do
  !-----------------------------------------------------------------------------
  ! For per-volume quantities, multiply with mass density
  !-----------------------------------------------------------------------------
  if (target%pervolume(iv) .and. k==2) then
    do i3=l(3),u(3)
      do i2=l(2),u(2)
        do i1=l(1),u(1)
          tmem(i1,i2,i3) = tmem(i1,i2,i3)*td(i1,i2,i3)
        end do
      end do
    end do
  end if
  deallocate (f, tot)
  if (detailed_timer) then
    wc(4) = wallclock()
    do i=1,3
      !$omp atomic update
      update_t(i) = update_t(i) + (wc(i+1)-wc(i))
    end do
  end if
  call trace%end()
END SUBROUTINE update

!===============================================================================
!> Interpolate in space and time, for a case where it is safe = the source point
!> is inside the valid domain (inside a distance of half a cell from the source
!> patch boundary for non-staggered points)
!===============================================================================
SUBROUTINE interpolate (source, target, ii, iv, jt, pt, pi, set)
  class(patch_t)           :: source
  class(patch_t), pointer  :: target
  integer                  :: ii(3), iv, jt(2)
  real                     :: pt(2)
  real, optional           :: pi(:)
  logical, optional        :: set
  !.............................................................................
  real(8)                  :: pos(3)
  integer                  :: j(3)
  real                     :: p(3)
  integer, save            :: itimer=0
  !-----------------------------------------------------------------------------
  if (present(set)) then
    pi(01) = pt(1)*(1.-p(3))*(1.-p(2))*(1.0-p(1))
    pi(02) = pt(1)*(1.-p(3))*(1.-p(2))*(    p(1))
    pi(03) = pt(1)*(1.-p(3))*(   p(2))*(1.0-p(1))
    pi(04) = pt(1)*(1.-p(3))*(   p(2))*(    p(1))
    pi(05) = pt(1)*(   p(3))*(1.-p(2))*(1.0-p(1))
    pi(06) = pt(1)*(   p(3))*(1.-p(2))*(    p(1))
    pi(07) = pt(1)*(   p(3))*(   p(2))*(1.0-p(1))
    pi(08) = pt(1)*(   p(3))*(   p(2))*(    p(1))
    pi(09) = pt(2)*(1.-p(3))*(1.-p(2))*(1.0-p(1))
    pi(10) = pt(2)*(1.-p(3))*(1.-p(2))*(    p(1))
    pi(11) = pt(2)*(1.-p(3))*(   p(2))*(1.0-p(1))
    pi(12) = pt(2)*(1.-p(3))*(   p(2))*(    p(1))
    pi(13) = pt(2)*(   p(3))*(1.-p(2))*(1.0-p(1))
    pi(14) = pt(2)*(   p(3))*(1.-p(2))*(    p(1))
    pi(15) = pt(2)*(   p(3))*(   p(2))*(1.0-p(1))
    pi(16) = pt(2)*(   p(3))*(   p(2))*(    p(1))
  else if (source%pervolume(iv)) then
    call interpolate_specific (source, target, ii, iv, jt, pt)
  else
    pos = target%myposition (ii, iv)
    call source%index_space (pos, iv, j, p)
    !if (io%verbose > 1) &
    !  write(io_unit%log,'(a,2i6,2x,3i4,2x,3f7.3)') &
    !    'interpolate: target, source, j, p =', target%id, source%id, j, p
    !call trace%begin ('patch_t%interpolate', itimer=itimer)
    if (source%unsigned(iv)) then
     if (present(pi)) then
      target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = exp( &
        pi(01)*log(source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(1),1)) + &
        pi(02)*log(source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(1),1)) + &
        pi(03)*log(source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(1),1)) + &
        pi(04)*log(source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(1),1)) + &
        pi(05)*log(source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(1),1)) + &
        pi(06)*log(source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(1),1)) + &
        pi(07)*log(source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(1),1)) + &
        pi(08)*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1)) + &
        pi(09)*log(source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(2),1)) + &
        pi(10)*log(source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(2),1)) + &
        pi(11)*log(source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(2),1)) + &
        pi(12)*log(source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(2),1)) + &
        pi(13)*log(source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(2),1)) + &
        pi(14)*log(source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(2),1)) + &
        pi(15)*log(source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(2),1)) + &
        pi(16)*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1)))
     else
      target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = exp( &
       pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(1),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(1),1)))   + &
                         (   p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(1),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(1),1))))  + &
              (   p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(1),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(1),1)))   + &
                         (   p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(1),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1))))) + &
       pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(2),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(2),1)))   + &
                         (   p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(2),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(2),1))))  + &
              (   p(3))*((1.-p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(2),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(2),1)))   + &
                         (   p(2))*((1.0-p(1))*log(source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(2),1))    + &
                                    (    p(1))*log(source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))))
     end if
    else
     if (present(pi)) then
      target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
        pi(01)*source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(1),1) + &
        pi(02)*source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(1),1) + &
        pi(03)*source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(1),1) + &
        pi(04)*source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(1),1) + &
        pi(05)*source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(1),1) + &
        pi(06)*source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(1),1) + &
        pi(07)*source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(1),1) + &
        pi(08)*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1) + &
        pi(09)*source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(2),1) + &
        pi(10)*source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(2),1) + &
        pi(11)*source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(2),1) + &
        pi(12)*source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(2),1) + &
        pi(13)*source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(2),1) + &
        pi(14)*source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(2),1) + &
        pi(15)*source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(2),1) + &
        pi(16)*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1)
     else
      target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
       pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(1),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(1),1))   + &
                         (   p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(1),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(1),1)))  + &
              (   p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(1),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(1),1))   + &
                         (   p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(1),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(1),1)))) + &
       pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)  ,j(3)  ,iv,jt(2),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)  ,j(3)  ,iv,jt(2),1))   + &
                         (   p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)+1,j(3)  ,iv,jt(2),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)+1,j(3)  ,iv,jt(2),1)))  + &
              (   p(3))*((1.-p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)  ,j(3)+1,iv,jt(2),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)  ,j(3)+1,iv,jt(2),1))   + &
                         (   p(2))*((1.0-p(1))*source%mem(j(1)  ,j(2)+1,j(3)+1,iv,jt(2),1)    + &
                                    (    p(1))*source%mem(j(1)+1,j(2)+1,j(3)+1,iv,jt(2),1))))
     end if
    end if
  end if
  if (target%id==io%id_debug .and. io%verbose>2) &
    write(io%output,1) target%id, ii, iv, source%id, j, jt, p, pt, &
      target%mem(ii(1),ii(2),ii(3),iv,target%it,1)
    1 format(i6,1x,3i3,i4,2x,"DBG patch_t%interpolate src:",i6, &
      "   j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2,"   out:",1p,e15.6)
  !call trace%end (itimer)
END SUBROUTINE interpolate

!===============================================================================
!> Interpolate in space and time, for a case where it is safe = the source point
!> is inside the valid domain (inside a distance of half a cell from the source
!> patch boundary for non-staggered points)
!===============================================================================
SUBROUTINE interpolate_specific (source, target, ii, iv, jt, pt)
  class(patch_t)           :: source
  class(patch_t), pointer  :: target
  integer                  :: ii(3), iv, jt(2)
  real                     :: pt(2)
  real                     :: tmp(2,2,2,2)
  !.............................................................................
  real(8)                  :: pos(3)
  integer                  :: j(3)
  real                     :: p(3)
  integer, save            :: itimer=0
  !-----------------------------------------------------------------------------
  pos = target%myposition (ii, iv)
  call source%index_space (pos, iv, j, p)
  !call trace%begin ('patch_t%interpolate', itimer=itimer)
  !-----------------------------------------------------------------------------
  ! Pick up 16 values for tri-linear interpolation, divided by density
  !-----------------------------------------------------------------------------
  associate (id=>source%idx%d)
  tmp(:,:,:,:) = source%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,iv,jt(1:2),1) &
               / source%mem(j(1):j(1)+1,j(2):j(2)+1,j(3):j(3)+1,id,jt(1:2),1)
  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
   pt(1)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,1)    + &
                                (    p(1))*tmp(2,1,1,1))   + &
                     (   p(2))*((1.0-p(1))*tmp(1,2,1,1)    + &
                                (    p(1))*tmp(2,2,1,1)))  + &
          (   p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,1)    + &
                                (    p(1))*tmp(2,1,2,1))   + &
                     (   p(2))*((1.0-p(1))*tmp(1,2,2,1)    + &
                                (    p(1))*tmp(2,2,2,1)))) + &
   pt(2)*((1.-p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,1,2)    + &
                                (    p(1))*tmp(2,1,1,2))   + &
                     (   p(2))*((1.0-p(1))*tmp(1,2,1,2)    + &
                                (    p(1))*tmp(2,2,1,2)))  + &
          (   p(3))*((1.-p(2))*((1.0-p(1))*tmp(1,1,2,2)    + &
                                (    p(1))*tmp(2,1,2,2))   + &
                     (   p(2))*((1.0-p(1))*tmp(1,2,2,2)    + &
                                (    p(1))*tmp(2,2,2,2))))
  !-----------------------------------------------------------------------------
  ! Multiply by target denssity
  !-----------------------------------------------------------------------------
  target%mem(ii(1),ii(2),ii(3),iv,target%it,1) = &
    target%mem(ii(1),ii(2),ii(3),iv,target%it,1) * &
    target%mem(ii(1),ii(2),ii(3),id,target%it,1)
  end associate
  if (target%id==io%id_debug .and. io%verbose>2) &
    write(io%output,1) target%id, ii, iv, source%id, j, jt, p, pt, &
      target%mem(ii(1),ii(2),ii(3),iv,target%it,1)
    1 format(i6,1x,3i3,i4,2x,"DBG patch_t%interpolate_specific src:", &
      i6,"   j,p;",3i4,1x,2i4,3f6.2,1x,2f6.2,"   out:",1p,e15.6)
  !call trace%end (itimer)
END SUBROUTINE interpolate_specific

SUBROUTINE info
  if (verbose > 0) then
    write (stdout,*) 'guard_zones_t%info: update-t =', update_t
  end if
END SUBROUTINE info

END MODULE guard_zones_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
Community Relations works closely with communities throughout the region on issues that affect the quality of life. These issues range from economic development and K12 education to diversity and leadership development.
Community Relations:

     Works with communitybased organizations to enhance access to UC Davis both as an educational institution and as an informational resource.
     Forges a closer working relationship between UC Davis campus leaders and local and regional community groups in addressing a myriad of social, economic, education, and environmental issues.
     Fosters partnerships and identifies areas of collaboration with a wide variety of community organizations.
     Highlights issues of emerging or growing concern in the region for UC Davis officials and faculty.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! Compile with the Intel compiler using
! ifort -coarray=shared -coarray-num-image=3 <file-name>
! to set the number of images at compile-time.
! Override the number of images with an environment variable FOR_COARRAY_NUM_IMAGES
program main
  implicit none
  real, allocatable :: array(:)[:]

  print *,"Image ", this_image()," before allocate"
  allocate(array(10)[*]) ! implicit synchronization
  print *,"Image ", this_image()," after allocate"
end program main

max_stars_repo_name: sourceryinstitute/miniFAVOR

max_issues_repo_name: sourceryinstitute/miniFAVOR

max_forks_repo_name: sourceryinstitute/miniFAVOR


text:
module initdata_module

  use amrex_paralleldescriptor_module, only: parallel_IOProcessor => amrex_pd_ioprocessor
  use network, only: nspec
  use amrex_fort_module, only : amrex_spacedim
  use base_state_geometry_module, only: nr_fine, max_radial_level
  use meth_params_module, only: nscal, rho_comp, rhoh_comp, temp_comp, spec_comp, pi_comp, &
                                prob_lo, prob_hi, grav_const
  use eos_module
  use eos_type_module
  use amrex_constants_module
  use probin_module, only : pert_amp, scale_height, pres_base, k_hoz, k_vert

  implicit none

  private

contains

  subroutine initdata(lev, time, lo, hi, &
       scal, scal_lo, scal_hi, nc_s, &
       vel, vel_lo, vel_hi, nc_v, &
       s0_init, p0_init, dx) bind(C, name="initdata")

    integer         , intent(in   ) :: lev, lo(3), hi(3)
    integer         , intent(in   ) :: scal_lo(3), scal_hi(3), nc_s
    integer         , intent(in   ) :: vel_lo(3), vel_hi(3), nc_v
    double precision, intent(in   ) :: time
    double precision, intent(inout) :: scal(scal_lo(1):scal_hi(1), &
         scal_lo(2):scal_hi(2), &
         scal_lo(3):scal_hi(3), 1:nc_s)
    double precision, intent(inout) :: vel(vel_lo(1):vel_hi(1), &
         vel_lo(2):vel_hi(2), &
         vel_lo(3):vel_hi(3), 1:nc_v)
    double precision, intent(in   ) :: s0_init(0:max_radial_level,0:nr_fine-1,1:nscal)
    double precision, intent(in   ) :: p0_init(0:max_radial_level,0:nr_fine-1)
    double precision, intent(in   ) :: dx(3)

    ! local variables

    integer          :: i,j,k,r
    type (eos_t) :: eos_state
    double precision :: x, y, rho0, rho_local, rho_base

!    ! abort program
!    call amrex_error()

    ! set velocity to zero
    vel(vel_lo(1):vel_hi(1),vel_lo(2):vel_hi(2),vel_lo(3):vel_hi(3),1:nc_v) = 0.d0

    do k=lo(3),hi(3)
       do j=lo(2),hi(2)
          do i=lo(1),hi(1)

             if (amrex_spacedim .eq. 1) then
                r = i
             else if (amrex_spacedim .eq. 2) then
                r = j
             else if (amrex_spacedim .eq. 3) then
                r = k
             end if

             ! set scalars using s0
             scal(i,j,k,rho_comp)  = s0_init(lev,r,rho_comp)
             scal(i,j,k,rhoh_comp) = s0_init(lev,r,rhoh_comp)
             scal(i,j,k,temp_comp) = s0_init(lev,r,temp_comp)
             scal(i,j,k,spec_comp:spec_comp+nspec-1) = &
                  s0_init(lev,r,spec_comp:spec_comp+nspec-1)

             ! initialize pi to zero for now
             scal(i,j,k,pi_comp) = 0.d0

          end do
       end do
    end do

    ! introduce density fluctuations - only implemented for 2D

    if ((amrex_spacedim == 1) .or. amrex_spacedim == 3) then
      call amrex_error("density perturbation only implemented for 2d")
    endif

    do k=lo(3),hi(3)
    do j=lo(2),hi(2)
    do i=lo(1),hi(1)

      r = j ! this loop is currently 2d only
      x = prob_lo(1) + (dble(i)+0.5d0) * dx(1)
      y = prob_lo(2) + (dble(j)+0.5d0) * dx(2)

      ! This seems to work ok with sealed box?

!      rho0 = s0_init(lev,r,rho_comp)
!      rho_local = rho0 * (1.0 + &
!        pert_amp * exp(-y/(2.0_num * scale_height)) * &
!        cos(x * k_hoz * M_PI / (prob_hi(1) - prob_lo(1)) ) * &
!        sin(y * k_vert * M_PI / (prob_hi(2) - prob_lo(2)) )  &
!      )

      ! What about this ?

      rho_base = pres_base / scale_height / abs(grav_const)
      rho0 = s0_init(lev,r,rho_comp)
      rho_local = rho0 &
        + rho_base * pert_amp * exp(-y/(2d0 * scale_height)) * &
        cos(x * k_hoz * m_pi / (prob_hi(1) - prob_lo(1)) ) * &
        sin(y * k_vert * m_pi / (prob_hi(2) - prob_lo(2)) )

      ! if k_vert is 0, dont multiply by sin(0) = 0
      ! (make this all more elegant later)

      if (k_vert == 0d0) then
        rho_local = rho0 &
          + rho_base * pert_amp * exp(-y/(2d0 * scale_height)) * &
          cos(x * k_hoz * m_pi / (prob_hi(1) - prob_lo(1)) )
      endif

      !
      eos_state%rho   = rho_local
      eos_state%p     = p0_init(lev,r)
      eos_state%T     = s0_init(lev,r,temp_comp)
      eos_state%xn(:) = 1.0d0 ! single fluid

      ! (rho,p) --> T, h
      call eos(eos_input_rp, eos_state)

      scal(i,j,k,rho_comp)  = eos_state%rho
      scal(i,j,k,rhoh_comp) = eos_state%rho * eos_state%h
      scal(i,j,k,temp_comp) = eos_state%T
      scal(i,j,k,spec_comp:spec_comp+nspec-1) = eos_state%rho * eos_state%xn ! single fluid

    end do
    end do
    end do

  end subroutine initdata

  subroutine initdata_sphr(time, lo, hi, &
       scal, scal_lo, scal_hi, nc_s, &
       vel, vel_lo, vel_hi, nc_v, &
       s0_init, p0_init, &
       dx, r_cc_loc, r_edge_loc, &
       cc_to_r, ccr_lo, ccr_hi) bind(C, name="initdata_sphr")

    integer         , intent(in   ) :: lo(3), hi(3)
    integer         , intent(in   ) :: scal_lo(3), scal_hi(3), nc_s
    integer         , intent(in   ) :: vel_lo(3), vel_hi(3), nc_v
    integer         , intent(in   ) :: ccr_lo(3), ccr_hi(3)
    double precision, intent(in   ) :: time
    double precision, intent(inout) :: scal(scal_lo(1):scal_hi(1), &
         scal_lo(2):scal_hi(2), &
         scal_lo(3):scal_hi(3), nc_s)
    double precision, intent(inout) :: vel(vel_lo(1):vel_hi(1), &
         vel_lo(2):vel_hi(2), &
         vel_lo(3):vel_hi(3), nc_v)
    double precision, intent(in   ) :: s0_init(0:max_radial_level,0:nr_fine-1,1:nscal)
    double precision, intent(in   ) :: p0_init(0:max_radial_level,0:nr_fine-1)
    double precision, intent(in   ) :: dx(3)
    double precision, intent(in   ) :: r_cc_loc  (0:max_radial_level,0:nr_fine-1)
    double precision, intent(in   ) :: r_edge_loc(0:max_radial_level,0:nr_fine)
    double precision, intent(in   ) :: cc_to_r   (ccr_lo(1):ccr_hi(1), &
         ccr_lo(2):ccr_hi(2),ccr_lo(3):ccr_hi(3))

    call amrex_error("error: called initdata_sphr, but is a stub")

  end subroutine initdata_sphr

end module initdata_module

max_stars_repo_name: pgrete/MAESTROeX

max_issues_repo_name: pgrete/MAESTROeX

max_forks_repo_name: pgrete/MAESTROeX


 MODULE module_ltng_cpmpr92z
 CONTAINS
 SUBROUTINE ltng_cpmpr92z ( &
                            dx, dy, xland, ht, z, t, &
                            kLNB, &
                            cldtop_adjustment, &
                            msft, &
                            cldfra_cup, &
                            htop, &
                            shall, &
                            ids, ide, jds, jde, kds, kde, &
                            ims, ime, jms, jme, kms, kme, &
                            ips, ipe, jps, jpe, kps, kpe, &
                            total_flashrate &
                          )
 USE module_state_description
 USE module_model_constants
 USE module_wrf_error
 IMPLICIT NONE
 REAL, INTENT(IN ) :: dx, dy
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: xland, ht
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: z, t
 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: kLNB
 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN ) :: cldfra_cup
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: htop
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: shall
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: msft
 REAL, INTENT(IN ) :: cldtop_adjustment
 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
 INTEGER, INTENT(IN ) :: ips,ipe, jps,jpe, kps,kpe
 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( OUT) :: total_flashrate
 REAL :: dA
 REAL :: zkm
 REAL :: cldfra_scaling_f
 REAL :: mapfac_scaling_f
 REAL, PARAMETER:: baseArea=1296.
 INTEGER :: i,k,j
 INTEGER :: kbtm, ktop
 CHARACTER (LEN=250) :: message
 dA = dx*dy/1E6
 cldfra_scaling_f=0.
 mapfac_scaling_f=1.
 total_flashrate( ips:ipe,jps:jpe ) = 0.
 jloop: DO j=jps,jpe
    iloop: DO i=ips,ipe
        mapfac_scaling_f = 1.0 / (msft(i,j) * msft(i,j))
        ktop = kLNB(i,j)
        IF (shall(i,j)>0.5) THEN
          ktop = -1
        ELSE
          ktop=nint(htop(i,j))
          k = kps
          DO WHILE( cldfra_cup(i,k,j) .le. 0.01 .and. k .le. ktop )
            k = k+1
          ENDDO
          kbtm = k
          IF ( (ktop - kbtm) > 0 ) THEN
            cldfra_scaling_f=SUM(cldfra_cup(i,kbtm:ktop,j))/REAL(ktop-kbtm+1)
            WRITE(message, * ) ' lightning_cpm: scaling by average cloud fraction (kbotom, ktop) = ', cldfra_scaling_f, kbtm, ktop
            CALL wrf_debug ( 15, message )
            IF ( t(i,ktop,j) .lt. 273.15 .and. &
              ktop .ge. kps .and. ktop .le. kpe ) THEN
              zkm = ( z(i,ktop,j) - ht(i,j) )/1E3 + cldtop_adjustment
              IF ( zkm .gt. 0. ) THEN
                IF ( xland(i,j) .lt. 1.5 ) THEN
                  total_flashrate(i,j) = 3.44E-5 * (zkm**4.9) /60.
                ELSE
                  total_flashrate(i,j) = 6.57E-6 * (zkm**4.9) /60.
                ENDIF
                total_flashrate(i,j) = total_flashrate(i,j) * cldfra_scaling_f * mapfac_scaling_f
              ENDIF
            ENDIF
          ENDIF
        ENDIF
    ENDDO iloop
 ENDDO jloop
 total_flashrate(ips:ipe,jps:jpe) = total_flashrate(ips:ipe,jps:jpe) * dA/baseArea
 END SUBROUTINE ltng_cpmpr92z
 END MODULE module_ltng_cpmpr92z


!-------------------------------------------------------------------------------
!> module Copy topography
!!
!! @par Description
!!          subroutines for preparing topography data (copy from parent domain)
!!
!! @author Team SCALE
!<
!-------------------------------------------------------------------------------
module mod_copytopo
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use gtool_file_h
  use gtool_file, only: &
     FileRead
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: COPYTOPO

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  integer, public            :: CNVTOPO_TYPE = -1

  integer, public, parameter :: I_IGNORE     =  0
  integer, public, parameter :: I_GTOPO30    =  1
  integer, public, parameter :: I_DEM50M     =  2
  integer, public, parameter :: I_GMTED2010  =  3

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: COPYTOPO_transgrid
  private :: COPYTOPO_setalpha
  private :: COPYTOPO_input_data
  private :: COPYTOPO_mix_data

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, private, parameter :: handle = 1
  integer, private            :: itp_nh = 4

  character(len=H_LONG), private :: COPYTOPO_IN_BASENAME  = ''

  real(RP), private :: COPYTOPO_TRANSITION_DX = -1.0_RP  !< thickness of transition region [m]: x
  real(RP), private :: COPYTOPO_TRANSITION_DY = -1.0_RP  !< thickness of transition region [m]: y
  real(RP), private :: COPYTOPO_TRANSFACT     = -1.0_RP  !< stretch factor of transition region
  real(RP), private :: COPYTOPO_FRACX         =  1.0_RP  !< fraction of transition region (x) [0-1]
  real(RP), private :: COPYTOPO_FRACY         =  1.0_RP  !< fraction of transition region (y) [0-1]
  real(RP), private :: COPYTOPO_taux          =  1.0_RP  !< maximum value for mixing tau (x) [s]
  real(RP), private :: COPYTOPO_tauy          =  1.0_RP  !< maximum value for mixing tau (y) [s]

  logical,  private :: COPYTOPO_ENTIRE_REGION = .false.  !< copy parent topo over an entire region
  logical,  private :: COPYTOPO_LINEAR_H      = .true.   !< linear or non-linear profile of relax region
  real(RP), private :: COPYTOPO_EXP_H         = 2.0_RP   !< factor of non-linear profile of relax region

  real(RP), private, allocatable :: CTRX(:)              !< center buffer factor [0-1]: x
  real(RP), private, allocatable :: CTRY(:)              !< center buffer factor [0-1]: y
  real(RP), private, allocatable :: COPYTOPO_alpha(:,:)  !> damping coefficient  [0-1]
  real(RP), private, allocatable :: topo_pd(:,:)         !> topography of parent domain

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup and Main
  subroutine COPYTOPO( &
      topo_cd          ) ! [inout]
    use scale_process, only: &
       PRC_MPIstop
    use scale_const, only: &
       D2R   => CONST_D2R
    use scale_grid, only: &
       GRID_CDZ,   &
       DX,         &
       DY,         &
       BUFFER_DX,  &
       BUFFER_DY,  &
       BUFFFACT
    use scale_grid_nest, only: &
       NEST_INTERP_LEVEL
    implicit none

    real(RP), intent(inout) :: topo_cd(:,:)

    NAMELIST / PARAM_COPYTOPO / &
       COPYTOPO_IN_BASENAME,    &
       COPYTOPO_TRANSITION_DX,  &
       COPYTOPO_TRANSITION_DY,  &
       COPYTOPO_TRANSFACT,      &
       COPYTOPO_FRACX,          &
       COPYTOPO_FRACY,          &
       COPYTOPO_taux,           &
       COPYTOPO_tauy,           &
       COPYTOPO_ENTIRE_REGION,  &
       COPYTOPO_LINEAR_H,       &
       COPYTOPO_EXP_H

    integer :: ierr
    !---------------------------------------------------------------------------

    if( IO_L ) write(IO_FID_LOG,*) ''
    if( IO_L ) write(IO_FID_LOG,*) '+++ Module[COPYTOPO]/Categ[COPYTOPO]'

    !--- read namelist
    rewind(IO_FID_CONF)
    read(IO_FID_CONF,nml=PARAM_COPYTOPO,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( IO_L ) write(IO_FID_LOG,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_COPYTOPO. Check!'
       call PRC_MPIstop
    endif
    if( IO_LNML ) write(IO_FID_LOG,nml=PARAM_COPYTOPO)

    if ( COPYTOPO_TRANSITION_DX < 0.0_RP ) then
       COPYTOPO_TRANSITION_DX = BUFFER_DX
    end if
    if ( COPYTOPO_TRANSITION_DY < 0.0_RP ) then
       COPYTOPO_TRANSITION_DY = BUFFER_DY
    end if
    if ( COPYTOPO_TRANSFACT < 0.0_RP ) then
       COPYTOPO_TRANSFACT = BUFFFACT
    end if

    allocate( CTRX(IA) )
    allocate( CTRY(JA) )
    allocate( COPYTOPO_alpha(IA,JA) )
    allocate( topo_pd       (IA,JA) )
    COPYTOPO_alpha(:,:) = 0.0_RP

    itp_nh = int( NEST_INTERP_LEVEL )

    ! copy topography from parent domain to transition region
    call COPYTOPO_transgrid

    call COPYTOPO_setalpha

    call COPYTOPO_input_data( topo_pd ) ! (out)

    call COPYTOPO_mix_data( topo_cd, &  ! (inout)
                            topo_pd  )  ! (in)

    return
  end subroutine COPYTOPO

  !-----------------------------------------------------------------------------
  !> Individual Procedures
  !-----------------------------------------------------------------------------

  !-----------------------------------------------------------------------------
  !> Generate transition grid
  subroutine COPYTOPO_transgrid
    use scale_process, only: &
       PRC_MPIstop, &
       PRC_myrank
    use scale_les_process, only: &
       PRC_2Drank,  &
       PRC_NUM_X,   &
       PRC_NUM_Y
    use scale_const, only: &
       RADIUS => CONST_RADIUS, &
       EPS    => CONST_EPS,    &
       D2R    => CONST_D2R
    use scale_grid, only: &
       DX,                  &
       DY,                  &
       CXG   => GRID_CXG,   &
       FXG   => GRID_FXG,   &
       CYG   => GRID_CYG,   &
       FYG   => GRID_FYG,   &
       CBFXG => GRID_CBFXG, &
       CBFYG => GRID_CBFYG, &
       BUFFFACT
    implicit none

    real(RP), allocatable :: CTRXG(:) !< center buffer factor [0-1]: x, global
    real(RP), allocatable :: CTRYG(:) !< center buffer factor [0-1]: y, global

    real(RP), allocatable :: buffx(:),  buffy(:)
    real(RP)              :: bufftotx,  bufftoty
    real(RP), allocatable :: transx(:), transy(:)
    real(RP)              :: transtotx, transtoty

    integer :: IAG, JAG
    integer :: imain, ibuff, itrans
    integer :: jmain, jbuff, jtrans
    integer :: copy_is, copy_ie, copy_js, copy_je
    integer :: i, j, ii, jj, t
    !---------------------------------------------------------------------------

    ! array size (global domain)
    IAG = IHALO + IMAX*PRC_NUM_X + IHALO
    JAG = JHALO + JMAX*PRC_NUM_Y + JHALO

    allocate( buffx (0:IAG) )
    allocate( buffy (0:JAG) )
    allocate( transx(0:IAG) )
    allocate( transy(0:JAG) )
    allocate( CTRXG (  IAG) )
    allocate( CTRYG (  JAG) )

    ! X-direction
    ! calculate buffer grid size
    buffx(:)  = DX
    transx(0) = DX
    bufftotx  = 0.0_RP
    transtotx = 0.0_RP

    do i = IHALO+1, IAG
       if( abs(CBFXG(i) - 0.0_RP) < EPS ) exit
       buffx(i) = buffx(i-1) * BUFFFACT
       bufftotx = bufftotx + buffx(i)
    enddo
    ibuff = i - (IHALO+1)

    do i = 1, IAG
       if( transtotx >= COPYTOPO_TRANSITION_DX ) exit
       transx(i) = transx(i-1) * COPYTOPO_TRANSFACT
       transtotx = transtotx + transx(i)
    enddo
    itrans = i - 1
    imain  = IAG - 2*ibuff - 2*itrans - 2*IHALO

    if ( imain < 1 ) then
       write(*,*) 'xxx Not appropriate transition width for global domain(X).', COPYTOPO_TRANSITION_DX
       write(*,*) '    # of buffer region (one side)', ibuff
       write(*,*) '    # of transion region (one side)', itrans
       call PRC_MPIstop
    endif

    ! calc transition factor (global domaim)
    CTRXG(:) = 0.0_RP
    do i = 1, IHALO+ibuff
       CTRXG(i) = 1.0_RP
    enddo

    if ( itrans > 0 ) then
       copy_is = IHALO+ibuff+1
       copy_ie = IHALO+ibuff+itrans
       do i = copy_is, copy_ie
          CTRXG(i) = (transtotx+bufftotx+FXG(IHALO    )-CXG(i)) / transtotx
       enddo
       copy_is = IHALO+ibuff+itrans+imain+1
       copy_ie = IHALO+ibuff+itrans+imain+itrans+ibuff
       do i = copy_is, copy_ie
          CTRXG(i) = (transtotx+bufftotx-FXG(IAG-IHALO)+CXG(i)) / transtotx
       enddo
    endif

    copy_is = IHALO+ibuff+itrans+imain+itrans+ibuff+1
    copy_ie = IHALO+ibuff+itrans+imain+itrans+ibuff+IHALO
    do i = copy_is, copy_ie
       CTRXG(i) = 1.0_RP
    enddo
    CTRXG(:) = max( min( CTRXG(:), 1.0_RP ), 0.0_RP )

    ! Y-direction
    ! calculate buffer grid size
    buffy(:)  = DY
    transy(0) = DY
    bufftoty  = 0.0_RP
    transtoty = 0.0_RP

    do j = JHALO+1, JAG
       if( abs(CBFYG(j) - 0.0_RP) < EPS ) exit
       buffy(j) = buffy(j-1) * BUFFFACT
       bufftoty = bufftoty + buffy(j)
    enddo
    jbuff = j - (JHALO+1)

    do j = 1, JAG
       if( transtoty >= COPYTOPO_TRANSITION_DY ) exit
       transy(j) = transy(j-1) * COPYTOPO_TRANSFACT
       transtoty = transtoty + transy(j)
    enddo
    jtrans = j - 1
    jmain  = JAG - 2*jbuff - 2*jtrans - 2*JHALO

    if ( jmain < 1 ) then
       write(*,*) 'xxx Not appropriate transition width for global domain(Y).', COPYTOPO_TRANSITION_DY
       write(*,*) '    # of buffer region (one side)', jbuff
       write(*,*) '    # of transion region (one side)', jtrans
       call PRC_MPIstop
    endif

    ! calc transition factor (global domaim)
    CTRYG(:) = 0.0_RP
    do j = 1, JHALO+jbuff
       CTRYG(j) = 1.0_RP
    enddo

    if ( jtrans > 0 ) then
       copy_js = JHALO+jbuff+1
       copy_je = JHALO+jbuff+jtrans
       do j = copy_js, copy_je
          CTRYG(j) = (transtoty+bufftoty+FYG(JHALO    )-CYG(j)) / transtoty
       enddo
       copy_js = JHALO+jbuff+jtrans+jmain+1
       copy_je = JHALO+jbuff+jtrans+jmain+jtrans+jbuff
       do j = copy_js, copy_je
          CTRYG(j) = (transtoty+bufftoty-FYG(JAG-JHALO)+CYG(j)) / transtoty
       enddo
    endif

    copy_js = JHALO+jbuff+jtrans+jmain+jtrans+jbuff+1
    copy_je = JHALO+jbuff+jtrans+jmain+jtrans+jbuff+JHALO
    do j = copy_js, copy_je
       CTRYG(j) = 1.0_RP
    enddo
    CTRYG(:) = max( min( CTRYG(:), 1.0_RP ), 0.0_RP )

    ! horizontal coordinate (local domaim)
    do i = 1, IA
       ii = i + PRC_2Drank(PRC_myrank,1) * IMAX
       CTRX(i) = CTRXG(ii)
    enddo
    do j = 1, JA
       jj = j + PRC_2Drank(PRC_myrank,2) * JMAX
       CTRY(j) = CTRYG(jj)
    enddo

    deallocate( transx )
    deallocate( transy )

    return
  end subroutine COPYTOPO_transgrid

  !-----------------------------------------------------------------------------
  !> Calc dumping coefficient alpha
  subroutine COPYTOPO_setalpha
    use scale_const, only: &
       EPS => CONST_EPS, &
       PI  => CONST_PI
    use scale_grid, only: &
       CBFX => GRID_CBFX, &
       CBFY => GRID_CBFY
    use scale_comm, only: &
       COMM_vars8, &
       COMM_wait

    real(RP) :: coef_x, alpha_x1
    real(RP) :: coef_y, alpha_y1
    real(RP) :: ee1

    integer :: i, j
    !---------------------------------------------------------------------------

    ! check invalid fraction
    COPYTOPO_FRACX = max( min( COPYTOPO_FRACX, 1.0_RP ), EPS )
    COPYTOPO_FRACY = max( min( COPYTOPO_FRACY, 1.0_RP ), EPS )

    if ( COPYTOPO_taux <= 0.0_RP ) then ! invalid tau
       coef_x = 0.0_RP
    else
       coef_x = 1.0_RP / COPYTOPO_taux
    endif

    if ( COPYTOPO_tauy <= 0.0_RP ) then ! invalid tau
       coef_y = 0.0_RP
    else
       coef_y = 1.0_RP / COPYTOPO_tauy
    endif

    do j = 1, JA
    do i = 1, IA
       ee1 = CTRX(i)
       if ( ee1 <= 1.0_RP - COPYTOPO_FRACX ) then
          ee1 = 0.0_RP
       else
          ee1 = ( ee1 - 1.0_RP + COPYTOPO_FRACX ) / COPYTOPO_FRACX
       endif

       if ( COPYTOPO_LINEAR_H ) then
          alpha_x1 = coef_x * ee1
       else
          alpha_x1 = coef_x * ee1 * exp( -(1.0_RP-ee1) * COPYTOPO_EXP_H )
       end if

       ee1 = CTRY(j)
       if ( ee1 <= 1.0_RP - COPYTOPO_FRACY ) then
          ee1 = 0.0_RP
       else
          ee1 = ( ee1 - 1.0_RP + COPYTOPO_FRACY ) / COPYTOPO_FRACY
       endif

       if ( COPYTOPO_LINEAR_H ) then
          alpha_y1 = coef_y * ee1
       else
          alpha_y1 = coef_y * ee1 * exp( -(1.0_RP-ee1) * COPYTOPO_EXP_H )
       end if

       COPYTOPO_alpha(i,j) = max( alpha_x1, alpha_y1 )
    enddo
    enddo

    call COMM_vars8( COPYTOPO_alpha(:,:), 1 )
    call COMM_wait ( COPYTOPO_alpha(:,:), 1 )

    return
  end subroutine COPYTOPO_setalpha

  !-----------------------------------------------------------------------------
  !> Calc dumping coefficient alpha
  subroutine COPYTOPO_input_data( &
      topo_pd      ) ! (out)
    use scale_const, only: &
       D2R => CONST_D2R
    use scale_comm, only: &
       COMM_vars8, &
       COMM_wait
    use scale_grid_nest, only: &
       PARENT_IMAX,     &
       PARENT_JMAX,     &
       NEST_TILE_NUM_X, &
       NEST_TILE_NUM_Y, &
       NEST_TILE_ID,    &
       NEST_domain_shape
    use scale_interpolation_nest, only: &
       INTRPNEST_domain_compatibility,  &
       INTRPNEST_interp_fact_latlon,    &
       INTRPNEST_interp_2d
    use scale_grid_real, only: &
       LAT => REAL_LAT, &
       LON => REAL_LON
    implicit none

    real(RP), intent(out) :: topo_pd(:,:)

    real(RP)              :: dummy (1,1,1)
    real(RP), allocatable :: read2D(:,:)
    real(RP), allocatable :: lon_org (:,:)
    real(RP), allocatable :: lat_org (:,:)
    real(RP), allocatable :: topo_org(:,:)
    real(RP), allocatable :: hfact(:,:,:)
    integer,  allocatable :: igrd (:,:,:)
    integer,  allocatable :: jgrd (:,:,:)

    integer :: IALL, JALL  ! number of grids for whole domain
    integer :: PTI, PTJ    ! number of grids for a tile
    integer :: tilei, tilej
    integer :: rank
    integer :: i, j
    integer :: cxs, cxe, cys, cye  ! for child domain
    integer :: pxs, pxe, pys, pye  ! for parent domain
    !---------------------------------------------------------------------------

    PTI  = PARENT_IMAX(handle)
    PTJ  = PARENT_JMAX(handle)
    IALL = PTI * NEST_TILE_NUM_X
    JALL = PTJ * NEST_TILE_NUM_Y

    allocate( hfact   ( IA,   JA,  itp_nh ) )
    allocate( igrd    ( IA,   JA,  itp_nh ) )
    allocate( jgrd    ( IA,   JA,  itp_nh ) )
    allocate( lon_org ( IALL, JALL        ) )
    allocate( lat_org ( IALL, JALL        ) )
    allocate( topo_org( IALL, JALL        ) )

    do i = 1, size( NEST_TILE_ID(:) )
       ! read data from split files
       rank = NEST_TILE_ID(i)

       call NEST_domain_shape ( tilei, tilej,   & ! [out]
                                cxs,   cxe,     & ! [out]
                                cys,   cye,     & ! [out]
                                pxs,   pxe,     & ! [out]
                                pys,   pye,     & ! [out]
                                i               ) ! [in ]
       allocate( read2D  ( tilei,tilej ) )

       call FileRead( read2D(:,:), COPYTOPO_IN_BASENAME, "lon",  1, rank )
       lon_org (cxs:cxe,cys:cye) = read2D(pxs:pxe,pys:pye) * D2R
       call FileRead( read2D(:,:), COPYTOPO_IN_BASENAME, "lat",  1, rank )
       lat_org (cxs:cxe,cys:cye) = read2D(pxs:pxe,pys:pye) * D2R
       call FileRead( read2D(:,:), COPYTOPO_IN_BASENAME, "TOPO", 1, rank )
       topo_org(cxs:cxe,cys:cye) = read2D(pxs:pxe,pys:pye)

       deallocate( read2D )
    end do

    call INTRPNEST_domain_compatibility( lon_org(:,:), lat_org(:,:), dummy(:,:,:), &
                                         LON(:,:),     LAT(:,:),     dummy(:,:,:), &
                                         skip_z=.true.                             )
    !call check_domain_compatibility( lon_org(:,:), lat_org(:,:), dummy(:,:,:), &
    !                                 LON(:,:),     LAT(:,:),     dummy(:,:,:), &
    !                                 skip_z=.true.                             )

    call INTRPNEST_interp_fact_latlon( hfact  (:,:,:),   & ! (out)
                                       igrd   (:,:,:),   & ! (out)
                                       jgrd   (:,:,:),   & ! (out)
                                       LAT    (:,:),     & ! (in)
                                       LON    (:,:),     & ! (in)
                                       IA,   JA,         & ! (in)
                                       lat_org(:,:),     & ! (in)
                                       lon_org(:,:),     & ! (in)
                                       IALL, JALL        ) ! (in)

    call INTRPNEST_interp_2d( topo_pd (:,:),    &
                              topo_org(:,:),    &
                              hfact   (:,:,:),  &
                              igrd    (:,:,:),  &
                              jgrd    (:,:,:),  &
                              IA, JA            )

    call COMM_vars8( topo_pd(:,:), 1 )
    call COMM_wait ( topo_pd(:,:), 1 )

    return
  end subroutine COPYTOPO_input_data

  !-----------------------------------------------------------------------------
  !> Mixing TOPO data using parent and child domains
  subroutine COPYTOPO_mix_data( &
      topo_cd,  &  ! (inout)
      topo_pd    ) ! (in)
    implicit none
    real(RP), intent(inout) :: topo_cd(:,:)  ! topography of child domain (mine)
    real(RP), intent(in)    :: topo_pd(:,:)  ! topography of parent domain

    real(RP) :: frac
    integer  :: i, j
    !---------------------------------------------------------------------------

    if ( COPYTOPO_ENTIRE_REGION ) then
       topo_cd(:,:) = topo_pd(:,:)
    else
       do j = 1, JA
       do i = 1, IA
          frac = COPYTOPO_alpha(i,j)
          topo_cd(i,j) = topo_cd(i,j) * ( 1.0_RP - frac ) &
                       + topo_pd(i,j) * frac
       end do
       end do
    endif

    return
  end subroutine COPYTOPO_mix_data

end module mod_copytopo


text:
!> Implementation of the installation configuration.
!>
!> An install table can currently have the following fields
!>
!>```toml
!>library = bool
!>```
module fpm_manifest_install
  use fpm_error, only : error_t, fatal_error, syntax_error
  use fpm_toml, only : toml_table, toml_key, toml_stat, get_value
  implicit none
  private

  public :: install_config_t, new_install_config

  !> Configuration data for installation
  type :: install_config_t

    !> Install library with this project
    logical :: library

  contains

    !> Print information on this instance
    procedure :: info

  end type install_config_t

contains

  !> Create a new installation configuration from a TOML data structure
  subroutine new_install_config(self, table, error)

    !> Instance of the install configuration
    type(install_config_t), intent(out) :: self

    !> Instance of the TOML data structure
    type(toml_table), intent(inout) :: table

    !> Error handling
    type(error_t), allocatable, intent(out) :: error

    call check(table, error)
    if (allocated(error)) return

    call get_value(table, "library", self%library, .false.)

  end subroutine new_install_config

  !> Check local schema for allowed entries
  subroutine check(table, error)

    !> Instance of the TOML data structure
    type(toml_table), intent(inout) :: table

    !> Error handling
    type(error_t), allocatable, intent(out) :: error

    type(toml_key), allocatable :: list(:)
    integer :: ikey

    call table%get_keys(list)
    if (size(list) < 1) return

    do ikey = 1, size(list)
      select case(list(ikey)%key)
      case default
        call syntax_error(error, "Key "//list(ikey)%key//" is not allowed in install table")
        exit
      case("library")
        continue
      end select
    end do
    if (allocated(error)) return

  end subroutine check

  !> Write information on install configuration instance
  subroutine info(self, unit, verbosity)

    !> Instance of the build configuration
    class(install_config_t), intent(in) :: self

    !> Unit for IO
    integer, intent(in) :: unit

    !> Verbosity of the printout
    integer, intent(in), optional :: verbosity

    integer :: pr
    character(len=*), parameter :: fmt = '("#", 1x, a, t30, a)'

    if (present(verbosity)) then
      pr = verbosity
    else
      pr = 1
    end if

    if (pr < 1) return

    write(unit, fmt) "Install configuration"
    write(unit, fmt) " - library install", &
      & trim(merge("enabled ", "disabled", self%library))

  end subroutine info

end module fpm_manifest_install

max_stars_repo_name: LKedward/fpm-haskell

max_issues_repo_name: LKedward/fpm-haskell

max_forks_repo_name: LKedward/fpm-haskell