! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download

!> A selection of sorting routines. The simplest version sorts a given array
!> of values and returns a copy of the array with the items in ascending sorted
!> order.  This works on integer and real(r8) arrays.
!>
!> Another version returns an integer index array.  Accessing the original
!> list using these indices will traverse the array in sorted order.
!>
!> The final version of the sort routine is similar to the index sort
!> but a compare routine is supplied by the user, allowing this code
!> to sort a list of any kinds of items.
!>
!> @todo FIXME - some routines work in place and some return copies.
!> this affects the calling code in a very visible way.  can we make
!> these consistent?  or add an optional argument that says copy or no?
!> (hard to do because some are functions that return the copy; some
!> are subroutines that work in place.  would need to change the calling
!> code to be all subroutines and pass in either a single array for
!> "in place" or a src and dst array for copy.

module sort_mod

use     types_mod, only : r8, i8

implicit none
private

! what this should be:
!public :: sort, index_sort

! for now, let code indicate what sort they want.
public :: sort, index_sort, insertion_sort, index_insertion_sort
!public :: simple_sort, simple_index_sort

logical, save :: module_initialized = .false.

!> single interface to sort real or integer arrays

interface sort
   module procedure rsort
   module procedure isort
end interface sort

interface insertion_sort
   module procedure insertion_sort_real
   module procedure insertion_sort_int
end interface insertion_sort

! the simple sorts are here so we can time them,
! but in all cases they are slower and a bad choice.

interface simple_sort
   module procedure simple_sort_real
   module procedure simple_sort_int
end interface simple_sort

!> single interface to return indices in sorted order,
!> leaving the original array in place.

interface index_sort
   module procedure index_sort_real_int
   module procedure index_sort_real_i8
   module procedure index_sort_int_int
!  module procedure index_sort_i8_i8
   module procedure index_sort_user
end interface

interface index_insertion_sort
   module procedure index_insertion_sort_real
   module procedure index_insertion_sort_int
end interface index_insertion_sort

interface simple_index_sort
   module procedure simple_index_sort_real
   module procedure simple_index_sort_int
end interface simple_index_sort

contains

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

subroutine initialize_module()

if (module_initialized) return

module_initialized = .true.

end subroutine initialize_module

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

!> Uses a heap sort algorithm on real(r8) array x(), returns sorted array.
!> x() must be allocated or declared to be exactly the intended size;
!> all items in x() are sorted.
!>
!> @param x the (real) array to sort
!>

function rsort(x)

real(r8), intent(in) :: x(:)
real(r8)             :: rsort(size(x))

integer  :: num, level, indx, i, j
real(r8) :: l_val

! Get the array size
num = size(x)

! Copy to output
rsort = x

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
indx = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val = rsort(level)
   else
      l_val = rsort(indx)
      rsort(indx) = rsort(1)
      indx = indx - 1
      if(indx == 1) then
         rsort(1) = l_val
         return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= indx)
      if(j < indx) then
         if(rsort(j) < rsort(j + 1)) j = j + 1
      endif
      if(l_val < rsort(j)) then
         rsort(i) = rsort(j)
         i = j
         j = 2 * j
      else
         j = indx + 1
      endif
   end do

   rsort(i) = l_val

end do

end function rsort

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

!> Uses a heap sort algorithm on integer array x(), returns sorted array.
!> x() must be allocated or declared to be exactly the intended size;
!> all items in x() are sorted.
!>
!> @param x the (integer) array to sort
!>

function isort(x)

integer, intent(in) :: x(:)
integer             :: isort(size(x))

integer :: num, level, indx, i, j
integer :: l_val

! Get the array size
num = size(x)

! Copy to output
isort = x

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
indx = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val = isort(level)
   else
      l_val = isort(indx)
      isort(indx) = isort(1)
      indx = indx - 1
      if(indx == 1) then
         isort(1) = l_val
         return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= indx)
      if(j < indx) then
         if(isort(j) < isort(j + 1)) j = j + 1
      endif
      if(l_val < isort(j)) then
         isort(i) = isort(j)
         i = j
         j = 2 * j
      else
         j = indx + 1
      endif
   end do

   isort(i) = l_val

end do

end function isort

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

!> Uses a heap sort algorithm on x(), returns integer array of sorted indices.
!> x(num) array returned unchanged; index array doesn't need to be
!> initialized before this call.
!> usage: do i=1,num;  x(indx(i)) is next item in sorted order; enddo
!> This differs from index_sort_int_int only in the data type of the x() array.
!>
!> @param x the (real) array to sort
!> @param indx the array of integer indices to be used to traverse the input array in sorted order
!> @param num the length of x

subroutine index_sort_real_int(x, indx, num)

integer,  intent(in)  :: num
real(r8), intent(in)  :: x(num)
integer,  intent(out) :: indx(num)

integer  :: ind, i, j, l_val_indx, level
real(r8) :: l_val

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
   indx(i) = i
end do

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
ind = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val = x(indx(level))
      l_val_indx = indx(level)
   else
      l_val = x(indx(ind))
      l_val_indx = indx(ind)

      indx(ind) = indx(1)
      ind = ind - 1
      if(ind == 1) then
         indx(1) = l_val_indx
         return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= ind)
      if(j < ind) then
         if(x(indx(j)) < x(indx(j + 1))) j = j + 1
      endif
      if(l_val < x(indx(j))) then
         indx(i) = indx(j)
         i = j
         j = 2 * j
      else
         j = ind + 1
      endif
   end do

   indx(i) = l_val_indx

end do

end subroutine index_sort_real_int

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

!> Uses a heap sort algorithm on x(), returns long-integer array of sorted indices.
!> x(num) array returned unchanged; index array doesn't need to be
!> initialized before this call.
!> usage: do i=1,num;  x(indx(i)) is next item in sorted order; enddo
!> This differs from index_sort_real_int only in the data type of the integer array
!> and the data type of the length of the integer array.
!>
!> @param x the (real) array to sort
!> @param indx the array of long-integer indices to be used to traverse the input array in sorted order
!> @param num the length of x as a long-integer

subroutine index_sort_real_i8(x, indx, num)

integer(i8),  intent(in)  :: num
real(r8),     intent(in)  :: x(num)
integer(i8),  intent(out) :: indx(num)

integer(i8) :: ind, i, j, l_val_indx, level
real(r8) :: l_val

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
   indx(i) = i
end do

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
ind = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val = x(indx(level))
      l_val_indx = indx(level)
   else
      l_val = x(indx(ind))
      l_val_indx = indx(ind)

      indx(ind) = indx(1)
      ind = ind - 1
      if(ind == 1) then
         indx(1) = l_val_indx
         return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= ind)
      if(j < ind) then
         if(x(indx(j)) < x(indx(j + 1))) j = j + 1
      endif
      if(l_val < x(indx(j))) then
         indx(i) = indx(j)
         i = j
         j = 2 * j
      else
         j = ind + 1
      endif
   end do

   indx(i) = l_val_indx

end do

end subroutine index_sort_real_i8

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

!> Uses a heap sort algorithm on x(), returns integer array of sorted indices.
!> x(num) array returned unchanged; index array doesn't need to be
!> initialized before this call.
!> usage: do i=1,num;  x(indx(i)) is next item in sorted order; enddo
!> This differs from index_sort_real_int only in the data type of the x() array.
!>
!> @param x the (integer) array to sort
!> @param indx the array of integer indices to be used to traverse the input array in sorted order
!> @param num the length of x

subroutine index_sort_int_int(x, indx, num)

integer, intent(in)  :: num
integer, intent(in)  :: x(num)
integer, intent(out) :: indx(num)

integer :: ind, i, j, l_val_indx, level
integer :: l_val

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
   indx(i) = i
end do

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
ind = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val = x(indx(level))
      l_val_indx = indx(level)
   else
      l_val = x(indx(ind))
      l_val_indx = indx(ind)

      indx(ind) = indx(1)
      ind = ind - 1
      if(ind == 1) then
         indx(1) = l_val_indx
      return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= ind)
      if(j < ind) then
         if(x(indx(j)) < x(indx(j + 1))) j = j + 1
      endif
      if(l_val < x(indx(j))) then
         indx(i) = indx(j)
         i = j
         j = 2 * j
      else
         j = ind + 1
      endif
   end do

   indx(i) = l_val_indx

end do

end subroutine index_sort_int_int

!-------------------------------------------------------------------------
!> Uses a heap sort algorithm on x(), returns long-integer array of sorted indices.
!> x(num) array returned unchanged; index array doesn't need to be
!> initialized before this call.
!> usage: do i=1,num;  x(indx(i)) is next item in sorted order; enddo
!> This differs from index_sort_int_int only in the data type of the x() array.
!>
!> @param x the (long-integer) array to sort
!> @param indx the array of long-integer indices to be used to traverse the input array in sorted order
!> @param num the length of x

subroutine index_sort_i8_i8(x, indx, num)

implicit none

integer(i8), intent(in)  :: num
integer(i8), intent(in)  :: x(num)
integer(i8), intent(out) :: indx(num)

integer(i8) :: ind, i, j, l_val_indx, level
integer(i8) :: l_val

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
  indx(i) = i
end do

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
ind = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val = x(indx(level))
      l_val_indx = indx(level)
   else
      l_val = x(indx(ind))
      l_val_indx = indx(ind)
      indx(ind) = indx(1)
      ind = ind - 1
      if(ind == 1) then
         indx(1) = l_val_indx
         return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= ind)
      if(j < ind) then
         if(x(indx(j)) < x(indx(j + 1))) j = j + 1
      endif
      if(l_val < x(indx(j))) then
         indx(i) = indx(j)
         i = j
         j = 2 * j
      else
         j = ind + 1
      endif
   end do

   indx(i) = l_val_indx

end do

end subroutine index_sort_i8_i8

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

!> Uses a heap sort algorithm to compute a sorted array of indices.  The
!> actual items being sorted are opaque to this routine; the user supplies
!> a sorting function which must have access to the original items.  This
!> routine only manipulates the integer index values which are then returned.
!>
!> The third argument is a user-supplied function which takes 2 integers (a,b)
!> as arguments and return an integer code:
!>     -1 if mything(a) <  mything(b)
!>      0 if mything(a) == mything(b)
!>      1 if mything(a) >  mything(b)
!> The function only gets the index numbers; the user code should use them
!> to compare the items corresponding to those index numbers and return the
!> ordering between them.
!> This seems to work best when the comparefunc() is a public routine inside
!> a module so the interface specification is known to the compiler.
!>
!> usage: do i=1,num; mything(indx(i)) is next item in sorted order; enddo
!>
!> @param indx the array of integers
!> @param num the length of x
!> @param comparefunc the name of the comparison function
!>

subroutine index_sort_user(indx, num, comparefunc)

integer,  intent(in)  :: num
integer,  intent(out) :: indx(num)
interface
  integer function comparefunc(a, b)
    integer, intent(in) :: a, b
  end function comparefunc
end interface

integer :: ind, i, j, l_val_indx, level
integer :: compval

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
   indx(i) = i
enddo

! Only one element, just send it back
if(num <= 1) return

level = num / 2 + 1
ind = num

! Keep looping until finished
do
   ! Keep going down levels until bottom
   if(level > 1) then
      level = level - 1
      l_val_indx = indx(level)
   else
      l_val_indx = indx(ind)

      indx(ind) = indx(1)
      ind = ind - 1
      if(ind == 1) then
        indx(1) = l_val_indx
        return
      endif
   endif

   i = level
   j = 2 * level

   do while(j <= ind)
      if(j < ind) then
         compval = comparefunc(indx(j), indx(j+1))
         if(compval < 0) j = j + 1
      endif
      compval = comparefunc(l_val_indx, indx(j))
      if(compval < 0) then
         indx(i) = indx(j)
         i = j
         j = 2 * j
      else
         j = ind + 1
      endif

   enddo
   indx(i) = l_val_indx

enddo

end subroutine index_sort_user

!-------------------------------------------------------------------------
! Insertion sort can be fast for small arrays and for arrays that are
! already nearly sorted.  works in place.  real(r8) version.

subroutine insertion_sort_real(x)

real(r8), intent(inout) :: x(:)

integer :: i, j, num
real(r8) :: v

num = size(x)
i = 2
do while (i <= num)
  v = x(i)
  j = i - 1
  do while (j >= 1)
     if (x(j) <= v) exit
     x(j+1) = x(j)
     j = j - 1
  end do
  x(j + 1) = v
  i = i + 1
end do

end subroutine insertion_sort_real

!-------------------------------------------------------------------------
! Insertion sort can be fast for small arrays and for arrays that are
! already nearly sorted.  works in place.  integer version.

subroutine insertion_sort_int(x)

integer, intent(inout) :: x(:)

integer :: i, j, num, v

num = size(x)
i = 2
do while (i <= num)
  v = x(i)
  j = i - 1
  do while (j >= 1)
     if (x(j) <= v) exit
     x(j+1) = x(j)
     j = j - 1
  end do
  x(j + 1) = v
  i = i + 1
end do

end subroutine insertion_sort_int

!-------------------------------------------------------------------------
! index version for real(r8) data.

subroutine index_insertion_sort_real(x, indx, num)

integer,  intent(in) :: num
integer, intent(out) :: indx(num)
real(r8), intent(in) :: x(num)

integer :: i, j, v_indx
real(r8) :: v

do i=1, num
   indx(i) = i
enddo

i = 2
do while (i <= num)
   v = x(indx(i))
   v_indx = indx(i)
   j = i - 1
   do while (j >= 1)
      if (x(indx(j)) <= v) exit
      indx(j+1) = indx(j)
      j = j - 1
   end do
   indx(j + 1) = v_indx
   i = i + 1
end do

end subroutine index_insertion_sort_real

!-------------------------------------------------------------------------
! index version for integer data.

subroutine index_insertion_sort_int(x, indx, num)

integer, intent(in)  :: num
integer, intent(out) :: indx(num)
integer, intent(in)  :: x(num)

integer :: i, j, v_indx
integer :: v

do i=1, num
   indx(i) = i
enddo

i = 2
do while (i <= num)
  v = x(indx(i))
  v_indx = indx(i)
  j = i - 1
  do while (j >= 1)
     if (x(indx(j)) <= v) exit
     indx(j+1) = indx(j)
     j = j - 1
  end do
  indx(j + 1) = v_indx
  i = i + 1
end do

end subroutine index_insertion_sort_int

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

! a simple sort for real(r8) array data.  do not use - slower in all cases.
! (is this a variant of a selection sort?). works in place.

subroutine simple_sort_real(x)

real(r8), intent(inout) :: x(:)

real(r8) :: tmp
integer :: j, k, num

if ( .not. module_initialized ) call initialize_module

num = size(x)

! O(N^2) sort
do j = 1, num - 1
   do k = j + 1, num
      ! Exchange two elements if they're in the wrong order
      if(x(j) > x(k)) then
         tmp = x(k)
         x(k) = x(j)
         x(j) = tmp
      end if
   end do
end do

end subroutine simple_sort_real

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

! a simple sort for integer array data.  do not use - slower in all cases.
! (is this a variant of a selection sort?). works in place.

subroutine simple_sort_int(x)

integer, intent(inout) :: x(:)

integer :: tmp
integer :: j, k, num

if ( .not. module_initialized ) call initialize_module

num = size(x)

! O(N^2) sort
do j = 1, num - 1
   do k = j + 1, num
      ! Exchange two elements if they're in the wrong order
      if(x(j) > x(k)) then
         tmp = x(k)
         x(k) = x(j)
         x(j) = tmp
      end if
   end do
end do

end subroutine simple_sort_int

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

!  real(r8) index sort

subroutine simple_index_sort_real(x, indx, num)

integer,  intent(in)  :: num
integer,  intent(out) :: indx(num)
real(r8), intent(in)  :: x(num)

integer :: i, j, k, itmp

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
   indx(i) = i
end do

! O(N^2) sort
do j = 1, num
   do k = 1, num - 1
      ! Exchange two elements if they're in the wrong order
      if(x(indx(k)) > x(indx(k+1))) then
         itmp = indx(k)
         indx(k) = indx(k+1)
         indx(k+1) = itmp
      endif
   end do
end do

end subroutine simple_index_sort_real

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

!  integer  index sort

subroutine simple_index_sort_int(x, indx, num)

integer, intent(in)  :: num
integer, intent(out) :: indx(num)
integer, intent(in)  :: x(num)

integer :: i, j, k, itmp

if ( .not. module_initialized ) call initialize_module

! Initialize the index array to input order
do i = 1, num
   indx(i) = i
end do

! O(N^2) sort
do j = 1, num
   do k = 1, num - 1
      ! Exchange two elements if they're in the wrong order
      if(x(indx(k)) > x(indx(k+1))) then
         itmp = indx(k)
         indx(k) = indx(k+1)
         indx(k+1) = itmp
      endif
   end do
end do

end subroutine simple_index_sort_int

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

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

end module sort_mod


text:
      double precision pi,pisq,pisqo6
      parameter(pi=3.14159265358979d0,pisq=pi*pi,pisqo6=pisq/6d0)

      double precision zip,half,one,two,three,four,eight
      parameter(zip=0d0,half=0.5d0,one=1d0,two=2d0)
      parameter(three=3d0,four=4d0,eight=8d0)

      double complex im,impi,czip,chalf,cone,ctwo,c2ipi
      parameter(im=(0d0,1d0),impi=(0d0,3.14159265358979d0),
     . czip=(0d0,0d0),chalf=(0.5d0,0d0),cone=(1d0,0d0),ctwo=(2d0,0d0),
     . c2ipi=2d0*impi)

      DOUBLE COMPLEX nan
      parameter (nan = (1D123, 1D123))

ext: f


MODULE dynnxt
   !!=========================================================================
   !!                       ***  MODULE  dynnxt  ***
   !! Ocean dynamics: time stepping
   !!=========================================================================
   !! History :  OPA  !  1987-02  (P. Andrich, D. L Hostis)  Original code
   !!                 !  1990-10  (C. Levy, G. Madec)
   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
   !!            8.0  !  1997-02  (G. Madec & M. Imbard)  opa, release 8.0
   !!            8.2  !  1997-04  (A. Weaver)  Euler forward step
   !!             -   !  1997-06  (G. Madec)  lateral boudary cond., lbc routine
   !!    NEMO    1.0  !  2002-08  (G. Madec)  F90: Free form and module
   !!             -   !  2002-10  (C. Talandier, A-M. Treguier) Open boundary cond.
   !!            2.0  !  2005-11  (V. Garnier) Surface pressure gradient organization
   !!            2.3  !  2007-07  (D. Storkey) Calls to BDY routines.
   !!            3.2  !  2009-06  (G. Madec, R.Benshila)  re-introduce the vvl option
   !!            3.3  !  2010-09  (D. Storkey, E.O'Dea) Bug fix for BDY module
   !!            3.3  !  2011-03  (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL
   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
   !!            3.6  !  2014-04  (G. Madec) add the diagnostic of the time filter trends
   !!            3.7  !  2015-11  (J. Chanut) Free surface simplification
   !!-------------------------------------------------------------------------

   !!-------------------------------------------------------------------------
   !!   dyn_nxt       : obtain the next (after) horizontal velocity
   !!-------------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers
   USE dom_oce        ! ocean space and time domain
   USE sbc_oce        ! Surface boundary condition: ocean fields
   USE sbcrnf         ! river runoffs
   USE sbcisf         ! ice shelf
   USE phycst         ! physical constants
   USE dynadv         ! dynamics: vector invariant versus flux form
   USE dynspg_ts      ! surface pressure gradient: split-explicit scheme
   USE domvvl         ! variable volume
   USE bdy_oce   , ONLY: ln_bdy
   USE bdydta         ! ocean open boundary conditions
   USE bdydyn         ! ocean open boundary conditions
   USE bdyvol         ! ocean open boundary condition (bdy_vol routines)
   USE trd_oce        ! trends: ocean variables
   USE trddyn         ! trend manager: dynamics
   USE trdken         ! trend manager: kinetic energy
   !
   USE in_out_manager ! I/O manager
   USE iom            ! I/O manager library
   USE lbclnk         ! lateral boundary condition (or mpp link)
   USE lib_mpp        ! MPP library
   USE prtctl         ! Print control
   USE timing         ! Timing

   IMPLICIT NONE
   PRIVATE

   PUBLIC    dyn_nxt   ! routine called by step.F90

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: dynnxt.F90 10068 2018-08-28 14:09:04Z nicolasmartin $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS

   SUBROUTINE dyn_nxt ( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dyn_nxt  ***
      !!
      !! ** Purpose :   Finalize after horizontal velocity. Apply the boundary
      !!             condition on the after velocity, achieve the time stepping
      !!             by applying the Asselin filter on now fields and swapping
      !!             the fields.
      !!
      !! ** Method  : * Ensure after velocities transport matches time splitting
      !!             estimate (ln_dynspg_ts=T)
      !!
      !!              * Apply lateral boundary conditions on after velocity
      !!             at the local domain boundaries through lbc_lnk call,
      !!             at the one-way open boundaries (ln_bdy=T),
      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
      !!
      !!              * Apply the time filter applied and swap of the dynamics
      !!             arrays to start the next time step:
      !!                (ub,vb) = (un,vn) + atfp [ (ub,vb) + (ua,va) - 2 (un,vn) ]
      !!                (un,vn) = (ua,va).
      !!             Note that with flux form advection and non linear free surface,
      !!             the time filter is applied on thickness weighted velocity.
      !!             As a result, dyn_nxt MUST be called after tra_nxt.
      !!
      !! ** Action :   ub,vb   filtered before horizontal velocity of next time-step
      !!               un,vn   now horizontal velocity of next time-step
      !!----------------------------------------------------------------------
      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
      !
      INTEGER  ::   ji, jj, jk   ! dummy loop indices
      INTEGER  ::   ikt          ! local integers
      REAL(wp) ::   zue3a, zue3n, zue3b, zuf, zcoef    ! local scalars
      REAL(wp) ::   zve3a, zve3n, zve3b, zvf, z1_2dt   !   -      -
      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zue, zve
      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ze3u_f, ze3v_f, zua, zva
      !!----------------------------------------------------------------------
      !
      IF( ln_timing    )   CALL timing_start('dyn_nxt')
      IF( ln_dynspg_ts )   ALLOCATE( zue(jpi,jpj)     , zve(jpi,jpj)     )
      IF( l_trddyn     )   ALLOCATE( zua(jpi,jpj,jpk) , zva(jpi,jpj,jpk) )
      !
      IF( kt == nit000 ) THEN
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dyn_nxt : time stepping'
         IF(lwp) WRITE(numout,*) '~~~~~~~'
      ENDIF

      IF ( ln_dynspg_ts ) THEN
         ! Ensure below that barotropic velocities match time splitting estimate
         ! Compute actual transport and replace it with ts estimate at "after" time step
         zue(:,:) = e3u_a(:,:,1) * ua(:,:,1) * umask(:,:,1)
         zve(:,:) = e3v_a(:,:,1) * va(:,:,1) * vmask(:,:,1)
         DO jk = 2, jpkm1
            zue(:,:) = zue(:,:) + e3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)
            zve(:,:) = zve(:,:) + e3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)
         END DO
         DO jk = 1, jpkm1
            ua(:,:,jk) = ( ua(:,:,jk) - zue(:,:) * r1_hu_a(:,:) + ua_b(:,:) ) * umask(:,:,jk)
            va(:,:,jk) = ( va(:,:,jk) - zve(:,:) * r1_hv_a(:,:) + va_b(:,:) ) * vmask(:,:,jk)
         END DO
         !
         IF( .NOT.ln_bt_fw ) THEN
            ! Remove advective velocity from "now velocities"
            ! prior to asselin filtering
            ! In the forward case, this is done below after asselin filtering
            ! so that asselin contribution is removed at the same time
            DO jk = 1, jpkm1
               un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:) + un_b(:,:) )*umask(:,:,jk)
               vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:) + vn_b(:,:) )*vmask(:,:,jk)
            END DO
         ENDIF
      ENDIF

      ! Update after velocity on domain lateral boundaries
      ! --------------------------------------------------

      !
      CALL lbc_lnk_multi( ua, 'U', -1., va, 'V', -1. )     !* local domain boundaries
      !
      !                                !* BDY open boundaries
      IF( ln_bdy .AND. ln_dynspg_exp )   CALL bdy_dyn( kt )
      IF( ln_bdy .AND. ln_dynspg_ts  )   CALL bdy_dyn( kt, dyn3d_only=.true. )

!!$   Do we need a call to bdy_vol here??
      !
      IF( l_trddyn ) THEN             ! prepare the atf trend computation + some diagnostics
         z1_2dt = 1._wp / (2. * rdt)        ! Euler or leap-frog time step
         IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1._wp / rdt
         !
         !                                  ! Kinetic energy and Conversion
         IF( ln_KE_trd  )   CALL trd_dyn( ua, va, jpdyn_ken, kt )
         !
         IF( ln_dyn_trd ) THEN              ! 3D output: total momentum trends
            zua(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) * z1_2dt
            zva(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) * z1_2dt
            CALL iom_put( "utrd_tot", zua )        ! total momentum trends, except the asselin time filter
            CALL iom_put( "vtrd_tot", zva )
         ENDIF
         !
         zua(:,:,:) = un(:,:,:)             ! save the now velocity before the asselin filter
         zva(:,:,:) = vn(:,:,:)             ! (caution: there will be a shift by 1 timestep in the
         !                                  !  computation of the asselin filter trends)
      ENDIF

      ! Time filter and swap of dynamics arrays
      ! ------------------------------------------
      IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
         DO jk = 1, jpkm1
            un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
            vn(:,:,jk) = va(:,:,jk)
         END DO
         IF( .NOT.ln_linssh ) THEN                          ! e3._b <-- e3._n
!!gm BUG ????    I don't understand why it is not : e3._n <-- e3._a
            DO jk = 1, jpkm1
!               e3t_b(:,:,jk) = e3t_n(:,:,jk)
!               e3u_b(:,:,jk) = e3u_n(:,:,jk)
!               e3v_b(:,:,jk) = e3v_n(:,:,jk)
               !
               e3t_n(:,:,jk) = e3t_a(:,:,jk)
               e3u_n(:,:,jk) = e3u_a(:,:,jk)
               e3v_n(:,:,jk) = e3v_a(:,:,jk)
            END DO
!!gm BUG end
         ENDIF
                                                            !

      ELSE                                             !* Leap-Frog : Asselin filter and swap
         !                                ! =============!
         IF( ln_linssh ) THEN             ! Fixed volume !
            !                             ! =============!
            DO jk = 1, jpkm1
               DO jj = 1, jpj
                  DO ji = 1, jpi
                     zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
                     zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                     !
                     ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
                     vb(ji,jj,jk) = zvf
                     un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
                     vn(ji,jj,jk) = va(ji,jj,jk)
                  END DO
               END DO
            END DO
            !                             ! ================!
         ELSE                             ! Variable volume !
            !                             ! ================!
            ! Before scale factor at t-points
            ! (used as a now filtered scale factor until the swap)
            ! ----------------------------------------------------
            DO jk = 1, jpkm1
               e3t_b(:,:,jk) = e3t_n(:,:,jk) + atfp * ( e3t_b(:,:,jk) - 2._wp * e3t_n(:,:,jk) + e3t_a(:,:,jk) )
            END DO
            ! Add volume filter correction: compatibility with tracer advection scheme
            ! => time filter + conservation correction (only at the first level)
            zcoef = atfp * rdt * r1_rau0

            e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1)

            IF ( ln_rnf ) THEN
               IF( ln_rnf_depth ) THEN
                  DO jk = 1, jpkm1 ! Deal with Rivers separetely, as can be through depth too
                     DO jj = 1, jpj
                        DO ji = 1, jpi
                           IF( jk <=  nk_rnf(ji,jj)  ) THEN
                               e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( - rnf_b(ji,jj) + rnf(ji,jj) ) &
                                      &          * ( e3t_n(ji,jj,jk) / h_rnf(ji,jj) ) * tmask(ji,jj,jk)
                           ENDIF
                        ENDDO
                     ENDDO
                  ENDDO
               ELSE
                  e3t_b(:,:,1) = e3t_b(:,:,1) - zcoef *  ( -rnf_b(:,:) + rnf(:,:))*tmask(:,:,1)
               ENDIF
            END IF

            IF ( ln_isf ) THEN   ! if ice shelf melting
               DO jk = 1, jpkm1 ! Deal with isf separetely, as can be through depth too
                  DO jj = 1, jpj
                     DO ji = 1, jpi
                        IF( misfkt(ji,jj) <= jk .AND. jk <= misfkb(ji,jj) ) THEN
                           e3t_b(ji,jj,jk) =   e3t_b(ji,jj,jk) - zcoef *  ( fwfisf_b(ji,jj) - fwfisf(ji,jj) ) &
                                &          * ( e3t_n(ji,jj,jk) * r1_hisf_tbl(ji,jj) ) * tmask(ji,jj,jk)
                        ENDIF
                     END DO
                  END DO
               END DO
            END IF
            !
            IF( ln_dynadv_vec ) THEN      ! Asselin filter applied on velocity
               ! Before filtered scale factor at (u/v)-points
               CALL dom_vvl_interpol( e3t_b(:,:,:), e3u_b(:,:,:), 'U' )
               CALL dom_vvl_interpol( e3t_b(:,:,:), e3v_b(:,:,:), 'V' )
               DO jk = 1, jpkm1
                  DO jj = 1, jpj
                     DO ji = 1, jpi
                        zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) )
                        zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) )
                        !
                        ub(ji,jj,jk) = zuf                      ! ub <-- filtered velocity
                        vb(ji,jj,jk) = zvf
                        un(ji,jj,jk) = ua(ji,jj,jk)             ! un <-- ua
                        vn(ji,jj,jk) = va(ji,jj,jk)
                     END DO
                  END DO
               END DO
               !
            ELSE                          ! Asselin filter applied on thickness weighted velocity
               !
               ALLOCATE( ze3u_f(jpi,jpj,jpk) , ze3v_f(jpi,jpj,jpk) )
               ! Before filtered scale factor at (u/v)-points stored in ze3u_f, ze3v_f
               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3u_f, 'U' )
               CALL dom_vvl_interpol( e3t_b(:,:,:), ze3v_f, 'V' )
               DO jk = 1, jpkm1
                  DO jj = 1, jpj
                     DO ji = 1, jpi
                        zue3a = e3u_a(ji,jj,jk) * ua(ji,jj,jk)
                        zve3a = e3v_a(ji,jj,jk) * va(ji,jj,jk)
                        zue3n = e3u_n(ji,jj,jk) * un(ji,jj,jk)
                        zve3n = e3v_n(ji,jj,jk) * vn(ji,jj,jk)
                        zue3b = e3u_b(ji,jj,jk) * ub(ji,jj,jk)
                        zve3b = e3v_b(ji,jj,jk) * vb(ji,jj,jk)
                        !
                        zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n  + zue3a ) ) / ze3u_f(ji,jj,jk)
                        zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n  + zve3a ) ) / ze3v_f(ji,jj,jk)
                        !
                        ub(ji,jj,jk) = zuf                     ! ub <-- filtered velocity
                        vb(ji,jj,jk) = zvf
                        un(ji,jj,jk) = ua(ji,jj,jk)            ! un <-- ua
                        vn(ji,jj,jk) = va(ji,jj,jk)
                     END DO
                  END DO
               END DO
               e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1)        ! e3u_b <-- filtered scale factor
               e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1)
               !
               DEALLOCATE( ze3u_f , ze3v_f )
            ENDIF
            !
         ENDIF
         !
         IF( ln_dynspg_ts .AND. ln_bt_fw ) THEN
            ! Revert "before" velocities to time split estimate
            ! Doing it here also means that asselin filter contribution is removed
            zue(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
            zve(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
            DO jk = 2, jpkm1
               zue(:,:) = zue(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
               zve(:,:) = zve(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
            END DO
            DO jk = 1, jpkm1
               ub(:,:,jk) = ub(:,:,jk) - (zue(:,:) * r1_hu_n(:,:) - un_b(:,:)) * umask(:,:,jk)
               vb(:,:,jk) = vb(:,:,jk) - (zve(:,:) * r1_hv_n(:,:) - vn_b(:,:)) * vmask(:,:,jk)
            END DO
         ENDIF
         !
      ENDIF ! neuler =/0
      !
      ! Set "now" and "before" barotropic velocities for next time step:
      ! JC: Would be more clever to swap variables than to make a full vertical
      ! integration
      !
      !
      IF(.NOT.ln_linssh ) THEN
         hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1)
         hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1)
         DO jk = 2, jpkm1
            hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk)
            hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk)
         END DO
         r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
         r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
      ENDIF
      !
      un_b(:,:) = e3u_a(:,:,1) * un(:,:,1) * umask(:,:,1)
      ub_b(:,:) = e3u_b(:,:,1) * ub(:,:,1) * umask(:,:,1)
      vn_b(:,:) = e3v_a(:,:,1) * vn(:,:,1) * vmask(:,:,1)
      vb_b(:,:) = e3v_b(:,:,1) * vb(:,:,1) * vmask(:,:,1)
      DO jk = 2, jpkm1
         un_b(:,:) = un_b(:,:) + e3u_a(:,:,jk) * un(:,:,jk) * umask(:,:,jk)
         ub_b(:,:) = ub_b(:,:) + e3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk)
         vn_b(:,:) = vn_b(:,:) + e3v_a(:,:,jk) * vn(:,:,jk) * vmask(:,:,jk)
         vb_b(:,:) = vb_b(:,:) + e3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk)
      END DO
      un_b(:,:) = un_b(:,:) * r1_hu_a(:,:)
      vn_b(:,:) = vn_b(:,:) * r1_hv_a(:,:)
      ub_b(:,:) = ub_b(:,:) * r1_hu_b(:,:)
      vb_b(:,:) = vb_b(:,:) * r1_hv_b(:,:)
      !
      IF( .NOT.ln_dynspg_ts ) THEN        ! output the barotropic currents
         CALL iom_put(  "ubar", un_b(:,:) )
         CALL iom_put(  "vbar", vn_b(:,:) )
      ENDIF
      IF( l_trddyn ) THEN                ! 3D output: asselin filter trends on momentum
         zua(:,:,:) = ( ub(:,:,:) - zua(:,:,:) ) * z1_2dt
         zva(:,:,:) = ( vb(:,:,:) - zva(:,:,:) ) * z1_2dt
         CALL trd_dyn( zua, zva, jpdyn_atf, kt )
      ENDIF
      !
      IF(ln_ctl)   CALL prt_ctl( tab3d_1=un, clinfo1=' nxt  - Un: ', mask1=umask,   &
         &                       tab3d_2=vn, clinfo2=' Vn: '       , mask2=vmask )
      !
      IF( ln_dynspg_ts )   DEALLOCATE( zue, zve )
      IF( l_trddyn     )   DEALLOCATE( zua, zva )
      IF( ln_timing    )   CALL timing_stop('dyn_nxt')
      !
   END SUBROUTINE dyn_nxt

   !!=========================================================================
END MODULE dynnxt


text:
! This file is part of mctc-lib.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!> File type support
module mctc_io_filetype
   implicit none
   private

   public :: filetype, get_filetype

   !> Possible file types
   type :: enum_filetype

      !> Unknown file type
      integer :: unknown = 0

      !> xyz-format
      integer :: xyz = 1

      !> Turbomole coordinate format
      integer :: tmol = 2

      !> mol-format
      integer :: molfile = 3

      !> Vasp coordinate input
      integer :: vasp = 4

      !> Protein database format
      integer :: pdb = 5

      !> Structure data format
      integer :: sdf = 6

      !> GenFormat of DFTB+
      integer :: gen = 7

      !> Gaussian external format
      integer :: gaussian = 8

   end type enum_filetype

   !> File type enumerator
   type(enum_filetype), parameter :: filetype = enum_filetype()

contains

elemental function get_filetype(file) result(ftype)

   !> Name of the file
   character(len=*), intent(in) :: file

   !> File type from extension
   integer :: ftype

   integer :: iext, isep

   ftype = filetype%unknown
   iext = index(file, '.', back=.true.)
   isep = scan(file, '\/', back=.true.)

   if (iext > isep .and. iext > 0) then
      select case(to_lower(file(iext+1:)))
      case('coord', 'tmol')
         ftype = filetype%tmol
      case('xyz', 'log')
         ftype = filetype%xyz
      case('mol')
         ftype = filetype%molfile
      case('sdf')
         ftype = filetype%sdf
      case('poscar', 'contcar', 'vasp')
         ftype = filetype%vasp
      case('pdb')
         ftype = filetype%pdb
      case('gen')
         ftype = filetype%gen
      case('ein')
         ftype = filetype%gaussian
      end select
      if (ftype /= filetype%unknown) return
   else
      iext = len(file) + 1
   end if

   if (iext > isep) then
      select case(to_lower(file(isep+1:iext-1)))
      case('coord')
         ftype = filetype%tmol
      case('poscar', 'contcar')
         ftype = filetype%vasp
      end select
   end if

end function get_filetype

!> Convert input string to lowercase
elemental function to_lower(str) result(lcstr)

   !> Input string
   character(len=*), intent(in) :: str

   !> Lowercase version of string
   character(len=len(str)):: lcstr

   integer :: ilen, iquote, i, iav, iqc
   integer, parameter :: offset = iachar('A') - iachar('a')

   ilen = len(str)
   iquote = 0
   lcstr = str

   do i = 1, ilen
      iav = iachar(str(i:i))
      if (iquote == 0 .and. (iav == 34 .or.iav == 39)) then
         iquote = 1
         iqc = iav
         cycle
      end if
      if (iquote == 1 .and. iav==iqc) then
         iquote=0
         cycle
      end if
      if (iquote == 1) cycle
      if (iav >= iachar('A') .and. iav <= iachar('Z')) then
         lcstr(i:i) = achar(iav - offset)
      else
         lcstr(i:i) = str(i:i)
      end if
   end do

end function to_lower

end module mctc_io_filetype

max_stars_repo_name: loriab/mctc-lib

max_issues_repo_name: loriab/mctc-lib

max_forks_repo_name: loriab/mctc-lib


text:
*DECK DGAMI
      DOUBLE PRECISION FUNCTION DGAMI (A, X)
C***BEGIN PROLOGUE  DGAMI
C***PURPOSE  Evaluate the incomplete Gamma function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      DOUBLE PRECISION (GAMI-S, DGAMI-D)
C***KEYWORDS  FNLIB, INCOMPLETE GAMMA FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Evaluate the incomplete gamma function defined by
C
C DGAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) .
C
C DGAMI is evaluated for positive values of A and non-negative values
C of X.  A slight deterioration of 2 or 3 digits accuracy will occur
C when DGAMI is very large or very small, because logarithmic variables
C are used.  The function and both arguments are double precision.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  DGAMIT, DLNGAM, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DGAMI
      DOUBLE PRECISION A, X, FACTOR, DLNGAM, DGAMIT
C***FIRST EXECUTABLE STATEMENT  DGAMI
      IF (A .LE. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI',
     +   'A MUST BE GT ZERO', 1, 2)
      IF (X .LT. 0.D0) CALL XERMSG ('SLATEC', 'DGAMI',
     +   'X MUST BE GE ZERO', 2, 2)
C
      DGAMI = 0.D0
      IF (X.EQ.0.0D0) RETURN
C
C THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW.
      FACTOR = EXP (DLNGAM(A) + A*LOG(X))
C
      DGAMI = FACTOR * DGAMIT (A, X)
C
      RETURN
      END

max_stars_repo_name: floryst/ITK

max_issues_repo_name: floryst/ITK

max_forks_repo_name: floryst/ITK


text:
ibex_pkg.sv
ibex_alu.sv
ibex_compressed_decoder.sv
ibex_controller.sv
ibex_counters.sv
ibex_cs_registers.sv
ibex_decoder.sv
ibex_ex_block.sv
ibex_id_stage.sv
ibex_if_stage.sv
ibex_load_store_unit.sv
ibex_multdiv_slow.sv
ibex_multdiv_fast.sv
ibex_prefetch_buffer.sv
ibex_fetch_fifo.sv
ibex_register_file_ff.sv
ibex_wb_stage.sv
ibex_core.sv


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE MPINIT()
      COMMON /MAP/MODEL(4,4),VIEW(4,4),PROJ(4,4),CPNEAR,CPFAR,VWPORT(4),
     *       MVP(4,4),VP(4,4),CPLINE(2,2,10),CPPLAN(2,3,10),PEYE(3),
     *       PLOOK(3),ETWIST,NCPLIN,NCPLAN,TMAT1(4,4),TMAT2(4,4),
     *       TMAT3(4,4),TVEC1(4),TVEC2(4),TVEC3(4),TVEC4(4),TARR1(32),
     *       TARR2(32),TARR3(32),TARR4(32),TARR5(32),TARR6(32),
     *       TARR7(32),TARR8(32)
      REAL MODEL,MVP
      REAL DEVCAP(23)
      REAL DEFOUT(7)
      COMMON /STATUS/DEVCAP,DEFOUT
      REAL DEVP(5)
      COMMON /DEVICE/DEVP
      REAL COLP(3)
      REAL PALETT(3,16)
      COMMON /COLOR/COLP,PALETT
      REAL TEXTP(40)
      COMMON /TEXT/TEXTP
      REAL VECTP(5)
      REAL XCUR
      REAL YCUR
      COMMON /VECTRC/VECTP,XCUR,YCUR
      INTEGER IDEX(200,2)
      INTEGER NVECT(200,2)
      REAL XSIZE(200,2)
      REAL YSIZE(200,2)
      REAL X0(2300,2)
      REAL Y0(2300,2)
      REAL X1(2300,2)
      REAL Y1(2300,2)
      COMMON /FONT/IDEX,NVECT,XSIZE,YSIZE,X0,Y0,X1,Y1
      REAL GRAPHP(100)
      COMMON /GRAPH/GRAPHP
      COMMON /MAPPAR/MAPP(11)
      REAL MAPP
      COMMON /STORAG/MEMORY(1000)

      CALL MXIDEN(4,MODEL)
      CALL MXIDEN(4,VIEW)
      CALL MXIDEN(4,PROJ)
      CALL MPVIEW(0.,DEVP(4),0.,DEVP(5))
      NCPLIN = 0
      NCPLAN = 0
      CALL MXIDEN(4,MVP)
      CALL MXIDEN(4,VP)
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
C======================================================================
C  SU_NFACTO, Subroutine
C
C  Determines N deficit in Sunflower
C----------------------------------------------------------------------
C  Revision history
C  xx/xx/19xx     Written
C  02/08/1993 PWW Header revision and minor changes
C  03/29/2001 WDB Converted to modular form
C  12/01/2001 WDB Further modular conversion
!  02/25/2005 CHP Check for NFAC < 0.
C----------------------------------------------------------------------
C
C  Called : SU_GROSUB
C
C  Calls  : None
C----------------------------------------------------------------------

      SUBROUTINE SU_NFACTO(DYNAMIC,               !Control
     %    TANC,TCNP,TMNC,xlanc,xlcnp,xlmnc,       !Inputs
     %    AGEFAC,NDEF3,NFAC,NSTRES)               !Outputs

      USE ModuleDefs
      IMPLICIT NONE
      SAVE
C----------------------------------------------------------------------
C                     Variable Declaration
C----------------------------------------------------------------------
      REAL        AGEFAC
      REAL        NDEF3
      REAL        NFAC
      REAL        NSTRES
      REAL        TANC
      REAL        TCNP
      REAL        TMNC
      INTEGER     DYNAMIC
      REAL        xlanc,xlcnp,xlmnc
C     ----------------------------------------------------------------

      IF(DYNAMIC.EQ.RUNINIT) THEN
          NSTRES = 1.0
          AGEFAC = 1.0
          NDEF3 = 1.0
          NFAC = 1.0

      ELSEIF(DYNAMIC.EQ.SEASINIT) THEN
          NSTRES = 1.0
          AGEFAC = 1.0
          NDEF3 = 1.0
          NFAC = 1.0

      ELSEIF(DYNAMIC.EQ.INTEGR) THEN

! If the actual nitrogen content in the above ground biomass (TANC)
! decreases below a critical level (TCNP), then compute N stress
! based on proportion nitrogen below critical level (TCNP-TANC) and
! total nitrogen between the critical level and minimum
! level (TCNP-TMNC).
        NFAC   = 1.0 - (TCNP-TANC)/(TCNP-TMNC)
c       write(*,*)'nfacto 62',tcnp,tanc,tmnc,NFAC
          NFAC   = AMIN1 (NFAC,1.0)
          NFAC   = AMAX1 (NFAC,0.001)

          !Compute nitrogen stress factor for reducing leaf expansion

          AGEFAC = NFAC
          AGEFAC = AMIN1 (AGEFAC,1.0)

          !Make grain per plant calculation less sensitive to N stress
        NDEF3  = 1.0
        IF (NFAC .LT. 0.8) THEN
          NDEF3 = AMIN1 (0.2 + NFAC,1.0)
        ENDIF

          !Compute nitrogen stress factor for PHOTOSYNTHESIS

        NSTRES =  1.0 - (XLCNP-XLANC)/(XLCNP-XLMNC)
        IF (NSTRES .GT. 0.5)  THEN
          NSTRES = NSTRES * 0.4  + 0.6
        ELSE
          IF (NSTRES .GT. 0.0) THEN
            NSTRES = NSTRES * 1.2  + 0.2
          ELSE
            NSTRES = 0.2
          ENDIF
        ENDIF
        NSTRES = AMIN1 (NSTRES,1.0)
      ENDIF

      RETURN

      END SUBROUTINE SU_NFACTO

C----------------------------------------------------------------------
C                         DEFINITIONS
C----------------------------------------------------------------------
! AGEFAC      !Nitrogen stress factor affecting expansion (0-1)
! ISTAGE      !Plant growth staqe
! NDEF3       !Nitrogen stress factor affecting grains per plant (0-1)
! NFAC        !Nitrogen stress factor based on actual and critical
!              nitrogen content in vegetative tisue
! NSTRES      !Nitrogen stress factor affecting growth (0-1)
! TANC        !Nitrogen content in above ground biomass, decimal
! TCNP        !Critical nitrogen concentration in tops, g N/g dry wt.
! TMNC        !Plant top minimum N concentration g N/g dry matter
! DYNAMIC     !Modular control variable
! RUNINIT     !DYNAMIC control variable
! SEASINIT    !DYNAMIC control variable
! RATE        !DYNAMIC control variable
! INTEGR      !DYNAMIC control variable
! OUTPUT      !DYNAMIC control variable
! SEASEND     !DYNAMIC control variable
C----------------------------------------------------------------------

ext: for

max_stars_repo_name: xiangyi-wang/dssat-csm-os

max_issues_repo_name: xiangyi-wang/dssat-csm-os

max_forks_repo_name: xiangyi-wang/dssat-csm-os


text:
!------------------------------------------------------------------------------------
! File:   HMC_Module_Phys_Convolution_Apps_HydraulicStructure.f90
!
! Author(s):    Fabio Delogu, Francesco Silvestro, Simone Gabellani
! Date:         20190410
!
! Convolution Apps HydraulicStructure subroutine(s) for HMC model
!------------------------------------------------------------------------------------

!------------------------------------------------------------------------------------
! Module Header
module HMC_Module_Phys_Convolution_Apps_HydraulicStructure

    !------------------------------------------------------------------------------------
    ! External module(s)
    use HMC_Module_Namelist,                only: oHMC_Namelist
    use HMC_Module_Vars_Loader,             only: oHMC_Vars

    use HMC_Module_Tools_Debug

    use HMC_Module_Phys_HydraulicStructure, only: HMC_Phys_Lake_Tank, HMC_Phys_Dam_Discharge

    ! Implicit none for all subroutines in this module
    implicit none
    !------------------------------------------------------------------------------------------

contains

    !------------------------------------------------------------------------------------------
    ! Subroutine to compute hydraulic structure effect(s)
    subroutine HMC_Phys_Convolution_Apps_HydraulicStructure(iID, iRows, iCols, &
                                                            dDtDataForcing, iNDam, iNLake)

        !------------------------------------------------------------------------------------------
        ! Variable(s)
        integer(kind = 4)           :: iID, iRows, iCols
        integer(kind = 4)           :: iNDam, iNLake

        real(kind = 4)              :: dDtDataForcing

        real(kind = 4), dimension (iRows, iCols)    :: a2dVarFlowDeep, a2dVarHydro

        real(kind = 4), dimension (iNLake)          :: a1dVarVLake, a1dVarQoutLake

        real(kind = 4), dimension (iNDam)           :: a1dVarVDam, a1dVarHDam, a1dVarLDam
        real(kind = 4), dimension (iNDam)           :: a1dVarCoeffDam, a1dVarQoutDam
        !------------------------------------------------------------------------------------------

        !------------------------------------------------------------------------------------------
        ! Get variable(s)
        a2dVarFlowDeep = oHMC_Vars(iID)%a2dFlowDeep
        a2dVarHydro = oHMC_Vars(iID)%a2dHydro

        a1dVarVLake = oHMC_Vars(iID)%a1dVLake
        a1dVarQoutLake = oHMC_Vars(iID)%a1dQoutLake

        a1dVarVDam = oHMC_Vars(iID)%a1dVDam
        a1dVarHDam = oHMC_Vars(iID)%a1dHDam
        a1dVarLDam = oHMC_Vars(iID)%a1dLDam
        a1dVarCoeffDam = oHMC_Vars(iID)%a1dCoeffDam

        ! Info start
        if (iDEBUG.gt.0) then
            call mprintf(.true., iINFO_Extra, ' Phys :: Convolution :: HydraulicStructure ... ' )
        endif
        !------------------------------------------------------------------------------------------

        !------------------------------------------------------------------------------------------
        ! Method to compute lake tank
        call HMC_Phys_Lake_Tank(iID, iRows, iCols, dDtDataForcing, iNLake, &
                                a1dVarVLake, a1dVarQoutLake, &
                                a2dVarHydro, a2dVarFlowDeep)
        !------------------------------------------------------------------------------------------

        !------------------------------------------------------------------------------------------
        ! Method to compute dam discharge --> Commentata per cambio calcolo di svaso delle dighe
        !call HMC_Phys_Dam_Discharge(iID, iNDam, dDtDataForcing, &
        !                                a1dVarVDam, a1dVarHDam, a1dVarLDam, a1dVarCoeffDam)
        !------------------------------------------------------------------------------------------

        !------------------------------------------------------------------------------------------
        ! Updating field(s) in global declaration
        oHMC_Vars(iID)%a2dFlowDeep = a2dVarFlowDeep
        oHMC_Vars(iID)%a2dHydro = a2dVarHydro

        oHMC_Vars(iID)%a1dVLake = a1dVarVLake
        oHMC_Vars(iID)%a1dQoutLake = a1dVarQoutLake

        oHMC_Vars(iID)%a1dVDam = a1dVarVDam
        oHMC_Vars(iID)%a1dHDam = a1dVarHDam
        oHMC_Vars(iID)%a1dLDam = a1dVarLDam
        oHMC_Vars(iID)%a1dCoeffDam = a1dVarCoeffDam

        ! Info end
        if (iDEBUG.gt.0) then
            call mprintf(.true., iINFO_Extra, ' Phys :: Convolution :: HydraulicStructure ... OK' )
        endif
        !------------------------------------------------------------------------------------------

    end subroutine HMC_Phys_Convolution_Apps_HydraulicStructure
    !------------------------------------------------------------------------------------------

end module HMC_Module_Phys_Convolution_Apps_HydraulicStructure
!------------------------------------------------------------------------------------------

max_stars_repo_name: c-hydro/hmc-dev

max_issues_repo_name: c-hydro/hmc-dev

max_forks_repo_name: c-hydro/hmc-dev


text:
C    @(#)comdrp.f	20.3 2/13/96
      subroutine comdrp (annote, iter, num1, deldrp, num2, delpku,
     1                   status, getpqv, numitr)

C              This subroutine computes the generation reallocation.
C
C     Parameters: ANNOTE = Character string identifying calls.
C                 ITER   = Current solution iteration number.
C                 NUM1   = Count of busses which dropped generation
C                          this iteration.
C        (Input)  DELDRP = MW value of unallocated generation
C                          which must be dropped.
C        (Output) DELDRP = MW value of unallocated generation
C                          remaining to be dropped.
C                 NUM2   = Count of busses which picked up generation
C                          this iteration.
C                 DELPKU = MW value of generation picked up this
C                          iteration.
C                 STATUS = 0/1 = no errors/errors encountered.
C
C                 GETPQV = Subroutine name for computing P, Q, V.
C
C                 NUMITR = Cumulative iteration count.
C

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/alpha2.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/comdrx.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/gendrp.inc'
      include 'ipfinc/intbus.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/prt.inc'
c
c***kln Single to double precision.
c
      double precision pgen, pnew, pmax, pmin, pnewmw, dropmw
      double precision totpku, pmaxmw, pminmw, poldmw, totdrp
      double precision delpct, pgenmw, pickup

      integer status
      character btyp * 1, commnt * 37, annote * (*)
      external getpqv

      num1 = 0
      num2 = 0

      if (idswb .ne. 0) then

         call forbtm
         write (outbuf,10) annote, iter
   10    format (t10, a, i3, ' Status of machines with ',
     &           'generation dropped or picked up ')
         call shdlod(1)

         write (outbuf,20)
   20    format('0 Generator', t19, 'Type', t25, 'Zone',
     1      t31, '   ----------------- Generation -----------------',
     2      t95,'Comments')
         call shdlod(2)

         write (outbuf,30)
   30    format(t34,'Minimum    Maximum    Initial    Final    Change')
         call shdlod(3)

         write (outbuf,40)
   40    format(t34,' (MW)       (MW)       (MW)      (MW)       (MW)')
         call shdlod(4)

         outbuf = ' '
         call shdlod(5)
         call fortop
         write (outbuf, 50) drptot, 'Initially dropped generation'
   50    format (t71, f10.1, t95, a)
         call prtout (1)
      endif
      itrtot = itrtot + 1
      numitr = itrtot
      status = 0
C         Compute the total generation to be dropped. The varying nature
C         of system and area slack busses require recomputation of this
C         quantity.
C         TOTDRP = total generation dropped in this case.
C         DELDRP = generation dropped (deficit) this iteration.
      if (oldrop .eq. -9.0e10) then
         oldrop = -drptot
      endif
      totdrp = oldrop
      deldrp = 0.0
      do 110 i = 1, numdrp
         nb = gndpno(i)
         kt = inp2opt(nb)
         pmaxmw = gndpmx(i)
         pminmw = gndpmn(i)
         pmin = pminmw / bmva
         pmax = pmaxmw / bmva
         if (annote .eq. 'DC Iteration') then
            if (gndpty(i) .ne. 0 .and. numitr .gt. 1) then
C                       Update P- and Q-injections for area slack buses.
               call getpqv (kt,pk,dpk,qk,dqk,vk)
               pnetu(kt) = pk
            endif
            pgen = pnetu(kt) + ploadu(kt)
         else
            if (gndpty(i) .ne. 0) then
C                        Update P- and Q-injections for area slack buses.
               call getpqv (kt,pk,dpk,qk,dqk,vk)
               pnetu(kt) = pk
            endif
            pgen = pnetu(kt) + ploadu(kt)
            pnew = pgen - ddim (pgen, pmax) + ddim (pmin, pgen)
            pnetu(kt) = pnew - ploadu(kt)
         endif
         pnew = pgen - ddim (pgen, pmax) + ddim (pmin, pgen)
         pnetu(kt) = pnew - ploadu(kt)
         pnewmw = pnew * bmva
         pgenmw = pgen * bmva
         dropmw = pnewmw - pgenmw
         totdrp = totdrp + dropmw
         deldrp = deldrp + dropmw
C
C        Update PGEN quantities on BUSDTA for non-slack busses.
C
         if (gndpty(i) .eq. 0) then
            busdta(8,nb) = busdta(8,nb) + dropmw
         endif
         if (dropmw .lt. -0.5) then
            num1 = num1 + 1
            commnt = 'Generation dropped'
         else
            commnt = ' '
         endif
         if (idswb .ne. 0) then
            call typno (btyp, ntypu(kt))
            write (outbuf, 100) bus(nb), base(nb), btyp, zone(nb),
     1         pminmw, pmaxmw, pgenmw, pnewmw, dropmw, commnt
  100       format (t3, a8, f6.1, t21, a1, t27, a2, t31, 5f10.1,
     1         t95, a)
            if (gndpty(i) .ne. 0) outbuf(122:) = '(Slack bus)'
            call prtout (1)
         endif

  110 continue

      if (itrtot .eq. 1) then
         target = totdrp
      else
         target = deldrp
      endif
      delpku = 0.0
C
C     "TARGET" is the additional generation pickup desired for this
C     subroutine call.
C
      do 170 itx = 1, 10

C
C        First pass: determine the number and amount of eligible
C        generators for pickup (NUM2, TOTPMX).
C
         num2 = 0
         spinpu = 0.0
         totpmx = 0.0
C
C        NUM2   = total number of eligible pickup generators
C                 (P_i < P_i_max)
C
C        SPINPU = total spinning reserve (p.u.).
C        TOTPMX = total P_max on eligible pickup generators.
C
         do 130 i = 1, numgen
            nb = gennum(i)
            kt = inp2opt(nb)
            pmax = busdta(7,nb) / bmva
C                     Update P- and Q-injections for area slack buses.
               if (gentyp(i) .ne. 0) then
C                     Update P- and Q-injections for area slack buses.
                  call getpqv (kt,pk,dpk,qk,dqk,vk)
                  pnetu(kt) = pk
               endif
               pgen = pnetu(kt) + ploadu(kt)
            if (itx .eq. 1) then
               pold(i) = pgen
C                    Compensate DELPKU for slack bus excursions since
C                    previous iteration.
               if (gentyp(i) .ne. 0) then
                  j = gentyp(i)
                  if (itrtot .eq. 1) slkgen(j) = genpol(i) / bmva
                  delpku = delpku + (pgen - slkgen(j)) * bmva
                  slkgen(j) = pgen
               endif
            endif
            if (pgen + 0.005 .lt. pmax .and.
     1          pmax .gt. 0.0 .and.
     2          gentyp(i) .eq. 0) then
               spinpu = spinpu + ddim (pmax, pgen)
               totpmx = totpmx + pmax
               num2 = num2 + 1
            endif
  130    continue
C
C        Compute percentage pickup PCTPKU to balance DELDRP with
C        TARGET. Note: TARGET <= 0 means pickup is => 0.
C
         if (totpmx .eq. 0) then
            delpct = 0.0
            go to 172
         else
            delpct = (-delpku - target) / bmva / totpmx
         endif
         pctpku = pctpku + delpct
C
C        Now apply this percentage pickup PCTPKU until TARGET is
C        reached.
C
         totpku = 0.0

         do 160 i = 1, numgen
            nb = gennum(i)
            kt = inp2opt(nb)
            pmax = busdta(7,nb) / bmva
            pgen = pnetu(kt) + ploadu(kt)
            pnew = dmin1 (genpol(i)/bmva + pctpku * pmax, pmax)
C                   Update PGEN quantities on BUSDTA for non-slack busses.
            if (gentyp(i) .eq. 0) then
               busdta(8,nb) = busdta(8,nb) + (pnew - pgen) * bmva
                  pnetu(kt) = pnew - ploadu(kt)
               pnetu(kt) = pnew - ploadu(kt)
               totpku = totpku + pnew * bmva - genpol(i)
               delpku = delpku + (pnew - pgen) * bmva
            endif

  160    continue

         if (gensum_flag .eq. 1) then
           if (itx .eq. 1) call space (1)
           write (outbuf, 162) itx, num1, target, num2, delpku,
     1        totpmx * bmva, pctpku * 100.0
  162      format ('  Iteration ',i3,  ' Total dropped ', i3, ' (',
     1        f10.1, ') Total pickup ', i3, ' (', f10.1, ') P_max (',
     2        f10.1, ') Pickup (', f8.2, ') % ')
           call prtout (1)
         endif
C
C        Allocate PCTPKU until DELDRP < DRPTOL
C
         if (abs(delpku + target) .le. drptol) go to 176

  170 continue

  172 write (errbuf(1), 174 ) 10, spinpu * bmva, target, delpku
  174 format('Insufficent spinning reserve after ',i2,
     1       ' iterations (',f10.1,'), amount dropped (', f10.1,
     2       ') amount pickup (',f10.1,')')
      call prterx ('W', 1)
      status = 1

  176 continue
C
C     Summarize new allocated generation.
C
      totpku = 0.0
      delpku = 0.0
      num2 = 0

      do 230 i  = 1, numgen
         nb = gennum(i)
         kt = inp2opt(nb)
         pmaxmw = busdta(7,nb)
         pminmw = 0.0
            pgen = pnetu(kt) + ploadu(kt)

         poldmw = genpol(i)
         pgenmw = pgen * bmva
         pickup = pgenmw - poldmw
         totpku = totpku + pickup
         j = gentyp(i)
         if (j .gt. 0) then
            delpku = delpku + pgenmw - slkgen(j) * bmva
         else
            delpku = delpku + pgenmw - pold(i) * bmva
         endif
         if (pgenmw - 0.5 .gt. poldmw) then
            num2 = num2 + 1
         endif

         if (idswb .ne. 0) then
            call typno (btyp, ntypu(kt))
            if (pmaxmw .gt. 0.0) then
               pct = 100.0 * pickup / pmaxmw
            else
               pct = 0.0
            endif
            write (outbuf, 220) bus(nb), base(nb), btyp, zone(nb),
     1         pminmw, pmaxmw, poldmw, pgenmw, pickup, pct
  220       format (t3, a8, f6.1, t21, a1, t27, a2, t31, 5f10.1,
     1         t95, 'Generation Pickup ',f6.2, ' %')
            if (gentyp(i) .ne. 0) outbuf(122:) = '(Slack bus)'
            call prtout (1)
         endif
  230 continue

      oldrop = totdrp + totpku

      if (idswb .ne. 0) then
         write (outbuf, 240) totdrp
  240    format ('0 Total dropped', t71, f10.1)
         call prtout (1)
         write (outbuf, 242) totpku
  242    format ('  Total pickup', t71, f10.1)
         call prtout (1)
      endif

      call space (1)

      do 260 i = 1, 5
         outbuf = ' '
         call shdlod(i)
  260 continue

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!   Copyright 2014 College of William and Mary
!
!   Licensed under the Apache License, Version 2.0 (the "License");
!   you may not use this file except in compliance with the License.
!   You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
!   Unless required by applicable law or agreed to in writing, software
!   distributed under the License is distributed on an "AS IS" BASIS,
!   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!   See the License for the specific language governing permissions and
!   limitations under the License.

! Generate [TEM,SAL]_nu.nc from gridded data (nc file)
! The section on nc read needs to be modified as appropriate- search for
! 'start11' and 'end11'
! Beware the order of vertical levels in the nc file!!!
! Assume elev=0 in interpolation of 3D variables
! The code will do all it can to extrapolate: below bottom/above surface.
! If a pt in hgrid.ll is outside the background nc grid, const. values will be filled.
!
! ifort -mcmodel=medium -CB -O2 -o gen_nudge_from_hycom.exe ../UtilLib/pt_in_poly_test.f90 ../UtilLib/compute_zcor.f90 gen_nudge_from_hycom.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff

!   Input:
!     (1) hgrid.gr3;
!     (2) hgrid.ll;
!     (3) vgrid.in (SCHISM R3000 and up);
!     (4) TEM_nudge.gr3: used to mark elem's that need nudging (assuming identical to SAL_nudge.gr3!)
!     (5) gen_nudge_from_nc.in:
!                     1st line: T,S values for pts outside bg grid in nc (make sure nudging zone is inside bg grid in nc)
!                     2nd line: time step in .nc in sec; output stride
!                     3rd line: # of nc files
!     (6) HYCOM files: TS_[1,2,..nfiles].nc (includes lon/lat; beware scaling etc)
!                      The extent the HYCOM files cover needs to be
!                      larger than the nudging region specified in TEM_nudge.gr3,
!                      and lon/lat coord monotonically increasing.
!   Output: [TEM,SAL]_nu.nc (reduced to within nudging zone only); include3.gr3 (nudging zone)
!   Debug outputs: fort.11 (fatal errors); fort.*

      program gen_hot
      use netcdf
      use compute_zcor
      use pt_in_poly_test

!      implicit none

!      integer, parameter :: debug=1
!      integer, parameter :: mnp=600000
!      integer, parameter :: mne=1200000
!      integer, parameter :: mnv=60
      integer, parameter :: nbyte=4
      real, parameter :: small1=1.e-2 !used to check area ratios

!     netcdf related variables
      integer :: sid, tid, ncids(100),one_dim ! Netcdf file IDs
      integer :: latvid, lonvid, zmvid, hvid ! positional variables
      integer :: xdid, xvid ! longitude index
      integer :: ydid, yvid ! latitude index
      integer :: ldid, lvid ! vertical level, 1 is top
      integer :: svid, tvid ! salt & temp variable IDs
      integer :: uvid, vvid ! vel variable IDs
      integer :: evid ! SSH variable IDs
      integer :: var1d_dims(1),var4d_dims(4)
      integer, dimension(nf90_max_var_dims) :: dids

!     Local variables for data
      real (kind = 4), allocatable :: xind(:), yind(:), lind(:)
!     Lat, lon, bathymetry
      real (kind = 4), allocatable :: lat(:,:), lon(:,:), hnc(:)
!     Vertical postion, salinity, and temperature
      real (kind = 4), allocatable :: zm(:,:,:),salt(:,:,:),temp(:,:,:), &
     &uvel(:,:,:),vvel(:,:,:),ssh(:,:)
      integer, allocatable :: kbp(:,:),ihope(:,:)
!     File names for netcdf files
      character(len=1024) :: ncfile1,ncfile2
!     Command line arguments
      character(len=1024) :: s, yr, md
      character(len=4) :: iyear_char
      character(len=1) :: char1,char2
      character(len=20) :: char3,char4
      character(len=2) :: hr_char(4)
!     external function for number of command line arguments
      integer :: iargc

      integer :: status ! netcdf local status variable
      integer :: ier ! allocate error return.
      integer :: ixlen, iylen, ilen ! sampled lengths in each coordinate direction
      integer :: ixlen1, iylen1,ixlen2, iylen2 !reduced indices for CORIE grid to speed up interpolation
      real, allocatable :: xl(:),yl(:),dp(:),tempout(:,:),saltout(:,:),ztot(:),sigma(:),arco(:,:)
      integer, allocatable :: i34(:),elnode(:,:),include2(:),ixy(:,:),imap(:)
      dimension wild(100),wild2(100,2)
      dimension nx(4,4,3),month_day(12)
      dimension ndays_mon(12)
      allocatable :: z(:,:),sigma_lcl(:,:),kbp2(:),iparen_of_dry(:,:)
      real*8 :: aa1(1),dtmp

!     First statement
!     Currently we assume rectangular grid in HYCOM
!     (interp_mode=0). interp_mode=1 uses generic UG search (splitting
!     quads) and is kept for more generic cases
      interp_mode=0

      ndays_mon=(/31,28,31,30,31,30,31,31,30,31,30,31/)  !# of days in each month for non-leap yr
!      hr_char=(/'03','09','15','21'/) !each day has 4 starting hours in ROMS

      open(10,file='gen_nudge_from_nc.in',status='old')
      read(10,*) tem_outside,sal_outside !T,S values for pts outside bg grid in nc or TEM_nudge.gr3
      read(10,*) dtout,nt_out !time step in .nc [sec], output stride
      !read(10,*) istart_year,istart_mon,istart_day
      read(10,*) nndays !# of nc files
      close(10)

!     Read in hgrid and vgrid
      open(16,file='hgrid.ll',status='old')
      open(15,file='TEM_nudge.gr3',status='old')
      open(14,file='hgrid.gr3',status='old') !only need depth info and connectivity
      open(19,file='vgrid.in',status='old')
      open(11,file='fort.11',status='replace')
      read(14,*)
      read(14,*)ne,np
      allocate(xl(np),yl(np),dp(np),i34(ne),elnode(4,ne),include2(np),ixy(np,2),arco(4,np),imap(np))
      read(16,*); read(16,*)
      read(15,*); read(15,*)
      do i=1,np
        read(14,*)j,xtmp,ytmp,dp(i)
        read(16,*)j,xl(i),yl(i) !,dp(i)
        read(15,*)j,xtmp,ytmp,dtmp
        !dtmp is double
        if(abs(dtmp)>1.d-14) then
          include2(i)=1 !nint(tmp)
        else
          include2(i)=0
        endif
      enddo !i
      do i=1,ne
        read(14,*)j,i34(i),(elnode(l,i),l=1,i34(i))
      enddo !i
!     Open bnds
!      read(14,*) nope
!      read(14,*) neta
!      ntot=0
!      if(nope>mnope) stop 'Increase mnope (2)'
!      do k=1,nope
!        read(14,*) nond(k)
!        if(nond(k)>mnond) stop 'Increase mnond'
!        do i=1,nond(k)
!          read(14,*) iond(k,i)
!        enddo
!      enddo

!     Expand nudging marker to neighbors (to account for elem)
      imap=include2 !temp save
      do i=1,ne
        if(maxval(include2(elnode(1:i34(i),i)))>0) then
          imap(elnode(1:i34(i),i))=1
        endif
      enddo !i
      include2=imap

      close(14)
      close(15)
      close(16)

      !Output
      open(15,file='include3.gr3',status='replace')
      write(15,*); write(15,*)ne,np
      icount=0
      do i=1,np
        write(15,*)i,xl(i),yl(i),include2(i)
        if(include2(i)>0) icount=icount+1
      enddo !i
      print*, icount,' nodes included in _nu.nc'
      do i=1,ne
        write(15,*)i,i34(i),(elnode(l,i),l=1,i34(i))
      enddo !i
      close(15)

!     V-grid
      read(19,*)ivcor
      read(19,*)nvrt
      rewind(19)
      allocate(ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),z(np,nvrt),kbp2(np),tempout(nvrt,np),saltout(nvrt,np))
      call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp2)

!     Compute z-coord.
      do i=1,np
        if(ivcor==2) then
          call zcor_SZ_single(max(0.11,dp(i)),0.,0.1,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma,z(i,:),idry,kbp2(i))
        else if(ivcor==1) then !impose min h
          z(i,kbp2(i):nvrt)=max(0.11,dp(i))*sigma_lcl(kbp2(i):nvrt,i)
        else
          write(11,*)'Unknown ivcor:',ivcor
          stop
        endif

        !Extend below bottom for interpolation later
        z(i,1:kbp2(i)-1)=z(i,kbp2(i))
      enddo !i

!start11
!     Define limits for variables for sanity checks
      tempmin=-15
      tempmax=50
      saltmin=0
      saltmax=45.2

!     Assume S,T have same dimensions (and time step) and grids do not change over time
      timeout=-dtout !sec
      irecout=0
      irecout2=0 !final time record # in the output
!      iyr_now=istart_year
!      imon_now=istart_mon
!      iday_now=istart_day-1

      do ifile=1,nndays
!--------------------------------------------------------------------
      write(char3,'(i20)')ifile
      char3=adjustl(char3); len_char3=len_trim(char3)

!     Open nc file
      ncfile1='TS_'//char3(1:len_char3)//'.nc'
      print*, 'doing ',trim(ncfile1)
      status = nf90_open(trim(ncfile1),nf90_nowrite, sid)
      call check(status)

!     Get dimnesions from first file, assumed same for all variables
      if(ifile==1) then
        status = nf90_inq_varid(sid, "salinity", svid)
        status = nf90_inq_varid(sid, "temperature", tvid)
        write(20,*)'Done reading variable ID'

!       WARNING: indices reversed from ncdump!
        status = nf90_Inquire_Variable(sid, svid, dimids = dids)
        status = nf90_Inquire_Dimension(sid, dids(1), len = ixlen)
        status = nf90_Inquire_Dimension(sid, dids(2), len = iylen)
        status = nf90_Inquire_Dimension(sid, dids(3), len = ilen)
        status = nf90_Inquire_Dimension(sid, dids(4), len = ntime)
        print*, 'ixlen,iylen,ilen,ntime= ',ixlen,iylen,ilen,ntime

!       allocate arrays
        allocate(xind(ixlen),stat=ier)
        allocate(yind(iylen),stat=ier)
        allocate(lind(ilen),stat=ier)
        allocate(lat(ixlen,iylen))
        allocate(lon(ixlen,iylen))
        allocate(zm(ixlen, iylen, ilen))
        allocate(hnc(ilen))
        allocate(kbp(ixlen,iylen))
        allocate(ihope(ixlen,iylen))
        allocate(salt(ixlen,iylen,ilen),stat=ier)
        allocate(temp(ixlen,iylen,ilen),stat=ier)
!        allocate(salt0(ixlen,iylen,ilen),stat=ier)
!        allocate(temp0(ixlen,iylen,ilen),stat=ier)

!       get static info (lat/lon grids etc)
        status = nf90_inq_varid(sid, "depth", hvid)
        status = nf90_get_var(sid, hvid, hnc)
        status = nf90_inq_varid(sid, "xlon", xvid)
        status = nf90_get_var(sid, xvid, xind)
        status = nf90_inq_varid(sid, "ylat", yvid)
        status = nf90_get_var(sid, yvid, yind)

!       processing static info
!       lon/lat as 2D arrays mostly for potential extension to
!       non-rectangular grids
        do i=1,ixlen
          lon(i,:)=xind(i)
          if(i<ixlen) then; if(xind(i)>=xind(i+1)) then
            write(11,*)'Lon must be increasing:',i,xind(i),xind(i+1)
            stop
          endif; endif;
        enddo !i
        do j=1,iylen
          lat(:,j)=yind(j)
          if(j<iylen) then; if(yind(j)>=yind(j+1)) then
            write(11,*)'Lat must be increasing:',j,yind(j),yind(j+1)
            stop
          endif; endif;
        enddo !j
!        lon=lon-360 !convert to our long.

!       Grid bounds for searching for parents  later
!       Won't work across dateline
        xlmin=minval(lon); xlmax=maxval(lon)
        ylmin=minval(lat); ylmax=maxval(lat)

!       Compute z-coord. (assuming eta=0)
!       WARNING: In zm(), 1 is bottom; ilen is surface (SCHISM convention)
        do i=1,ixlen
          do j=1,iylen
            do k=1,ilen
              zm(i,j,k)=-hnc(1+ilen-k)
            enddo !k
          enddo !j
        enddo !i

!       Arrays no longer used after this: hnc,lind
        deallocate(hnc,lind)
      endif !ifile==1

!     Vertical convention follows SCHISM from now on; i.e., 1 is at bottom
      if(ifile==1) then
        ilo=1 !align time origin
      else
        ilo=1
      endif !ifile

      do it2=ilo,ntime !time records in each nc file
        irecout=irecout+1
        timeout=timeout+dtout !sec
        write(20,*)'Time out (days)=',timeout/86400

        !Make sure it2=ilo is output (some init. work below)
        if(mod(irecout-1,nt_out)/=0) cycle
        irecout2=irecout2+1

!       Read T,S
!       WARNING! Make sure the order of vertical indices is 1 at bottom, ilen at surface; revert if necessary!
        call cpu_time(tt0)
        status = nf90_get_var(sid,svid,salt(:,:,ilen:1:-1),start=(/1,1,1,it2/),count=(/ixlen,iylen,ilen,1/))
        status = nf90_get_var(sid,tvid,temp(:,:,ilen:1:-1),start=(/1,1,1,it2/),count=(/ixlen,iylen,ilen,1/))

        !Scaling etc
        salt=salt*1.e-3+20
        temp=temp*1.e-3+20
        !Define junk value for sid; the test is
        !salt<rjunk+0.1
        rjunk=-3.e4*1.e-3+20

!       done read nc this step
        call cpu_time(tt1)
        write(20,*)'reading nc this step took (sec):',tt1-tt0
        call flush(20)

!       Do some prep for 1st step only
        if(ifile==1.and.it2==ilo) then
!         Assume these won't change over time iteration!!
!         Find bottom index and extend
          ndrypt=0 !# of dry nodes in nc
          do i=1,ixlen
            do j=1,iylen
              if(salt(i,j,ilen)<rjunk+0.1) then
                kbp(i,j)=-1 !dry
                ndrypt=ndrypt+1
              else !wet
                !Extend near bottom
                klev0=-1 !flag
                do k=1,ilen
                  if(salt(i,j,k)>rjunk) then
                    klev0=k; exit
                  endif
                enddo !k
                if(klev0<=0) then
                  write(11,*)'Impossible (1):',i,j,salt(i,j,ilen)
                  stop
                endif !klev0
                salt(i,j,1:klev0-1)=salt(i,j,klev0)
                temp(i,j,1:klev0-1)=temp(i,j,klev0)
                kbp(i,j)=klev0 !>0; <=ilen; may change over time (but /=-1?)

                !Check
                do k=1,ilen
                  if(salt(i,j,k)<saltmin.or.salt(i,j,k)>saltmax.or. &
                    &temp(i,j,k)<tempmin.or.temp(i,j,k)>tempmax) then
                    write(11,*)'Fatal: no valid values:',it2,i,j,k,salt(i,j,k),temp(i,j,k)
                    stop
                  endif
                enddo !k
              endif !salt
            enddo !j
          enddo !i
!          print*, 'ndrypt=',ndrypt

!         Compute S,T etc@ invalid pts based on nearest neighbor
!         Search around neighborhood of a pt
          allocate(iparen_of_dry(2,max(1,ndrypt)))
          call cpu_time(tt0)

          icount=0
          do i=1,ixlen
            do j=1,iylen
              if(kbp(i,j)==-1) then !invalid pts
                icount=icount+1
                if(icount>ndrypt) stop 'overflow'
                !Compute max possible tier # for neighborhood
                mmax=max(i-1,ixlen-i,j-1,iylen-j)

                m=0 !tier #
                loop6: do
                  m=m+1
                  do ii=max(-m,1-i),min(m,ixlen-i)
                    i3=max(1,min(ixlen,i+ii))
                    do jj=max(-m,1-j),min(m,iylen-j)
                      j3=max(1,min(iylen,j+jj))
                      if(kbp(i3,j3)>0) then !found
                        i1=i3; j1=j3
                        exit loop6
                      endif
                    enddo !jj
                  enddo !ii

                  if(m==mmax) then
                    write(11,*)'Max. exhausted:',i,j,mmax
                    write(11,*)'kbp'
                    do ii=1,ixlen
                      do jj=1,iylen
                        write(11,*)ii,jj,kbp(ii,jj)
                      enddo !jj
                    enddo !ii
                    stop
                  endif
                end do loop6

                salt(i,j,1:ilen)=salt(i1,j1,1:ilen)
                temp(i,j,1:ilen)=temp(i1,j1,1:ilen)

                !Save for other steps
                iparen_of_dry(1,icount)=i1
                iparen_of_dry(2,icount)=j1

                !Check
                do k=1,ilen
                  if(salt(i,j,k)<saltmin.or.salt(i,j,k)>saltmax.or. &
                    &temp(i,j,k)<tempmin.or.temp(i,j,k)>tempmax) then
                    write(11,*)'Fatal: no valid values after searching:',it2,i,j,k,salt(i,j,k),temp(i,j,k)
                    stop
                  endif
                enddo !k

              endif !kbp(i,j)==-1
            enddo !j=iylen1,iylen2
          enddo !i=ixlen1,ixlen2
          call cpu_time(tt1)
          if(ndrypt/=icount) stop 'mismatch(7)'
          write(20,*)'extending took (sec):',tt1-tt0,ndrypt,icount
          call flush(20)

          !Save T/S for abnormal cases later
!          salt0=salt
!          temp0=temp

          call cpu_time(tt0)
!         Test outputs
          icount=0
          do i=1,ixlen
            do j=1,iylen
              icount=icount+1
              write(98,*)icount,lon(i,j),lat(i,j),salt(i,j,1) !,i,j
              write(99,*)icount,lon(i,j),lat(i,j),temp(i,j,ilen)
              write(100,*)icount,lon(i,j),lat(i,j),temp(i,j,1)
            enddo !j
          enddo !i
          print*, 'done outputting test outputs for nc'

!         Find parent elements for hgrid.ll
          ixy=0
          loop4: do i=1,np
            !won't work across dateline
            if(xl(i)<xlmin.or.xl(i)>xlmax.or.yl(i)<ylmin.or.yl(i)>ylmax) cycle loop4
            if(include2(i)==0) cycle loop4

            if(interp_mode==0) then !SG search
              do ix=1,ixlen-1
                if(xl(i)>=xind(ix).and.xl(i)<=xind(ix+1)) then
                  !Lower left corner index
                  ixy(i,1)=ix
                  xrat=(xl(i)-xind(ix))/(xind(ix+1)-xind(ix))
                  exit
                endif
              enddo !ix
              do iy=1,iylen-1
                if(yl(i)>=yind(iy).and.yl(i)<=yind(iy+1)) then
                  !Lower left corner index
                  ixy(i,2)=iy
                  yrat=(yl(i)-yind(iy))/(yind(iy+1)-yind(iy))
                  exit
                endif
              enddo !ix

              if(ixy(i,1)==0.or.ixy(i,2)==0) then
                write(11,*)'Did not find parent:',i,ixy(i,1:2)
                stop
              endif
              if(xrat<0.or.xrat>1.or.yrat<0.or.yrat>1) then
                write(11,*)'Ratio out of bound:',i,xrat,yrat
                stop
              endif

              !Bilinear shape function
              arco(1,i)=(1-xrat)*(1-yrat)
              arco(2,i)=xrat*(1-yrat)
              arco(4,i)=(1-xrat)*yrat
              arco(3,i)=xrat*yrat
            else !interp_mode=1; generic search with UG
              do ix=1,ixlen-1
                do iy=1,iylen-1
                  x1=lon(ix,iy); x2=lon(ix+1,iy); x3=lon(ix+1,iy+1); x4=lon(ix,iy+1)
                  y1=lat(ix,iy); y2=lat(ix+1,iy); y3=lat(ix+1,iy+1); y4=lat(ix,iy+1)
                  a1=abs(signa_single(xl(i),x1,x2,yl(i),y1,y2))
                  a2=abs(signa_single(xl(i),x2,x3,yl(i),y2,y3))
                  a3=abs(signa_single(xl(i),x3,x4,yl(i),y3,y4))
                  a4=abs(signa_single(xl(i),x4,x1,yl(i),y4,y1))
                  b1=abs(signa_single(x1,x2,x3,y1,y2,y3))
                  b2=abs(signa_single(x1,x3,x4,y1,y3,y4))
                  rat=abs(a1+a2+a3+a4-b1-b2)/(b1+b2)
                  if(rat<small1) then
                    ixy(i,1)=ix; ixy(i,2)=iy
!                   Find a triangle
                    in=0 !flag
                    do l=1,2 !split quad
                      ap=abs(signa_single(xl(i),x1,x3,yl(i),y1,y3))
                      if(l==1) then !nodes 1,2,3
                        bb=abs(signa_single(x1,x2,x3,y1,y2,y3))
                        wild(l)=abs(a1+a2+ap-bb)/bb
                        if(wild(l)<small1*5) then
                          in=1
                          arco(1,i)=max(0.,min(1.,a2/bb))
                          arco(2,i)=max(0.,min(1.,ap/bb))
                          arco(3,i)=max(0.,min(1.,1-arco(1,i)-arco(2,i)))
                          arco(4,i)=0.
                          exit
                        endif
                      else !nodes 1,3,4
                        bb=abs(signa_single(x1,x3,x4,y1,y3,y4))
                        wild(l)=abs(a3+a4+ap-bb)/bb
                        if(wild(l)<small1*5) then
                          in=2
                          arco(1,i)=max(0.,min(1.,a3/bb))
                          arco(3,i)=max(0.,min(1.,a4/bb))
                          arco(4,i)=max(0.,min(1.,1-arco(1,i)-arco(3,i)))
                          arco(2,i)=0.
                          exit
                        endif
                      endif
                    enddo !l=1,2
                    if(in==0) then
                      write(11,*)'Cannot find a triangle:',(wild(l),l=1,2)
                      stop
                    endif
                    !ixy(i,3)=in
                    cycle loop4
                  endif !rat<small1
                enddo !iy=iylen1,iylen2-1
              enddo !ix=ixlen1,ixlen2-1
            endif !interp_mode
          end do loop4 !i=1,np

          call cpu_time(tt1)
          write(20,*)'computing weights took:',tt1-tt0
          call flush(20)

        else !skip to extend directly
          do i=1,ixlen
            do j=1,iylen
              if(kbp(i,j)>0) then !valid
                !Extend bottom - kbp changes over time
                klev0=-1 !flag
                do k=1,ilen
                  if(salt(i,j,k)>rjunk) then
                    klev0=k; exit
                  endif
                enddo !k
                if(klev0<=0) then
                  write(11,*)'Impossible (4):',i,j,salt(i,j,ilen)
                  stop
                endif !klev0
                salt(i,j,1:klev0-1)=salt(i,j,klev0)
                temp(i,j,1:klev0-1)=temp(i,j,klev0)

                !Check
                do k=1,ilen
                  if(salt(i,j,k)<saltmin.or.salt(i,j,k)>saltmax.or. &
                    &temp(i,j,k)<tempmin.or.temp(i,j,k)>tempmax) then
                    write(11,*)'Warning: no valid values-use last value(3):', &
     &it2,i,j,k,salt(i,j,k),temp(i,j,k),';',salt(i,j,:),temp(i,j,:)
!'                   salt(i,j,k)=salt0(i,j,k)
!                    temp(i,j,k)=temp0(i,j,k)
                    stop
                  endif
                enddo !k
              endif !kbp
            enddo !j
          enddo !i

          icount=0
          do i=1,ixlen
            do j=1,iylen
              if(kbp(i,j)==-1) then !dry pt
                icount=icount+1
                !Use extended parents
                i1=iparen_of_dry(1,icount); j1=iparen_of_dry(2,icount)
                salt(i,j,1:ilen)=salt(i1,j1,1:ilen)
                temp(i,j,1:ilen)=temp(i1,j1,1:ilen)

                !Check
                do k=1,ilen
                  if(salt(i,j,k)<saltmin.or.salt(i,j,k)>saltmax.or. &
                    &temp(i,j,k)<tempmin.or.temp(i,j,k)>tempmax) then
                    write(11,*)'Warning: no valid values-use last value(4):', &
     &it2,i,j,k,salt(i,j,k),temp(i,j,k),i1,j1,';', &
     &salt(i1,j1,:),temp(i1,j1,:)
!'                   salt(i,j,k)=salt0(i,j,k)
!                    temp(i,j,k)=temp0(i,j,k)
                    stop
                  endif
                enddo !k
              endif !kbp
            enddo !j
          enddo !i
          if(icount/=ndrypt) stop 'mimatch (8)'

          !Save for later steps
!          salt0=salt
!          temp0=temp
        endif !ifile==1.and.it2==

!       Do interpolation
!        print*, 'Time out (days)=',timeout/86400,it2,ilo,ntime
        tempout=tem_outside; saltout=sal_outside !init
        npout=0 !# of _valid_ output points
        do i=1,np
          if(ixy(i,1)==0.or.ixy(i,2)==0) then
          else !found parent
            npout=npout+1
            imap(npout)=i
            ix=ixy(i,1); iy=ixy(i,2) !in=ixy(i,3)
            !Find vertical level
            do k=1,nvrt
              if(kbp(ix,iy)==-1) then
                lev=ilen-1; vrat=1
              else if(z(i,k)<=zm(ix,iy,kbp(ix,iy))) then
                lev=kbp(ix,iy); vrat=0
              else if(z(i,k)>=zm(ix,iy,ilen)) then !above f.s.
                lev=ilen-1; vrat=1
              else
                lev=-99 !flag
                do kk=1,ilen-1
                  if(z(i,k)>=zm(ix,iy,kk).and.z(i,k)<=zm(ix,iy,kk+1)) then
                    lev=kk
                    vrat=(zm(ix,iy,kk)-z(i,k))/(zm(ix,iy,kk)-zm(ix,iy,kk+1))
                    exit
                  endif
                enddo !kk
                if(lev==-99) then
                  write(11,*)'Cannot find a level:',i,k,z(i,k),(zm(ix,iy,l),l=1,kbp(ix,iy))
                  stop
                endif
              endif

!              write(18,*)i,k,ix,iy,lev,vrat,kbp(ix,iy)
              !Impose bounds for odd cases
              lev2=lev+1
              lev=max(1,min(ilen,lev))
              lev2=max(1,min(ilen,lev2))

              wild2(1,1)=temp(ix,iy,lev)*(1-vrat)+temp(ix,iy,lev2)*vrat
              wild2(1,2)=salt(ix,iy,lev)*(1-vrat)+salt(ix,iy,lev2)*vrat
              wild2(2,1)=temp(ix+1,iy,lev)*(1-vrat)+temp(ix+1,iy,lev2)*vrat
              wild2(2,2)=salt(ix+1,iy,lev)*(1-vrat)+salt(ix+1,iy,lev2)*vrat
              wild2(3,1)=temp(ix+1,iy+1,lev)*(1-vrat)+temp(ix+1,iy+1,lev2)*vrat
              wild2(3,2)=salt(ix+1,iy+1,lev)*(1-vrat)+salt(ix+1,iy+1,lev2)*vrat
              wild2(4,1)=temp(ix,iy+1,lev)*(1-vrat)+temp(ix,iy+1,lev2)*vrat
              wild2(4,2)=salt(ix,iy+1,lev)*(1-vrat)+salt(ix,iy+1,lev2)*vrat

              tempout(k,i)=dot_product(wild2(1:4,1),arco(1:4,i))
              saltout(k,i)=dot_product(wild2(1:4,2),arco(1:4,i))
!              if(in==1) then
!                tempout(k,i)=wild2(1,1)*arco(1,i)+wild2(2,1)*arco(2,i)+wild2(3,1)*arco(3,i)
!                saltout(k,i)=wild2(1,2)*arco(1,i)+wild2(2,2)*arco(2,i)+wild2(3,2)*arco(3,i)
!              else
!                tempout(k,i)=wild2(1,1)*arco(1,i)+wild2(3,1)*arco(2,i)+wild2(4,1)*arco(3,i)
!                saltout(k,i)=wild2(1,2)*arco(1,i)+wild2(3,2)*arco(2,i)+wild2(4,2)*arco(3,i)
!              endif

              !Check
              if(tempout(k,i)<tempmin.or.tempout(k,i)>tempmax.or. &
               &saltout(k,i)<saltmin.or.saltout(k,i)>saltmax) then
                write(11,*)'Interpolated values invalid:',i,k,tempout(k,i),saltout(k,i)
                stop
              endif

!             Enforce lower bound for temp. for eqstate
              tempout(k,i)=max(0.,tempout(k,i))
            enddo !k=1,nvrt
          endif !ixy(i,1)==0.or.
        enddo !i=1,np

!       Output (junks outside nudging zone)
        print*, 'outputting _nu.nc at day ',timeout/86400
!       For debug
        if(ifile==1.and.it2==ilo) then
          do i=1,np
            write(101,*)i,xl(i),yl(i),tempout(nvrt,i)
            write(102,*)i,xl(i),yl(i),tempout(1,i)
            write(103,*)i,xl(i),yl(i),saltout(nvrt,i)
          enddo !i

          !write(*,*)'# of output nodes=',npout
          write(104,*); write(104,*)npout
          do i=1,npout
            nd=imap(i)
            write(104,*)i,xl(nd),yl(nd),0.
          enddo !i
        endif !ifile
        close(101); close(102); close(103)

        if(ifile==1.and.it2==ilo) then
          iret=nf90_create('TEM_nu.nc',OR(NF90_CLOBBER,NF90_NETCDF4),ncid_T)
          iret=nf90_def_dim(ncid_T,'node',npout,node_dim)
          iret=nf90_def_dim(ncid_T,'nLevels',nvrt,nv_dim)
          iret=nf90_def_dim(ncid_T,'one',1,one_dim)
          iret=nf90_def_dim(ncid_T,'time', NF90_UNLIMITED,itime_dim)

          var1d_dims(1)=itime_dim
          iret=nf90_def_var(ncid_T,'time',NF90_DOUBLE,var1d_dims,itime_id_T)
          var1d_dims(1)=node_dim
          iret=nf90_def_var(ncid_T,'map_to_global_node',NF90_INT,var1d_dims,id_map_T)
          var4d_dims(1)=one_dim; var4d_dims(2)=nv_dim;
          var4d_dims(3)=node_dim; var4d_dims(4)=itime_dim
          iret=nf90_def_var(ncid_T,'tracer_concentration',NF90_FLOAT,var4d_dims,ivar_T)
          iret=nf90_enddef(ncid_T)

          iret=nf90_create('SAL_nu.nc',OR(NF90_CLOBBER,NF90_NETCDF4),ncid_S)
          iret=nf90_def_dim(ncid_S,'node',npout,node_dim)
          iret=nf90_def_dim(ncid_S,'nLevels',nvrt,nv_dim)
          iret=nf90_def_dim(ncid_S,'one',1,one_dim)
          iret=nf90_def_dim(ncid_S,'time', NF90_UNLIMITED,itime_dim)

          var1d_dims(1)=itime_dim
          iret=nf90_def_var(ncid_S,'time',NF90_DOUBLE,var1d_dims,itime_id_S)
          var1d_dims(1)=node_dim
          iret=nf90_def_var(ncid_S,'map_to_global_node',NF90_INT,var1d_dims,id_map_S)
          iret=nf90_def_var(ncid_S,'tracer_concentration',NF90_FLOAT,var4d_dims,ivar_S)
          iret=nf90_enddef(ncid_S)
        endif !ifile

        aa1(1)=timeout/86400
        iret=nf90_put_var(ncid_T,itime_id_T,aa1,(/irecout2/),(/1/))
        iret=nf90_put_var(ncid_S,itime_id_S,aa1,(/irecout2/),(/1/))
        iret=nf90_put_var(ncid_T,id_map_T,imap(1:npout),(/1/),(/npout/))
        iret=nf90_put_var(ncid_S,id_map_S,imap(1:npout),(/1/),(/npout/))
        iret=nf90_put_var(ncid_T,ivar_T,tempout(:,imap(1:npout)),(/1,1,1,irecout2/),(/1,nvrt,npout,1/))
        iret=nf90_put_var(ncid_S,ivar_S,saltout(:,imap(1:npout)),(/1,1,1,irecout2/),(/1,nvrt,npout,1/))
      enddo !it2=ilo,ntime

      status = nf90_close(sid)
!end11
      print*, 'done reading nc for file ',ifile
!--------------------------------------------------------------------
      enddo !ifile=1,

      iret=nf90_close(ncid_T)
      iret=nf90_close(ncid_S)
!      deallocate(lat,lon,zm,h,kbp,ihope,xind,yind,lind,salt,temp)

      print*, 'Finished'
!     End of main

!     Subroutines
      contains
!     Internal subroutine - checks error status after each netcdf,
!     prints out text message each time an error code is returned.
      subroutine check(status)
      integer, intent ( in) :: status

      if(status /= nf90_noerr) then
        print *, trim(nf90_strerror(status))
        print *, 'failed to open nc files'
        stop
      end if
      end subroutine check

      end program gen_hot

max_stars_repo_name: cseaton/schism

max_issues_repo_name: cseaton/schism

max_forks_repo_name: cseaton/schism


text: #redirect Wiscombes Davis Funeral Chapel

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search