text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1995  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ################################################################
c     ##                                                            ##
c     ##  program gda  --  simulated annealing on gaussian density  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "gda" implements Gaussian Density Annealing (GDA) algorithm
c     for global optimization via simulated annealing
c
c     literature reference:
c
c     J. Ma and J. E. Straub, "Simulated Annealing using the
c     Classical Density Distribution", Journal of Chemical Physics,
c     101, 533-541 (1994)
c
c
      program gda
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'files.i'
      include 'iounit.i'
      include 'minima.i'
      include 'potent.i'
      include 'vdwpot.i'
      include 'warp.i'
      integer i,igda,itrial,ntrial
      integer nstep,nvar,nok,nbad
      integer lext,next,freeunit
      real*8 bstart,bstop
      real*8 random,boxsize
      real*8 eps,h1,hmin,gda2
      real*8 minimum,grdmin
      real*8 xcm,ycm,zcm
      real*8, allocatable :: m2init(:)
      real*8, allocatable :: xx(:)
      logical exist,randomize
      character*1 answer
      character*6 mode,method
      character*7 ext,status
      character*120 gdafile
      character*120 record
      character*120 string
      external gda1,gda2,gda3
      external optsave
c
c
c     set up the structure, mechanics calculation and smoothing
c
      call initial
      call getxyz
      use_smooth = .true.
      use_gda = .true.
      call mechanic
c
c     get the number of optimized structures to be constructed
c
      ntrial = 0
      call nextarg (string,exist)
      if (exist)  read (string,*,err=10,end=10)  ntrial
   10 continue
      if (ntrial .le. 0) then
         write (iout,20)
   20    format (/,' Enter Number of Annealing Trials [1] :  ',$)
         read (input,30)  ntrial
   30    format (i10)
      end if
      if (ntrial .le. 0)  ntrial = 1
c
c     see if random coordinates are desired as starting structures
c
      randomize = .true.
      if (ntrial .eq. 1) then
         randomize = .false.
         call nextarg (answer,exist)
         if (.not. exist) then
            write (iout,40)
   40       format (/,' Use Randomized Initial Coordinates [N] :  ',$)
            read (input,50)  record
   50       format (a120)
            next = 1
            call gettext (record,answer,next)
         end if
         call upcase (answer)
         if (answer .eq. 'Y')  randomize = .true.
      end if
      if (randomize)  boxsize = 10.0d0 * (dble(n))**(1.0d0/3.0d0)
c
c     get initial and final values of inverse temperature
c
      bstart = -1.0d0
      bstop = -1.0d0
      call nextarg (string,exist)
      if (exist)  read (string,*,err=60,end=60)  bstart
      call nextarg (string,exist)
      if (exist)  read (string,*,err=60,end=60)  bstop
   60 continue
      if (bstart.le.0.0d0 .or. bstop.le.0.0d0) then
         write (iout,70)
   70    format (/,' Enter Initial and Final Beta [0.01, 10**10] :  ',$)
         read (input,80)  record
   80    format (a120)
         read (record,*,err=90,end=90)  bstart,bstop
   90    continue
      end if
      if (bstart .le. 0.0d0)  bstart = 0.01d0
      if (bstop .le. 0.0d0)  bstop = 1.0d10
c
c     perform dynamic allocation of some local arrays
c
      allocate (m2init(n))
      allocate (xx(4*n))
c
c     store the initial values of the squared mean Gaussian width
c
      do i = 1, n
         m2init(i) = m2(1)
      end do
c
c     write out a copy of coordinates for later update
c
      do itrial = 1, ntrial
         lext = 3
         call numeral (itrial,ext,lext)
         gdafile = filename(1:leng)//'.'//ext(1:lext)
         call version (gdafile,'new')
         igda = freeunit ()
         open (unit=igda,file=gdafile,status='new')
         call prtxyz (igda)
         close (unit=igda)
         outfile = gdafile
c
c     set an initial box size and generate random coordinates
c
         if (randomize) then
            do i = 1, n
               x(i) = boxsize * random ()
               y(i) = boxsize * random ()
               z(i) = boxsize * random ()
            end do
         end if
c
c     translate coordinates and M2's to optimization parameters
c
         nvar = 0
         do i = 1, n
            nvar = nvar + 1
            xx(nvar) = x(i)
            nvar = nvar + 1
            xx(nvar) = y(i)
            nvar = nvar + 1
            xx(nvar) = z(i)
         end do
         do i = 1, n
            nvar = nvar + 1
            xx(nvar) = m2init(i)
         end do
c
c     make changes to the potential to use potential smoothing
c
         use_smooth = .true.
         use_geom = .true.
         vdwtyp = 'GAUSSIAN'
c
c     make the call to the Bulirsch-Stoer integration routine
c
         nstep = 0
         status = '       '
         eps = 1.0d-8
         h1 = 0.01d0
         hmin = 0.0d0
         write (iout,100)
  100    format (//,' Gaussian Density Annealing Global Optimization :',
     &           //,' BS Step',5x,'Log(Beta)',6x,'Energy',
     &              9x,'Rg',8x,'Log(M2)',7x,'Status',/)
         call gdastat (nstep,bstart,xx,status)
         call diffeq (nvar,xx,bstart,bstop,eps,h1,hmin,nok,nbad,gda1)
         nstep = nok + nbad
c
c     make changes to the potential for standard optimization
c
         use_smooth = .false.
         use_geom = .false.
         vdwtyp = 'LENNARD-JONES'
c
c     make the call to the energy minimization routine
c
         mode = 'DTNCG'
         method = 'AUTO'
         nvar = 3 * n
         grdmin = 0.0001d0
         nextiter = nstep + 1
         call tncg (mode,method,nvar,xx,minimum,grdmin,
     &                     gda2,gda3,optsave)
c        call lbfgs (nvar,xx,minimum,grdmin,gda2,optsave)
         write (iout,110)  itrial,minimum
  110    format (/,' Global Energy Minimum for Trial',i4,' :',f15.4)
c
c     translate optimization parameters into atomic coordinates
c
         nvar = 0
         do i = 1, n
            nvar = nvar + 1
            x(i) = xx(nvar)
            nvar = nvar + 1
            y(i) = xx(nvar)
            nvar = nvar + 1
            z(i) = xx(nvar)
         end do
c
c     move the center of mass to the origin
c
         xcm = 0.0d0
         ycm = 0.0d0
         zcm = 0.0d0
         do i = 1, n
            xcm = xcm + x(i)
            ycm = ycm + y(i)
            zcm = zcm + z(i)
         end do
         xcm = xcm / dble(n)
         ycm = ycm / dble(n)
         zcm = zcm / dble(n)
         do i = 1, n
            x(i) = x(i) - xcm
            y(i) = y(i) - xcm
            z(i) = z(i) - xcm
         end do
c
c     write the final coordinates into a file
c
         igda = freeunit ()
         open (unit=igda,file=gdafile,status='old')
         rewind (unit=igda)
         call prtxyz (igda)
         close (unit=igda)
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (m2init)
      deallocate (xx)
c
c     perform any final tasks before program exit
c
      call final
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine gda1  --  gaussian density annealing gradient  ##
c     ##                                                            ##
c     ################################################################
c
c
      subroutine gda1 (beta,xx,g)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'iounit.i'
      include 'warp.i'
      integer i,nvar
      integer, allocatable :: hinit(:)
      integer, allocatable :: hstop(:)
      integer, allocatable :: hindex(:)
      real*8 beta,e,sum
      real*8, allocatable :: hdiag(:)
      real*8, allocatable :: h(:)
      real*8 xx(*)
      real*8 g(*)
      real*8, allocatable :: derivs(:,:)
c
c
c     translate optimization parameters to coordinates and M2's
c
      nvar = 0
      do i = 1, n
         nvar = nvar + 1
         x(i) = xx(nvar)
         nvar = nvar + 1
         y(i) = xx(nvar)
         nvar = nvar + 1
         z(i) = xx(nvar)
      end do
      do i = 1, n
         nvar = nvar + 1
         m2(i) = xx(nvar)
         if (m2(i) .lt. 0.0d0) then
            write (iout,10)  i,m2(i)
   10       format (' GDA1  --  Warning, Negative M2 at Atom',i6,
     &                  ' with Value',d12.4)
            m2(i) = -m2(i)
         end if
      end do
c
c     perform dynamic allocation of some local arrays
c
      allocate (derivs(3,n))
c
c     compute and store the Cartesian energy gradient vector
c
      call gradient (e,derivs)
c
c     translate the energy gradient into a dr/dbeta vector
c
      nvar = 0
      do i = 1, n
         nvar = nvar + 1
         g(nvar) = -(m2(i)/3.0d0) * derivs(1,i)
         nvar = nvar + 1
         g(nvar) = -(m2(i)/3.0d0) * derivs(2,i)
         nvar = nvar + 1
         g(nvar) = -(m2(i)/3.0d0) * derivs(3,i)
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (derivs)
c
c     perform dynamic allocation of some local arrays
c
      allocate (hinit(3*n))
      allocate (hstop(3*n))
      allocate (hindex((3*n*(3*n-1))/2))
      allocate (hdiag(3*n))
      allocate (h((3*n*(3*n-1))/2))
c
c     compute and store the Hessian elements
c
      call hessian (h,hinit,hstop,hindex,hdiag)
c
c     translate the Hessian diagonal into a dM2/dbeta vector
c
      do i = 1, n
         nvar = nvar + 1
         sum = hdiag(3*i-2) + hdiag(3*i-1) + hdiag(3*i)
         g(nvar) = -(m2(i)/3.0d0)**2 * sum
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (hinit)
      deallocate (hstop)
      deallocate (hindex)
      deallocate (hdiag)
      deallocate (h)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  function gda2  --  energy/gradient for TNCG optimization  ##
c     ##                                                            ##
c     ################################################################
c
c
      function gda2 (xx,g)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      integer i,nvar
      real*8 gda2,e
      real*8 xx(*)
      real*8 g(*)
      real*8, allocatable :: derivs(:,:)
c
c
c     translate optimization parameters to atomic coordinates
c
      nvar = 0
      do i = 1, n
         nvar = nvar + 1
         x(i) = xx(nvar)
         nvar = nvar + 1
         y(i) = xx(nvar)
         nvar = nvar + 1
         z(i) = xx(nvar)
      end do
c
c     perform dynamic allocation of some local arrays
c
      allocate (derivs(3,n))
c
c     compute and store the energy and gradient
c
      call gradient (e,derivs)
      gda2 = e
c
c     store Cartesian gradient as optimization gradient
c
      nvar = 0
      do i = 1, n
         nvar = nvar + 1
         g(nvar) = derivs(1,i)
         nvar = nvar + 1
         g(nvar) = derivs(2,i)
         nvar = nvar + 1
         g(nvar) = derivs(3,i)
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (derivs)
      return
      end
c
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine gda3  --  Hessian values for TNCG optimization  ##
c     ##                                                             ##
c     #################################################################
c
c
      subroutine gda3 (mode,xx,h,hinit,hstop,hindex,hdiag)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      integer i,nvar
      integer hinit(*)
      integer hstop(*)
      integer hindex(*)
      real*8 xx(*)
      real*8 hdiag(*)
      real*8 h(*)
      character*4 mode
c
c
c     translate optimization parameters to atomic coordinates
c
      if (mode .eq. 'NONE')  return
      nvar = 0
      do i = 1, n
         nvar = nvar + 1
         x(i) = xx(nvar)
         nvar = nvar + 1
         y(i) = xx(nvar)
         nvar = nvar + 1
         z(i) = xx(nvar)
      end do
c
c     compute and store the Hessian elements
c
      call hessian (h,hinit,hstop,hindex,hdiag)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
C    @(#)fvmov.f	20.3 2/13/96
      subroutine fvmov (a,i,c,k,n)
      real a(*),c(*)

      ia = 1
      ic = 1
      do 100 l=1,n
         c(ic) = a(ia)
         ia = ia + i
         ic = ic + k
  100 continue
      return
      end

ext: f

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 (c) 2004, NVIDIA CORPORATION.  All rights reserved.
!
! 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.
!
! test host's memory arguments in an internal procedure
!

module foo
  contains
    subroutine sub(result,d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14)
      integer result(:)
      integer d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12,d13,d14
      integer l1
      l1 = 25
      call csub
    contains
      subroutine csub()
	result(1) = d1
	result(2) = d2
	result(3) = d3
	result(4) = d4
	result(5) = d5
	result(6) = d6
	result(7) = d7
	result(8) = d8
	result(9) = d9
	result(10) = d10
	result(11) = d11
	result(12) = d12
	result(13) = d13
	result(14) = d14
	result(15) = l1
      end subroutine
    end subroutine
end module

program p
use foo
  integer result(15)
  integer expect(15)
  data expect/ 11,12,13,14,15,16,17,18,19,20,21,22,23,24,25/

  call sub(result,11,12,13,14,15,16,17,18,19,20,21,22,23,24)
  call check(result,expect,15)
end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
module boundary_interface
  integer :: npatch
  integer :: ntrifaces,nquadfaces
  integer, allocatable :: faceInfo(:)
  integer, allocatable :: cell2cell(:)

  integer, dimension(3,4)  :: tetf
  integer, dimension(4,5)  :: pyrf
  integer, dimension(4,5)  :: prizmf
  integer, dimension(4,6)  :: hexf
  integer, dimension(6,4)  :: numverts
  integer, dimension(4)    :: iface
  integer :: tetmap(4)
  integer :: pyrmap(5)
  integer :: prizmap(5)
  integer :: hexmap(6)
  !
  data tetf /1,2,3,1,4,2,2,4,3,1,3,4/                        ! ugrid/mcell face connectivity (tet)
  !data tetmap /4,1,2,3/                                     ! map to exodus face ordering
  data tetmap /1,2,3,4/
  data pyrf /1,2,3,4,1,5,2,2,2,5,3,3,4,3,5,5,1,4,5,5/        ! ugrid/mcell face connectivity (pyramid)
  !data pyrmap /5,1,2,3,4/                                     ! map to exodus face ordering
  data pyrmap /1,2,3,4,5/
  data prizmf /1,2,3,3,1,4,5,2,2,5,6,3,1,3,6,4,4,6,5,5/       ! prizm connectivity
  !data prizmap /4,1,2,3,5/                                    ! map to exodus face ordering
  data prizmap /1,2,3,4,5/
  data hexf /1,2,3,4,1,5,6,2,2,6,7,3,3,7,8,4,1,4,8,5,5,8,7,6/ ! ugrid/mcell connectivity (hex)
  data hexmap /1,2,3,4,5,6/
  !data hexmap /5,1,2,3,4,6/                                   ! map to exodus face ordering
  data iface/4,5,5,6/
  data numverts/3,3,3,3,0,0,&
       4,3,3,3,3,0,&
       3,4,4,4,3,0,&
       4,4,4,4,4,4/
contains
  subroutine deletedata
    if (allocated(faceInfo)) deallocate(faceInfo)
    if (allocated(cell2cell)) deallocate(cell2cell)
  end subroutine deletedata
  !>
  !> find boundary faces in a given graph, using the information of
  !> boundary nodes that are provided in iflag.
  !> If all nodes are set to 1, then the
  !> routine will find all the faces of a given graph.
  !>
  !> include 'findboundaryface.h' for the explicit interface
  !> since this subroutine allocates the data it finds
  !>
  !> Get the connectivity of the faces and also
  !> [cellId,facenum] pair for each face
  !>
  !> Jay Sitaraman
  !> 12/15/2017
  !>
  subroutine findnewboundaryface(cell2cell,faceInfo,iflag,ntrifaces,nquadfaces,ndc4,ndc5,ndc6,ndc8, &
       nnode,ntetra,npyra,nprizm,nhexa)
    !
    implicit none
    !
    integer, allocatable, intent(out) :: cell2cell(:)      !> cell2cell graph
    integer, allocatable, intent(out) :: faceInfo(:)       !> face information (faceconnectivity and cell connections
    integer, intent(inout) :: ntrifaces,nquadfaces         !> number of triangle and quad faces
    integer, intent(in) :: nnode,ntetra,npyra,nprizm,nhexa !> number of nodes and elements in the graph
    integer, intent(in) :: iflag(nnode)                    !> flag to say which one are boundary nodes
    integer, intent(in) :: ndc4(4,ntetra)                  !> tetra graph
    integer, intent(in) :: ndc5(5,npyra)                   !> pyra graph
    integer, intent(in) :: ndc6(6,nprizm)                  !> prizm graph
    integer, intent(in) :: ndc8(8,nhexa)                   !> hex graph
    !
    ! local variables
    !
    !
    integer :: i,j,k,l,nsum,ncells,maxfaces,nfaces,csum,nbface,icptr
    integer :: flocal(4)
    integer, allocatable :: face(:,:),ipointer(:),iflag_cell(:)
    integer, parameter :: base=1
    real*8 :: t1,t2
    logical :: verbose
    !
    ! find all exposed cell faces
    !
    verbose=.false.
    call cpu_time(t1)
    allocate(iflag_cell(ntetra+npyra+nprizm+nhexa))
    !
    maxfaces=0
    k=0
    do i=1,ntetra
       k=k+1
       iflag_cell(k)=0
       csum=0
       do j=1,4
          csum=csum+iflag(ndc4(j,i)+base)
       enddo
       if (csum > 0 ) then
          iflag_cell(k)=1
          maxfaces=maxfaces+4
       endif
    enddo
    !
    do i=1,npyra
       k=k+1
       iflag_cell(k)=0
       csum=0
       do j=1,5
          csum=csum+iflag(ndc5(j,i)+base)
       enddo
       if (csum > 0 ) then
          iflag_cell(k)=1
          maxfaces=maxfaces+5
       endif
    enddo
    !
    do i=1,nprizm
       k=k+1
       iflag_cell(k)=0
       csum=0
       do j=1,6
          csum=csum+iflag(ndc6(j,i)+base)
       enddo
       if (csum > 0 ) then
          iflag_cell(k)=1
          maxfaces=maxfaces+6
       endif
    enddo
    !
    do i=1,nhexa
       k=k+1
       iflag_cell(k)=0
       csum=0
       do j=1,8
          csum=csum+iflag(ndc8(j,i)+base)
       enddo
       if (csum >0 ) then
          iflag_cell(k)=1
          maxfaces=maxfaces+8
       endif
    enddo
    !
    ! do a hash insert of the faces now
    !
    allocate(ipointer(nnode))
    allocate(face(9,maxfaces))
    !
    face=0
    ipointer=0
    nfaces=0
    ncells=0
    !
    l=0
    do i=1,ntetra
       l=l+1
       if (iflag_cell(l)==0) cycle
       do j=1,iface(1)
          nsum=0
          flocal=0
          do k=1,numverts(j,1)
             flocal(k)=ndc4(tetf(k,j),i)+base
          enddo
          call insert_face(flocal,face,ipointer,i+ncells,&
               tetmap(j),nfaces,maxfaces,nnode,numverts(j,1))
       enddo
    enddo
    !
    ncells=ncells+ntetra
    !
    do i=1,npyra
       l=l+1
       if (iflag_cell(l)==0) cycle
       do j=1,iface(2)
          nsum=0
          flocal=0
          do k=1,numverts(j,2)
             flocal(k)=ndc5(pyrf(k,j),i)+base
          enddo
          call insert_face(flocal,face,ipointer,i+ncells,&
               pyrmap(j),nfaces,maxfaces,nnode,numverts(j,2))
       enddo
    enddo
    !
    ncells=ncells+npyra
    !
    do i=1,nprizm
       l=l+1
       if (iflag_cell(l)==0) cycle
       do j=1,iface(3)
          nsum=0
          flocal=0
          do k=1,numverts(j,3)
             flocal(k)=ndc6(prizmf(k,j),i)+base
          enddo
          call insert_face(flocal,face,ipointer,i+ncells,&
               prizmap(j),nfaces,maxfaces,nnode,numverts(j,3))
       enddo
    enddo
    !
    ncells=ncells+nprizm
    !
    do i=1,nhexa
       l=l+1
       if (iflag_cell(l)==0) cycle
       do j=1,iface(4)
          nsum=0
          flocal=0
          do k=1,numverts(j,4)
             flocal(k)=ndc8(hexf(k,j),i)+base
          enddo
          call insert_face(flocal,face,ipointer,i+ncells,&
               hexmap(j),nfaces,maxfaces,nnode,numverts(j,4))
       enddo
    enddo
    !
    ! find outer boundary face indices
    !
    ntrifaces=0
    nquadfaces=0
    nbface=0
    !
    allocate(cell2cell(4*ntetra+5*npyra+5*nprizm+6*nhexa))
    allocate(faceInfo(nfaces*8))
    !write(6,*) 'ntetra,npyra,nprizm,nhexa,nfaces=',ntetra,npyra,nprizm,nhexa,nfaces
    do i=1,nfaces
       if (face(4,i)==0) then
          faceInfo((i-1)*8+1)=face(3,i)-base
          faceInfo((i-1)*8+2)=face(2,i)-base
          faceInfo((i-1)*8+3)=face(1,i)-base
          faceInfo((i-1)*8+4)=face(1,i)-base
          ntrifaces=ntrifaces+1
       else
          faceInfo((i-1)*8+1)=face(4,i)-base
          faceInfo((i-1)*8+2)=face(3,i)-base
          faceInfo((i-1)*8+3)=face(2,i)-base
          faceInfo((i-1)*8+4)=face(1,i)-base
          nquadfaces=nquadfaces+1
       endif
       faceInfo((i-1)*8+5)=face(5,i)-base
       faceInfo((i-1)*8+6)=face(8,i)-base
       faceInfo((i-1)*8+7)=face(6,i)-base
       faceInfo((i-1)*8+8)=face(9,i)-base
       cell2cell(cell_index(face(5,i))+face(8,i))=face(6,i)-base
       if (face(6,i) > 0) then
          cell2cell(cell_index(face(6,i))+face(9,i))=face(5,i)-base
       else
          nbface=nbface+1
       endif
    enddo
    !
    !write(6,*) 'nbface=',nbface
    ! nbface=0
    ! do i=1,ncells
    !    do j=1,5
    !       if (cell2cell(5*(i-1)+j) .lt. 0) then
    !          nbface=nbface+1
    !       endif
    !    enddo
    ! enddo
    ! write(6,*) 'nbface=',nbface
    !
    deallocate(face)
    deallocate(ipointer)
    deallocate(iflag_cell)
    call cpu_time(t2)
    !
    if (verbose) write(6,*) 'Face calculation time :',t2-t1
    return
  contains
    function cell_index(cin) result (indx)
      integer, intent(in) :: cin
      integer :: indx
      if (cin <= ntetra) then
         indx=4*(cin-1)
      elseif (cin <= ntetra+npyra) then
         indx=5*(cin-(ntetra)-1)+4*ntetra
      elseif (cin <= ntetra+npyra+nprizm) then
         indx=5*(cin-(ntetra+npyra)-1)+4*ntetra+5*npyra
      elseif (cin <= ntetra+npyra+nprizm+nhexa) then
         indx=6*(cin-(ntetra+npyra+nprizm)-1)+4*ntetra+5*npyra+5*nprizm
      endif
    end function cell_index
    !
  end subroutine findnewboundaryface
  !
  !===========================================================================
  subroutine sort_face(a,n)
    implicit none
    integer :: n
    integer :: a(n)
    integer i,j,tmp

    do i=1,n
       do j=i+1,n
          if (a(j) < a(i)) then
             tmp=a(j)
             a(j)=a(i)
             a(i)=tmp
          endif
       enddo
    enddo

    return
  end subroutine sort_face
  !
  ! insert the face into a hashed list
  !
  !===========================================================================
  subroutine insert_face(flocal,face,ipointer,cellIndex,faceIndex, &
       nfaces,maxfaces,nnode,numpts)
    implicit none
    integer :: nfaces,maxfaces,nnode,numpts
    integer :: face(9,maxfaces)
    integer :: cellIndex,faceIndex
    integer :: flocal(4)
    integer :: ipointer(nnode)
    !
    ! local variables
    !
    integer :: ip,ip0
    integer :: f1(4),f2(4)
    !
    !
    f1=flocal
    call sort_face(f1,numpts)
    !
    ip0=ipointer(f1(1))
    ip=ip0
    checkloop: do while(ip > 0)
       f2=face(1:4,ip)
       call sort_face(f2,numpts)
       if (sum(abs(f1-f2))==0) then
          face(6,ip)=cellIndex
          face(9,ip)=faceIndex
          exit checkloop
       endif
       ip=face(7,ip)
    enddo checkloop
    if (ip==0) then
       nfaces=nfaces+1
       face(1:4,nfaces)=flocal
       face(5,nfaces)=cellIndex
       face(6,nfaces)=0
       face(7,nfaces)=ipointer(f1(1))
       face(8,nfaces)=faceIndex
       face(9,nfaces)=0
       ipointer(f1(1))=nfaces
    endif
    return
  end subroutine insert_face

end module boundary_interface

subroutine get_exposed_faces_prizms(ndc6,nprizm)
  use boundary_interface
  implicit none
  integer, intent(in) :: nprizm
  integer, intent(in) :: ndc6(6*nprizm)
  integer, allocatable :: ndc4(:,:),ndc5(:,:),ndc8(:,:),iflag(:)
  integer :: nnode,ntetra,npyra,nhexa

  nnode=maxval(ndc6)+1
  allocate(iflag(nnode))
  iflag=1
  ntetra=0
  npyra=0
  nhexa=0

  call findnewboundaryface(cell2cell,faceInfo,iflag,ntrifaces,nquadfaces,&
                           ndc4,ndc5,ndc6,ndc8, &
                           nnode,ntetra,npyra,nprizm,nhexa)
  deallocate(iflag)

end subroutine get_exposed_faces_prizms

subroutine get_face_count(ntriout,nquadout)
  use boundary_interface
  implicit none
  integer, intent(inout) :: ntriout,nquadout
  ntriout=ntrifaces
  nquadout=nquadfaces
end subroutine get_face_count

subroutine get_graph(cellout,faceout,csize,fsize)
  use boundary_interface
  implicit none
  integer, intent (in) :: fsize,csize
  integer, intent (out):: cellout(csize)
  integer, intent(out) :: faceout(fsize)
  integer :: m,i

  do i=1,fsize
     faceout(i)=faceInfo(i)
  enddo
  do i=1,csize
     cellout(i)=cell2cell(i)
  enddo

  call deleteData

end subroutine get_graph

max_stars_repo_name: jsitaraman/fvsand

max_issues_repo_name: jsitaraman/fvsand

max_forks_repo_name: jsitaraman/fvsand


text:
#if 0
Hello World! program
#endif
program hello
  print *, "Hello World!" ! output
end program hello

! commented line

max_stars_repo_name: idoo/regEx-comment

max_issues_repo_name: idoo/regEx-comment

max_forks_repo_name: idoo/regEx-comment


!WARNING : this code is slightly overbuilt for the actual problem
! but is a handy illustration of types etc

MODULE numerical_rootfinder

  !Implicit none should be in all modules and the main program.

  IMPLICIT NONE

  ! We probably do need Dbl to get these right. This is the
  ! most portable way to get it at the moment
  INTEGER, PARAMETER :: dbl = SELECTED_REAL_KIND(15, 307)

  ! Type to hold polynomial info
  ! I have made mine work for any order
  ! I assume coefficients go largest power to smallest
  TYPE polynomial_wrapper
    ! Since this holds order+1, but it useful to have as the size of coeffs
    INTEGER :: order_plus
    REAL(KIND=dbl), DIMENSION(:), ALLOCATABLE :: coeffs
  END TYPE

  CONTAINS

  SUBROUTINE define_poly(poly, coeffs, err)

    TYPE(polynomial_wrapper), INTENT(INOUT) :: poly
    REAL(KIND=dbl), DIMENSION(:), INTENT(IN) :: coeffs
    ! I make this optional because maybe I have already checked validity of coeffs
    LOGICAL, OPTIONAL :: err
    LOGICAL :: err_in

    err_in = .FALSE.

    poly%order_plus = SIZE(coeffs)

    ! 0th order poly has no roots, and I arbitrarily limit to max 8
    IF(poly%order_plus < 1 .OR. poly%order_plus > 8) THEN
      PRINT*, "Polynomial order should be between 1 and 8"
      err_in = .TRUE.
    ELSE
      ! In case we're re-defining
      IF(ALLOCATED(poly%coeffs)) DEALLOCATE(poly%coeffs)
      ALLOCATE(poly%coeffs(poly%order_plus))
      poly%coeffs = coeffs
    END IF

    ! Put value in before returning
    IF(PRESENT(err)) err = err_in

  END SUBROUTINE

  ! f is a terrible function name in general, but it matches the equation I am coding
  FUNCTION f(poly, x)
    TYPE(polynomial_wrapper) :: poly
    REAL(KIND=dbl) :: f, x, x_term
    INTEGER :: i

    f = 0.0
    x_term = 1.0

    ! I go backwards so I can accumulate the x powers
    DO i = poly%order_plus, 1,  -1
      f = f + x_term * poly%coeffs(i)
      x_term = x_term * x
    END DO

  END FUNCTION f

  FUNCTION f_prime(poly, x)
    TYPE(polynomial_wrapper) :: poly
    REAL(KIND=dbl) :: f_prime, x, x_term
    INTEGER :: i

    f_prime = 0.0
    x_term = 1.0

    ! I go backwards so I can accumulate the x powers
    DO i = poly%order_plus-1, 1,  -1
      f_prime = f_prime + x_term * poly%coeffs(i) * REAL(poly%order_plus-i)
      x_term = x_term * x
    END DO

  END FUNCTION f_prime

  FUNCTION root_near(poly, guess)
    ! Find a root of a function near a given guess
    ! We can return any valid root

    ! The polynomial
    TYPE(polynomial_wrapper) :: poly
    ! Initial guess for the root
    REAL(KIND=dbl), INTENT(IN) :: guess

    REAL(KIND=dbl):: change
    REAL(KIND=dbl) :: root_near
    REAL(KIND=dbl), PARAMETER :: threshold = 0.001_dbl
    INTEGER, PARAMETER :: max_iter = 50
    INTEGER :: i

    ! See numerical_root for discussion of termination condition

    root_near = guess
    ! Start with a number for change that definitely
    ! does not trip the fractional change limit
    change = guess * 100.0
    DO i = 1, max_iter
      !I chose to print the progress in this
      PRINT*, "Iteration", i, "Fractional change", change, "Guess", root_near
      ! I am using a fractional change limit
      IF(ABS(change/root_near) .LT. threshold) EXIT

      ! I am being a bit "clever" to minimise my temporaries here
      ! by storing the change before applying it
      ! I could also store the old and new values, or other
      ! approaches
      change = f(poly,root_near)/f_prime(poly,root_near)
      root_near = root_near - change
    END DO

    ! If we hit max_iter, we probably don't have a root
    ! This is a dumb way to create a NaN - we should really
    ! use the IEEE module, but I am keeping this simple
    IF(i >= max_iter) root_near = 0.0*1.0/(root_near-root_near)

  END FUNCTION root_near

END MODULE numerical_rootfinder

PROGRAM main

  USE numerical_rootfinder

  IMPLICIT NONE

  ! Select some bounds for finding
  REAL(KIND=dbl) :: guess = 0.555555555_dbl
  REAL(KIND=dbl) :: res
  TYPE(polynomial_wrapper) :: poly

  CALL define_poly(poly, (/1.0_dbl, -0.31428571428_dbl, &
      -0.21549295774_dbl, 0.20845070422_dbl/))

  ! Calculate and print the answer
  ! The Fortran standard lets compilers do something unhelpful
  ! if you PRINT the result of a function which has a print in,
  ! so I generally always capture a result, then print it
  res = root_near(poly, guess)

  PRINT*, "Root is ", res

END PROGRAM


text:
!	Scale function
!	dimtrans is the dimension of the transformation.
!	theta is the argument.
function transform(dimtrans,theta)
use transset
implicit none
integer,intent(in)::dimtrans
real (kind=8),intent(in)::theta(:)
real(kind=8)::transform(dimtrans)
!	elementin is the source element.
!	elementout is the output element.
!	Component counter is i.
!	j counts degrees.
!	k and kk are used to find coefficients of polynomials.
integer::elementin,elementout,i,j,k,kk
!	bas is the base.
!	pol is the polynomial value.
!	th is the element of theta used in the component.
!	wt is the component weight.
real(kind=8)::bas,pol,th,wt
transform=0.0_8
kk=1

do i=1,components
	k=kk
	kk=kk+degrees(i)+1

	elementout=elementsout(i)
	th=dot_product(inweights(:,i),theta)
	wt=1.0_8
	!	See if th is within bounds.
	if(lbmasks(i)) then
		if(lbs(i)>th) cycle
	end if
	if(ubmasks(i))then
		if(ubs(i)<th)cycle
	end if
	if(th==lbs(i).or.th==ubs(i))wt=0.5_8
!	Polynomial computation.
	pol=0.0_8
	bas=bases(i)
	do j=0,degrees(i)
		pol=pol*(th-bas)+coeff(kk-j-1)
	end do
!	Add component.
	transform(elementout)=transform(elementout)+wt*pol
end do

end function transform

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


!
!     ------------------------------------------------------------------
!              H L
!     ------------------------------------------------------------------
!
!       Returns the value of <i|L|j>, using a special formula to
!  preserve symmetry.
!
      DOUBLE PRECISION FUNCTION HL (EL, I, J, REL)
      use param_C
      USE radial_C
      USE nel_C, EL_=>EL
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
!      PARAMETER(NOD=220)
      PARAMETER (IWRITE = 6)
!...Translated by Pacific-Sierra Research 77to90  4.3E  00:33:16  11/16/01
!...Switches:
      CHARACTER*1 :: EL(*)*3
!      COMMON/PARAM/H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,NCFG,IB,IC,ID
!     :   ,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD,RMASS
!      COMMON/RADIAL/R(NOD),RR(NOD),R2(NOD),YK(NOD),YR(NOD),X(NOD)
!      POINTER(IQP,P(NOD,1)),(IQN,N(1)),(IQL,L(1)),(IQAZ,AZ(1)),
!     :       (IQMAX,MAX_(1))
!      COMMON/NEL/IQP,IQN,IQL,IQAZ,IQMAX,IQ(7)
      DIMENSION :: DI(NOD,6), DJ(NOD,6), G(NOD)
      LOGICAL :: REL
!
      IF (IABS(L(I)-L(J)) .NE. 0) THEN
         WRITE (IWRITE, 4) EL(I), L(I), EL(J), L(J)
    4    FORMAT(10X,'UNALLOWED L VALUES OCCURRED IN HL SUBROUTINE'/,2(10X,A3,&
            ' HAS L = ',I3))
      ENDIF
      LI = L(I)
      C = 2*LI + 1
      A1 = -D2/(C*(LI + 1))
      A2 = A1/((C + D2)*(LI + 1))
      A3 = A2/((LI + 2)*(LI + 1))
      ZR = Z*R(1)
      HL = H*C*P(1,I)*P(1,J)*(D1 + ZR*(A1 + ZR*(A2 + ZR*A3)))
      MM = MIN0(MAX_(I) + 3,MAX_(J) + 3,ND - 1)
      K = 2
      C = D4/D3
      DI(K,1) = P(K+1,I) - P(K-1,I)
      DI(K,2) = P(K+1,I) - D2*P(K,I) + P(K-1,I)
      DJ(K,1) = P(K+1,J) - P(K-1,J)
      DJ(K,2) = P(K+1,J) - D2*P(K,J) + P(K-1,J)
      HL = HL + DI(K,1)*DJ(K,1) + C*DI(K,2)*DJ(K,2)
      DO K = 4, MM, 2
         DI(K,1) = P(K+1,I) - P(K-1,I)
         DI(K,2) = P(K+1,I) - D2*P(K,I) + P(K-1,I)
         DI(K,4) = P(K+2,I) - D4*(P(K+1,I)+P(K-1,I)) + D6*P(K,I) + P(K-2,I)
         DI(K,3) = P(K+2,I) - P(K-2,I) - D2*DI(K,1)
         DI(K,5) = P(K+3,I) - P(K-3,I) - D4*(P(K+2,I)-P(K-2,I)) + 5.D0*(P(K+1,I&
            )-P(K-1,I))
         DI(K,6) = P(K+3,I) + P(K-3,I) - D6*(P(K+2,I)+P(K-2,I)) + 15.D0*(P(K+1,&
            I)+P(K-1,I)) - 20.D0*P(K,I)
         DJ(K,1) = P(K+1,J) - P(K-1,J)
         DJ(K,2) = P(K+1,J) - D2*P(K,J) + P(K-1,J)
         DJ(K,4) = P(K+2,J) - D4*(P(K+1,J)+P(K-1,J)) + D6*P(K,J) + P(K-2,J)
         DJ(K,3) = P(K+2,J) - P(K-2,J) - D2*DJ(K,1)
         DJ(K,5) = P(K+3,J) - P(K-3,J) - D4*(P(K+2,J)-P(K-2,J)) + 5.D0*(P(K+1,J&
            )-P(K-1,J))
         DJ(K,6) = P(K+3,J) + P(K-3,J) - D6*(P(K+2,J)+P(K-2,J)) + 15.D0*(P(K+1,&
            J)+P(K-1,J)) - 20.D0*P(K,J)
      END DO
      DO K = 4, MM, 2
         G(K) = DI(K,1)*DJ(K,1) + C*DI(K,2)*DJ(K,2) + (DI(K,3)*DJ(K,3)+DI(K,2)*&
            DJ(K,4)+DI(K,4)*DJ(K,2))/45.D0 - (DI(K,3)*DJ(K,5)+DI(K,5)*DJ(K,3))/&
            252.D0 - (DI(K,2)*DJ(K,6)+DI(K,6)*DJ(K,2)-1.1*DI(K,4)*DJ(K,4))/&
            378.D0
      END DO
      DO K = 4, MM, 2
         HL = HL + G(K)
      END DO
      TZ = Z + Z
      C = (LI + D5)**2
      HL2 = D5*(TZ*R(1)-C)*P(1,I)*P(1,J)
      DO K = 2, MM, 2
         G(K) = D2*(TZ*R(K)-C)*P(K,I)*P(K,J) + (TZ*R(K+1)-C)*P(K+1,I)*P(K+1,J)
      END DO
      DO K = 2, MM, 2
         HL2 = HL2 + G(K)
      END DO
      HL = (-HL/(D2*H)) + HL2*H1
      IF (REL) HL = HL - D2*RLSHFT(I,J)
      RETURN
      END FUNCTION HL


text:
! { dg-do compile }
! PR 24643
! substring references on implicitly typed CHARACTER variables didn't work
        PROGRAM P
        IMPLICIT CHARACTER*8 (Y)
        YLOCAL='A'
        YBTABLE=YLOCAL(1:2)
        END

max_stars_repo_name: vidkidz/crossbridge

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
address(Howard IM Field, 38.545444, 121.749643) is one of the intramural fields on the UC Davis campus. It is located on the east side of Howard Way just north of the North Entry Parking Structure. The field is bounded on the east by Toomey Field and on the north by Russell Boulevard.

Howard IM Field is equipped with lights. One of the light towers was replaced in 2007 to include a cell phone towers cell phone tower for Verizon.

During the fall, the UC Davis football UCD football team practices regularly on this field. During football season, it is usually lined with two football fields. At other times of the year, the field is sometimes lined for other field sports. Occasionally during football season the California Aggie Marching Banduh Banduh practices here when their other practice sites are unavailable. Regardless of when said practice takes place, it is at that exact time that the sprinklers generally decide to turn on.

When football games were still held at Toomey Field, the associated pregame festivities typically took place here. The Aggie Pack entrance to Toomey Field also faced Howard IM Field.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Describe Users/ChristineHsueh here.

20130128 10:05:21 nbsp Hi CHsueh!  Welcome to the Wiki.  I still havent decided whether Im going to do the run, but it sounds fun. Users/TomGarberson


text:
C     Checking C prettyprint of Fortran code: how about labelless
C     continue?
      subroutine continue_test
      logical l1, l2, l3, l4

      read *,x
      if(x.gt.0.) then
         continue
      elseif(x.lt.-1.) then
         continue
      else
c        x=x+1
      endif

      end


text:
!===============================================================================
function pressure(depth) result(press)

! !DESCRIPTION:
!  This function computes pressure in bars from depth in meters
!  using a mean density derived from depth-dependent global
!  average temperatures and salinities from Levitus 1994, and
!  integrating using hydrostatic balance.
!
!  References:
!
!     Levitus, S., R. Burgett, and T.P. Boyer, World Ocean Atlas
!          1994, Volume 3: Salinity, NOAA Atlas NESDIS 3, US Dept. of
!          Commerce, 1994.
!
!     Levitus, S. and T.P. Boyer, World Ocean Atlas 1994,
!          Volume 4: Temperature, NOAA Atlas NESDIS 4, US Dept. of
!          Commerce, 1994.
!
!     Dukowicz, J. K., 2000: Reduction of Pressure and Pressure
!          Gradient Errors in Ocean Simulations, J. Phys. Oceanogr.,
!          submitted.
!
! !REVISION HISTORY:
!  same as module

! !INPUT PARAMETERS:

real, intent(in) :: depth    ! depth in meters

! !OUTPUT PARAMETERS:

real :: press   ! pressure [bars]

! !LOCAL PARAMETERS:
!  adjustment from original code which was written as 0.1_r8
!  which the compiler disliked so constants are now pre-defined here
double precision, parameter ::              &
   pc1 = 0.059808,   &
   pc2 = 0.025,      &
   pc3 = 0.100766,   &
   pc4 = 2.28405e-7, &
   c1  = 1.0
!EOP
!BOC
!-----------------------------------------------------------------------
!  convert depth in meters to pressure in bars
!-----------------------------------------------------------------------

press = pc1*(exp(-pc2*depth) - c1) + pc3*depth + pc4*depth**2

!return

end function pressure
!===============================================================================

!===============================================================================
subroutine state(S_orig, T, P_orig, RHO)

!  ADJUSTED/SLIMMED DOWN FROM state_mod.F90 June 2016

! !DESCRIPTION:
!  Returns the density of water at level k from equation of state
!  $\rho = \rho(d,\theta,S)$ where $d$ is depth, $\theta$ is
!  potential temperature, and $S$ is salinity. the density can be
!  returned as a perturbation (RHOOUT) or as the full density
!  (RHOFULL). Note that only the polynomial EOS choice will return
!  a perturbation density; in other cases the full density is returned
!  regardless of which argument is requested.
!
!  This routine also computes derivatives of density with respect
!  to temperature and salinity at level k from equation of state
!  if requested (ie the optional arguments are present).
!
!  If $k = kk$ are equal the density for level k is returned.
!  If $k \neq kk$ the density returned is that for a parcel
!  adiabatically displaced from level k to level kk.
!
! !REVISION HISTORY:
!  same as module

double precision, parameter :: rho0=4.1/3.996*1000

! !INPUT PARAMETERS:

double precision, intent(in) :: &
  T,               &! pot. temperature at level k  [degC]
  S_orig,          &! salinity at level k         [kg/kg]
  P_orig            ! pressure at reference level   [bar]

! !OUTPUT PARAMETERS:

double precision, intent(out) :: &
  RHO               ! density of water [kg m^-3]

!-----------------------------------------------------------------------
!  local variables:
!-----------------------------------------------------------------------

double precision :: &
  S,                &! [PSU]
  P,                &! [dbar]
  SQR,DENOMK,       &! work arrays
  WORK1, WORK2, T2

double precision, parameter :: &
  c10   = 1.0E01,              &
  c1000 = 1.0E03

!-----------------------------------------------------------------------
!  first check for valid range if requested
!-----------------------------------------------------------------------

! THIS IS SKIPPED HERE

!-----------------------------------------------------------------------
!  MWJF EOS coefficients
!-----------------------------------------------------------------------

!*** these constants will be used to construct the numerator

double precision, parameter ::   &
  mwjfnp0s0t0 =   9.99843699e+2, &
  mwjfnp0s0t1 =   7.35212840e+0, &
  mwjfnp0s0t2 =  -5.45928211e-2, &
  mwjfnp0s0t3 =   3.98476704e-4, &
  mwjfnp0s1t0 =   2.96938239e+0, &
  mwjfnp0s1t1 =  -7.23268813e-3, &
  mwjfnp0s2t0 =   2.12382341e-3, &
  mwjfnp1s0t0 =   1.04004591e-2, &
  mwjfnp1s0t2 =   1.03970529e-7, &
  mwjfnp1s1t0 =   5.18761880e-6, &
  mwjfnp2s0t0 =  -3.24041825e-8, &
  mwjfnp2s0t2 =  -1.23869360e-11

!*** these constants will be used to construct the denominator

double precision, parameter ::    &
  mwjfdp0s0t0 =   1.0e+0,         &
  mwjfdp0s0t1 =   7.28606739e-3,  &
  mwjfdp0s0t2 =  -4.60835542e-5,  &
  mwjfdp0s0t3 =   3.68390573e-7,  &
  mwjfdp0s0t4 =   1.80809186e-10, &
  mwjfdp0s1t0 =   2.14691708e-3,  &
  mwjfdp0s1t1 =  -9.27062484e-6,  &
  mwjfdp0s1t3 =  -1.78343643e-10, &
  mwjfdp0sqt0 =   4.76534122e-6,  &
  mwjfdp0sqt2 =   1.63410736e-9,  &
  mwjfdp1s0t0 =   5.30848875e-6,  &
  mwjfdp2s0t3 =  -3.03175128e-16, &
  mwjfdp3s0t1 =  -1.27934137e-17

!*** MWJF numerator coefficients including pressure

double precision ::&
  mwjfnums0t0, mwjfnums0t1, mwjfnums0t2, mwjfnums0t3,                          &
  mwjfnums1t0, mwjfnums1t1, mwjfnums2t0,                                       &
  mwjfdens0t0, mwjfdens0t1, mwjfdens0t2, mwjfdens0t3, mwjfdens0t4,             &
  mwjfdens1t0, mwjfdens1t1, mwjfdens1t3,                                       &
  mwjfdensqt0, mwjfdensqt2
!-----------------------------------------------------------------------
!  McDougall, Wright, Jackett, and Feistel EOS
!  test value : rho = 1.033213242 for
!  S = 35.0 PSU, theta = 20.0, pressz = 200.0
!
!  from paper:
!  [ (PSU, degC, dbar) -> kg m^-3 ]
!  (35, 25, 2000) -> 1031.654 229
!  ( 0, 20, 1000) -> 1017.726 743
!  (40, 12, 8000) -> 1062.928 258
!-----------------------------------------------------------------------

!   unit conversion
P   = c10*P_orig   ! [bar] -> [dbar]
S   = c1000*S_orig ! [kg/kg] -> [PSU]
SQR = sqrt(S)      ! square root

!      write (*,*) 'S =', S, 'PSU'
!      write (*,*) 'T =', T, 'degC'
!      write (*,*) 'P =', P, 'dbar'

!*** first calculate numerator of MWJF density [P_1(S,T,p)]
mwjfnums0t0 = mwjfnp0s0t0 + P*(mwjfnp1s0t0 + P*mwjfnp2s0t0)
mwjfnums0t1 = mwjfnp0s0t1
mwjfnums0t2 = mwjfnp0s0t2 + P*(mwjfnp1s0t2 + P*mwjfnp2s0t2)
mwjfnums0t3 = mwjfnp0s0t3
mwjfnums1t0 = mwjfnp0s1t0 + P*mwjfnp1s1t0
mwjfnums1t1 = mwjfnp0s1t1
mwjfnums2t0 = mwjfnp0s2t0

WORK1 = mwjfnums0t0 + T * (mwjfnums0t1 + T * (mwjfnums0t2 +                    &
        mwjfnums0t3 * T)) + S * (mwjfnums1t0 +                                 &
        mwjfnums1t1 * T + mwjfnums2t0 * S)

!*** now calculate denominator of MWJF density [P_2(S,T,p)]
mwjfdens0t0 = mwjfdp0s0t0 + P*mwjfdp1s0t0
mwjfdens0t1 = mwjfdp0s0t1 + P**3 * mwjfdp3s0t1
mwjfdens0t2 = mwjfdp0s0t2
mwjfdens0t3 = mwjfdp0s0t3 + P**2 * mwjfdp2s0t3
mwjfdens0t4 = mwjfdp0s0t4
mwjfdens1t0 = mwjfdp0s1t0
mwjfdens1t1 = mwjfdp0s1t1
mwjfdens1t3 = mwjfdp0s1t3
mwjfdensqt0 = mwjfdp0sqt0
mwjfdensqt2 = mwjfdp0sqt2

WORK2 = mwjfdens0t0 + T * (mwjfdens0t1 + T * (mwjfdens0t2 +                    &
        T   * (mwjfdens0t3 + mwjfdens0t4 * T))) +                              &
        S   * (mwjfdens1t0 + T * (mwjfdens1t1 + T*T*mwjfdens1t3)+              &
        SQR * (mwjfdensqt0 + T*T*mwjfdensqt2))

RHO   = WORK1/WORK2-rho0

! SKIPPED RHOFULL, DRHODT, DRHODS CALCULATIONS HERE
! SKIPPED OTHER EOS FUNCTIONS

end subroutine state
!===============================================================================

subroutine hydrostatic_pressure(g,dz,RHO,SSH,HSP)
implicit none
!
!  calculate hydrostatic pressure with in-situ density
!
integer                                     :: k
integer, parameter                          :: imt=3600, jmt=2400, km=42
double precision,               intent(in)  :: g
real,    dimension(km),         intent(in)  :: dz
real,    dimension(imt,jmt,km), intent(in)  :: RHO
real,    dimension(imt,jmt),    intent(in)  :: SSH
real,    dimension(imt,jmt,km), intent(out) :: HSP
real,    parameter                          :: p5 = 0.5

do k=1,km
  if (k==1) then
    HSP(:,:,k) = g*(RHO(:,:,k)*1.0E3)*(dz(k)/2.0 + SSH(:,:)*1.0E-02)
  else
    HSP(:,:,k) = HSP(:,:,k-1) + p5*g*((RHO(:,:,k-1)+RHO(:,:,k))*1.0E03)&
                                 *(dz(k-1)+dz(k))/2.0
  endif
enddo !k
end subroutine hydrostatic_pressure

max_stars_repo_name: AJueling/LEC

max_issues_repo_name: AJueling/LEC

max_forks_repo_name: AJueling/LEC


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!!   MIT License
!!!!
!!!!   ParaMonte: plain powerful parallel Monte Carlo library.
!!!!
!!!!   Copyright (C) 2012-present, The Computational Data Science Lab
!!!!
!!!!   This file is part of the ParaMonte library.
!!!!
!!!!   Permission is hereby granted, free of charge, to any person obtaining a
!!!!   copy of this software and associated documentation files (the "Software"),
!!!!   to deal in the Software without restriction, including without limitation
!!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!!!!   and/or sell copies of the Software, and to permit persons to whom the
!!!!   Software is furnished to do so, subject to the following conditions:
!!!!
!!!!   The above copyright notice and this permission notice shall be
!!!!   included in all copies or substantial portions of the Software.
!!!!
!!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
!!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
!!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
!!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
!!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
!!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
!!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
!!!!
!!!!   ACKNOWLEDGMENT
!!!!
!!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
!!!!   As per the ParaMonte library license agreement terms, if you use any parts of
!!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
!!!!   work (education/research/industry/development/...) by citing the ParaMonte
!!!!   library as described on this page:
!!!!
!!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!> \brief
!> This module contains procedures and routines for the Kmeans the statistical classification and clustering methods.
!> \author
!> Amir Shahmoradi

module KmeansOOP_mod

    use Constants_mod, only: IK, RK
    use Err_mod, only: Err_type

    implicit none

    character(*), parameter :: MODULE_NAME = "@KmeansOOP_mod"

    real(RK) :: ndKuniform, btFactor

    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> The [KmeansTry_type](@ref kmeanstry_type) type containing information about one instance of Kmeans clustering.
    type :: KmeansTry_type
        integer(IK) , allocatable   :: nzsci                !< The number of Zero-Sized Clusters Iterations at each Kmeans try.
        integer(IK) , allocatable   :: niter                !< The number of iterations to reach convergence at each Kmeans try.
        integer(IK) , allocatable   :: Size(:)              !< An array of size `(ncMax)` containing the sizes of the individual clusters identified.
        real(RK)    , allocatable   :: potential            !< The sum of the distances-squared of all points from their corresponding group centers.
        real(RK)    , allocatable   :: Center(:,:)          !< An array of size `(ndMax,ncMax)` containing the centers of the individual clusters identified.
        real(RK)    , allocatable   :: NormedPoint(:,:,:)   !< An array of size `(ndMax,npMax,ncMax)` containing the normalized `Point(nd,np)` w.r.t. each cluster.
        real(RK)    , allocatable   :: MinDistanceSq(:)     !< An array of size `(npMax)` containing the minimum distances of the points from the cluster centers.
        integer(IK) , allocatable   :: Membership(:)        !< An array of size `(npMax)` each element of which represents the cluster ID to which the point belongs.
        type(Err_type)              :: Err                  !< An object of class [Err_type](@ref err_mod::err_type) containing information about error occurrence.
    end type KmeansTry_type

    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> The `KmeansOOP_type` class.
    type :: KmeansOOP_type
        logical                             :: failed               !< A logical value. If `.true.`, then none of the Kmeans clustering tries have succeeded.
        real(RK)                            :: relTol = 1.e-4_RK    !< The current requested relative tolerance below which should the potential fall, convergence is assumed.
        integer(IK)                         :: nd = 0_IK            !< The current number of data attributes to be used in clustering. The condition `nd<=ndMax` must hold at all times.
        integer(IK)                         :: np = 0_IK            !< The current number of input data points to be clustered. The condition `np<=npMax` must hold at all times.
        integer(IK)                         :: nc = 0_IK            !< The current number of Kmeans clusters to be requested. The condition `nc<=ncMax` must hold at all times.
        integer(IK)                         :: nt = 0_IK            !< The current number of Kmeans clustering trials. The condition `nt<=ntMax` must hold at all times.
        integer(IK)                         :: itopt = 0_IK         !< The Kmeans clustering try which has yielded the least potential.
        integer(IK)                         :: ndMax = 0_IK         !< The maximum possible number of data attributes to be used in clustering. Used for array allocation upon object instantiation.
        integer(IK)                         :: npMax = 0_IK         !< The maximum possible number of input data points to be clustered. Used for array allocation upon object instantiation.
        integer(IK)                         :: ncMax = 0_IK         !< The maximum possible number of Kmeans clusters to be requested. Used for array allocation upon object instantiation.
        integer(IK)                         :: ntMax = 0_IK         !< The maximum possible number of Kmeans clustering trials. Used for array allocation upon object instantiation.
        integer(IK)                         :: niterMax = 1000_IK   !< The maximum number of iterations beyond which it a lack of converge is assumed.
        integer(IK)                         :: nzsciMax = 10_IK     !< The maximum number of zero-sized cluster iterations beyond which it a lack of converge is assumed.
        type(KmeansTry_type), allocatable   :: Try(:)               !< An array of size `(ntMax)` of object of type [KmeansTry_type](@ref kmeanstry_type).
    contains
        procedure, pass :: run => runKmeans
        procedure, pass, private :: allocate => allocateKmeansComponents
    end type KmeansOOP_type

    interface KmeansOOP_type
        module procedure :: constructKmeans
    end interface KmeansOOP_type

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

contains

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Construct and return an object of [KmeansOOP_type](@ref KmeansOOP_type) class.
    !> Note that this constructor solely allocates the components of the output
    !> `Kmeans` object to be later used for Kmeans clustering.
    !>
    !> \param[in]       ndMax       :   The maximum possible number of dimensions (features) of the input data array `Point`.
    !> \param[in]       npMax       :   The maximum possible number of observations in the input data array `Point`.
    !> \param[in]       ncMax       :   The maximum possible number of clusters to be found via the Kmeans algorithm.
    !> \param[in]       ntMax       :   The maximum possible number of observations in the input data array `Point`.
    !> \param[in]       niterMax    :   The maximum number of iterations beyond which it a lack of converge is assumed (**optional**, default = 1000).
    !> \param[in]       nszciMax    :   The maximum number of zero-sized cluster iterations beyond which it a lack of converge is assumed (**optional**, default = 10).
    !> \param[in]       relTol      :   The current requested relative tolerance below which should the potential fall, convergence is assumed (**optional**, default = 1.e-8).
    !>
    !> \return
    !> `Kmeans` :   An object of type [KmeansCluster_type](@ref kmeanscluster_type) containing the properties of the identified clusters.
    !>              The `Err%stat` component of the output `Kmeans` will be set to `1` if the number of iterations reaches the input
    !>              `niterMax` or its default value before convergence occurs.
    !>              The `Err%stat` component of the output `Kmeans` will be set to `2` if more than `nzsciMax` clustering tries
    !>              result in clusters with zero members.
    !>              In both cases in the above, `Err%occurred = .true.` upon exiting the procedure.
    !>
    !> \warning
    !> The input value for `nt` must be at least 1, otherwise it will lead to a runtime error or at best, the output will be undefined.
    !>
    !> \warning
    !> The input value for `niterMax` must be larger than 0 since at least one iteration will be performed by the procedure.
    !>
    !> \warning
    !> This algorithm does not check for the consistency of the input values for `nc`, `nd`, and `np`.
    !> It is currently the user's responsibility to ensure that these input values are logical and sensible.
    !>
    !> \author
    ! Amir Shahmoradi, April 03, 2017, 2:16 PM, ICES, UTEXAS
    function constructKmeans(ndMax, npMax, ncMax, ntMax, niterMax, nzsciMax, relTol) result(Kmeans)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: constructKmeans
#endif
        use Constants_mod, only: IK, RK
        implicit none
        type(KmeansOOP_type)                :: Kmeans
        integer(IK) , intent(in)            :: ndMax, npMax, ncMax, ntMax
        integer(IK) , intent(in), optional  :: niterMax, nzsciMax
        real(RK)    , intent(in), optional  :: relTol
        integer(IK)                         :: it
        if (present(relTol)) Kmeans%relTol = relTol
        if (present(niterMax)) Kmeans%niterMax = niterMax
        if (present(nzsciMax)) Kmeans%nzsciMax = nzsciMax
        call Kmeans%allocate(ndMax, npMax, ncMax, ntMax)
    end function constructKmeans

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    subroutine allocateKmeansComponents(Kmeans, ndMax, npMax, ncMax, ntMax)
        use Constants_mod, only: IK, RK
        implicit none
        class(KmeansOOP_type), intent(inout)    :: Kmeans
        integer(IK) , intent(in)                :: ndMax, npMax, ncMax, ntMax
        integer(IK)                             :: it
        Kmeans%nd = ndMax
        Kmeans%np = npMax
        Kmeans%nc = ncMax
        Kmeans%nt = ntMax
        Kmeans%ndMax = ndMax
        Kmeans%npMax = npMax
        Kmeans%ncMax = ncMax
        Kmeans%ntMax = ntMax
        if (allocated(Kmeans%Try)) deallocate(Kmeans%Try)
        allocate(Kmeans%Try(ntMax))
        do it = 1, Kmeans%ntMax
            allocate(Kmeans%Try(it)%Size(ncMax), source = 0_IK)
            allocate(Kmeans%Try(it)%Membership(npMax), source = 0_IK)
            allocate(Kmeans%Try(it)%Center(ndMax,ncMax), source = 0._RK)
            allocate(Kmeans%Try(it)%MinDistanceSq(npMax), source = 0._RK)
            allocate(Kmeans%Try(it)%NormedPoint(ndMax,npMax,ncMax), source = 0._RK)
        end do
    end subroutine allocateKmeansComponents

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> This procedure is a method of the class [KmeansOOP_type](@ref KmeansOOP_type).
    !> Perform the Kmeans clustering on the input data set represented by the array `Point(nd, np)`.
    !> Upon performing all iterations of clustering, the index `Kmeans%itopt` will point to the final optimal
    !> iteration of Kmeans clustering which resulted in the smallest potential. If only one clustering try has been requested,
    !> then obviously, `Kmeans%itopt = 1`.
    !>
    !> \param[inout]    Kmeans                  :   An object of class [KmeansOOP_type](@ref KmeansOOP_type) which will contain the clustering results.
    !> \param[in]       Point                   :   The input array of size `(nd,np)` of the points to be clustered via the Kmeans algorithm.
    !> \param[in]       nd                      :   The current number of data attributes to be used in clustering. The condition `nd<=ndMax` must hold at all times.
    !> \param[in]       np                      :   The current number of input data points to be clustered. The condition `np<=npMax` must hold at all times.
    !> \param[in]       nc                      :   The current number of Kmeans clusters to be requested. The condition `nc<=ncMax` must hold at all times.
    !> \param[in]       nt                      :   The current number of Kmeans clustering trials. The condition `nt<=ntMax` must hold at all times.
    !> \param[in]       niterMax                :   The maximum number of iterations beyond which it a lack of converge is assumed (**optional**).
    !> \param[in]       nszciMax                :   The maximum number of zero-sized cluster iterations beyond which it a lack of converge is assumed (**optional**).
    !> \param[in]       relTol                  :   The current requested relative tolerance below which should the potential fall, convergence is assumed (**optional**).
    !> \param[in]       InitCenter              :   An input array of size `(nd,nc,nt)` representing the initial centers of the Kmeans clusters (**optional**).
    !>                                              If not provided, the cluster centers will be initialized via the Kmeans++ algorithm.
    !>
    !> \return
    !> `Kmeans` :   An object of type [KmeansCluster_type](@ref kmeanscluster_type) containing the properties of the identified clusters.
    !>              The `Err%stat` component of the output `Kmeans` will be set to `1` if the number of iterations reaches the input
    !>              `niterMax` or its default value before convergence occurs.
    !>              The `Err%stat` component of the output `Kmeans` will be set to `2` if more than `nzsciMax` clustering tries
    !>              result in clusters with zero members.
    !>              In both cases in the above, `Err%occurred = .true.` upon exiting the procedure.
    !>
    !> \warning
    !> The input value for `niterMax` must be larger than 0 since at least one iteration will be performed by the procedure.
    !>
    !> \warning
    !> This algorithm does not check for the consistency of the input values for `nc`, `nd`, and `np`.
    !> It is currently the user's responsibility to ensure that these input values are logical and sensible.
    !>
    !> \author
    ! Amir Shahmoradi, April 03, 2017, 2:16 PM, ICES, UTEXAS
    subroutine runKmeans(Kmeans, Point, nd, np, nc, nt, niterMax, nzsciMax, relTol, InitCenter)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: runKmeans
#endif
        use Constants_mod, only: IK, RK, NEGINF_IK

        implicit none

        real(RK)    , intent(in)            :: Point(nd,np)
        integer(IK) , intent(in)            :: nd, np, nc, nt
        integer(IK) , intent(in), optional  :: niterMax, nzsciMax
        real(RK)    , intent(in), optional  :: relTol, InitCenter(nd,nc,nt)
        class(KmeansOOP_type)               :: Kmeans

        character(*), parameter :: PROCEDURE_NAME = MODULE_NAME//"@runKmeans()"

        real(RK)                :: potential                ! The sum of all minimum distances squared.
        real(RK)                :: distanceSq               ! distance of each point ip from a given cluster center.
        real(RK)                :: minPotential             ! The minimum of sum of all minimum distances squared across all Kmeans tries.
        real(RK)                :: SumPoint(nd,nc)          ! The sum of the points in each cluster.
        real(RK)                :: MinDistanceSq(np)        ! minimum distance of each point ip from all centers of the clusters.
        integer(IK)             :: Membership(np)           ! current (old) Memberships of the clusters.

        logical                 :: isPresentInitCenter
        logical                 :: convergenceOccurred
        integer                 :: i, j, ip, ic, it

        Kmeans%itopt = 0_IK
        Kmeans%failed = .true.

        if (present(relTol)) Kmeans%relTol = relTol
        if (present(niterMax)) Kmeans%niterMax = niterMax
        if (present(nzsciMax)) Kmeans%nzsciMax = nzsciMax

        if (nd > Kmeans%ndMax .or. np > Kmeans%npMax .or. nc > Kmeans%ncMax .or. nt > Kmeans%ntMax) then
            call Kmeans%allocate(nd, np, nc, nt)
        else
            Kmeans%nd = nd
            Kmeans%np = np
            Kmeans%nc = nc
            Kmeans%nt = nt
        end if

        minPotential = huge(minPotential)
        isPresentInitCenter = present(InitCenter)

        loopMultipleKmeansTry: do it = 1, Kmeans%nt

            Kmeans%Try(it)%Err%occurred = .false.
            Kmeans%Try(it)%Err%stat = NEGINF_IK

            !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            ! Find the distances of all points to all cluster centers. If not given, initialize the centers via the Kmeans++ algorithm.
            !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            if(isPresentInitCenter) then
                !%%%%%%%%%%%%%%%%%%%
                ic = 1
                do ip = 1, np
                    MinDistanceSq(ip) = sum( (InitCenter(1:nd,ic,it) - Point(1:nd,ip))**2 )
                    Membership(ip) = ic
                end do
                !%%%%%%%%%%%%%%%%%%%
                do ic = 2, nc - 1
                    do ip = 1, np
                        distanceSq = sum( (InitCenter(1:nd,ic,it) - Point(1:nd,ip))**2 )
                        if (distanceSq < MinDistanceSq(ip)) then
                            MinDistanceSq(ip) = distanceSq
                            Membership(ip) = ic
                        end if
                    end do
                end do
                !%%%%%%%%%%%%%%%%%%%
                ic = nc
                SumPoint = 0._RK
                potential = 0._RK
                do ip = 1, np
                    distanceSq = sum( (InitCenter(1:nd,ic,it) - Point(1:nd,ip))**2 )
                    if (distanceSq < MinDistanceSq(ip)) then
                        MinDistanceSq(ip) = distanceSq
                        Membership(ip) = ic
                    end if
                    SumPoint(1:nd,Membership(ip)) = SumPoint(1:nd,Membership(ip)) + Point(1:nd,ip)
                    Kmeans%Try(it)%Size(Membership(ip)) = Kmeans%Try(it)%Size(Membership(ip)) + 1_IK
                    potential = potential + MinDistanceSq(ip)
                end do
                !%%%%%%%%%%%%%%%%%%%
            else
                call runKPP(nc, nd, np, Point, SumPoint, Membership, Kmeans%Try(it)%Size, potential)
            end if

            !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            ! Refine all centers until convergence.
            !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            Kmeans%Try(it)%nzsci = 0_IK
            Kmeans%Try(it)%niter = 1_IK

            loopCenterRefinement: do

                ! compute the new cluster centers and create the new memberships.

                !%%%%%%%%%%%%%%%%%%%
                ic = 1
                Kmeans%Try(it)%Center(1:nd,ic) = SumPoint(1:nd,ic) / real(Kmeans%Try(it)%Size(ic), kind = RK)
                do ip = 1, np
                    Kmeans%Try(it)%NormedPoint(1:nd,ip,ic) = Kmeans%Try(it)%Center(1:nd,ic) - Point(1:nd,ip)
                    Kmeans%Try(it)%MinDistanceSq(ip) = sum( Kmeans%Try(it)%NormedPoint(1:nd,ip,ic)**2 )
                    Kmeans%Try(it)%Membership(ip) = ic
                end do
                !%%%%%%%%%%%%%%%%%%%
                do ic = 2, nc - 1
                    Kmeans%Try(it)%Center(1:nd,ic) = SumPoint(1:nd,ic) / real(Kmeans%Try(it)%Size(ic), kind = RK)
                    do ip = 1, np
                        Kmeans%Try(it)%NormedPoint(1:nd,ip,ic) = Kmeans%Try(it)%Center(1:nd,ic) - Point(1:nd,ip)
                        distanceSq = sum( Kmeans%Try(it)%NormedPoint(1:nd,ip,ic)**2 )
                        if (distanceSq < Kmeans%Try(it)%MinDistanceSq(ip)) then
                            Kmeans%Try(it)%MinDistanceSq(ip) = distanceSq
                            Kmeans%Try(it)%Membership(ip) = ic
                        end if
                    end do
                end do
                !%%%%%%%%%%%%%%%%%%%
                ic = nc
                convergenceOccurred = .true.
                Kmeans%Try(it)%potential = 0._RK
                Kmeans%Try(it)%Center(1:nd,ic) = SumPoint(1:nd,ic) / real(Kmeans%Try(it)%Size(ic), kind = RK)
                do ip = 1, np
                    Kmeans%Try(it)%NormedPoint(1:nd,ip,ic) = Kmeans%Try(it)%Center(1:nd,ic) - Point(1:nd,ip)
                    distanceSq = sum( Kmeans%Try(it)%NormedPoint(1:nd,ip,ic)**2 )
                    if (distanceSq < Kmeans%Try(it)%MinDistanceSq(ip)) then
                        Kmeans%Try(it)%MinDistanceSq(ip) = distanceSq
                        Kmeans%Try(it)%Membership(ip) = ic
                    end if
                    if (Kmeans%Try(it)%Membership(ip)/=Membership(ip)) then
                        convergenceOccurred = .false.
                        Kmeans%Try(it)%Size(Membership(ip)) = Kmeans%Try(it)%Size(Membership(ip)) - 1_IK
                        Kmeans%Try(it)%Size(Kmeans%Try(it)%Membership(ip)) = Kmeans%Try(it)%Size(Kmeans%Try(it)%Membership(ip)) + 1_IK
                    end if
                    Kmeans%Try(it)%potential = Kmeans%Try(it)%potential + Kmeans%Try(it)%MinDistanceSq(ip)
                end do
                !%%%%%%%%%%%%%%%%%%%

                if (any(Kmeans%Try(it)%Size==0_IK)) then

                    Kmeans%Try(it)%niter = 1_IK
                    Kmeans%Try(it)%nzsci = Kmeans%Try(it)%nzsci + 1_IK

                    if (Kmeans%Try(it)%nzsci > Kmeans%nzsciMax) then
                        !Kmeans%Try(it)%Err%msg = "Kmeans%Try(it)%nzsci > Kmeans%nzsciMax"
                        Kmeans%Try(it)%Err%occurred = .true.
                        Kmeans%Try(it)%Err%stat = 2_IK
                        exit loopCenterRefinement
                    end if

                    ! find new random Centers and start over.

                    call runKPP(nc, nd, np, Point, SumPoint, Membership, Kmeans%Try(it)%Size, potential)

                    cycle loopCenterRefinement

#if defined DEBUG_ENABLED || defined TESTING_ENABLED
                elseif (any(Kmeans%Try(it)%Size<0_IK)) then
                    block
                        use Err_mod, only: abort
                        Kmeans%Try(it)%Err%msg = PROCEDURE_NAME//": Internal error occurred - Kmeans%Try(it)%Size < 0. Please report this issue along with circumstances to the developers."
                        call abort(Kmeans%Try(it)%Err)
                        return
                    end block
#endif

                end if

                ! if `niterMax` has reached, quit.

                if (Kmeans%Try(it)%niter < Kmeans%niterMax) then

                    ! If convergence has occurred, exit loopCenterRefinement.

                    if (convergenceOccurred .or. abs(Kmeans%Try(it)%potential-potential)/Kmeans%Try(it)%potential < Kmeans%relTol) exit loopCenterRefinement

                    Kmeans%Try(it)%niter = Kmeans%Try(it)%niter + 1_IK

                    do ip = 1, np
                        SumPoint(1:nd,Membership(ip)) = SumPoint(1:nd,Membership(ip)) - Point(1:nd,ip)
                        SumPoint(1:nd,Kmeans%Try(it)%Membership(ip)) = SumPoint(1:nd,Kmeans%Try(it)%Membership(ip)) + Point(1:nd,ip)
                        Membership(ip) = Kmeans%Try(it)%Membership(ip)
                    end do
                    potential = Kmeans%Try(it)%potential

                else

                    !Kmeans%Try(it)%Err%msg = "Kmeans%Try(it)%niter > Kmeans%niterMax"
                    Kmeans%Try(it)%Err%occurred = .true.
                    Kmeans%Try(it)%Err%stat = 1_IK
                    exit loopCenterRefinement

                end if

            end do loopCenterRefinement

            if (.not. Kmeans%Try(it)%Err%occurred) then
                if (Kmeans%Try(it)%potential < minPotential) Kmeans%itopt = it
                Kmeans%failed = .false.
            end if

        end do loopMultipleKmeansTry

    end subroutine runKmeans

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    !> \brief
    !> Find nc initial cluster centers for the Kmeans algorithm, using Kmeans++ recipe.
    !>
    !> \param[in]       nc              :   The number of clusters to be found via the Kmeans algorithm (must be > 1, though it is not verified).
    !> \param[in]       nd              :   The number of dimensions (features) of the input data array `Point`.
    !> \param[in]       np              :   The number of observations in the input data array `Point`.
    !> \param[in]       Point           :   The input array of size `(nd,np)` of the points to be clustered via the Kmeans algorithm.
    !> \param[out]      SumPoint        :   An output array of size `(nd,nc)` containing the sum of the points in each cluster.
    !>                                      When divided by the corresponding output cluster size, it will yield the final updated cluster centers.
    !> \param[out]      Membership      :   An output array of size `(np)` representing the IDs of the clusters to which individual points belong.
    !> \param[out]      potential       :   The sum of the minimum-distances-squared of all points from their corresponding group centers.
    !>
    !> \warning
    !> There is no benefit in calling this routine for `nc = 1`.
    !> As such, calling this routine with `nc = 1` will lead to undefined behavior.
    !>
    !> \todo
    !> All `minloc` queries in this procedure can be replaced with more efficient bin search algorithms.
    !>
    !> \todo
    !> The variable `CumSumDistSq` does not have to be defined as `real64`.
    !>
    !> \author
    ! Amir Shahmoradi, April 03, 2017, 2:16 PM, ICES, UTEXAS
    subroutine runKPP(nc, nd, np, Point, SumPoint, Membership, Size, potential)
#if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
        !DEC$ ATTRIBUTES DLLEXPORT :: runKPP
#endif
        use Statistics_mod, only: getRandInt

        implicit none

        integer(IK) , intent(in)    :: np, nd, nc
        real(RK)    , intent(in)    :: Point(nd,np)
        real(RK)    , intent(out)   :: SumPoint(nd,nc)
        integer(IK) , intent(out)   :: Membership(np)
        integer(IK) , intent(out)   :: Size(nc)
        real(RK)    , intent(out)   :: potential

        real(RK)                    :: dummy
        real(RK)                    :: distanceSq
        real(RK)                    :: Center(nd,nc)
        real(RK)                    :: MinDistanceSq(np)    ! An array representing the minimum distances of individual points from all clusters.
        real(RK)                    :: CumSumDistSq(0:np)   ! sum of distance squared of points up to the given index of the vector.
        integer(IK)                 :: randomInteger
        integer(IK)                 :: ic, ip, i

        CumSumDistSq(0) = 0._RK ! This element must always remain zero.

        ic = 1_IK
        Center(1:nd,ic) = Point( 1:nd , getRandInt(1,np) ) ! perhaps this random assignment is unnecessary, after all, Point is supposed to be of random order.
        do ip = 1, np
            MinDistanceSq(ip) = sum( (Center(1:nd,ic) - Point(1:nd,ip))**2 ) ! find the distance-squared of each point to the # ic cluster Center.
            CumSumDistSq(ip) = CumSumDistSq(ip-1) + MinDistanceSq(ip)
            Membership(ip) = ic
        end do

        ic = 2_IK
        call random_number(dummy)
        dummy = dummy * CumSumDistSq(np)
        Center(1:nd,ic) = Point( 1:nd, minloc(CumSumDistSq, dim = 1, mask = CumSumDistSq > dummy) - 1 ) ! -1 is essential to compensate for the 0 starting index of CumSumDistSq

        ! find the rest of cluster centers

        do ic = 2, nc - 1

            do ip = 1, np
                distanceSq = sum( (Center(1:nd,ic) - Point(1:nd,ip))**2 )
                if (distanceSq < MinDistanceSq(ip)) then
                    MinDistanceSq(ip) = distanceSq
                    Membership(ip) = ic
                end if
                CumSumDistSq(ip) = CumSumDistSq(ip-1) + MinDistanceSq(ip)
            end do

            call random_number(dummy)
            dummy = dummy * CumSumDistSq(np)
            Center(1:nd,ic+1) = Point( 1:nd , minloc(CumSumDistSq, dim = 1, mask = CumSumDistSq > dummy) - 1 )

        end do

        Size = 0_IK
        SumPoint = 0._RK
        potential = 0._RK
        do ip = 1, np
            distanceSq = sum( (Center(1:nd,nc) - Point(1:nd,ip))**2 )
            if (distanceSq < MinDistanceSq(ip)) then
                MinDistanceSq(ip) = distanceSq
                Membership(ip) = nc
            end if
            SumPoint(1:nd,Membership(ip)) = SumPoint(1:nd,Membership(ip)) + Point(1:nd,ip)
            Size(Membership(ip)) = Size(Membership(ip)) + 1_IK
            potential = potential + MinDistanceSq(ip)
        end do

    end subroutine runKPP

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end module KmeansOOP_mod ! LCOV_EXCL_LINE

max_stars_repo_name: cdslaborg/paramonte

max_issues_repo_name: cdslaborg/paramonte

max_forks_repo_name: cdslaborg/paramonte


text: #redirect wiki:Sacramento:Sky Bar


subroutine getsiga ()
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    getsiga
!   prgmmr: todling
!
! abstract:  Calculate analysis errors from Lanczos-CG results
!
! program history log:
!   2010-03-16  todling  - initial code
!   2010-05-14  todling  - update to use gsi_bundle
!   2010-05-27  todling  - gsi_4dcoupler; remove all user-specific TL-related references
!   2010-08-19  lueken   - add only to module use;no machine code, so use .f90
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

use kinds, only: i_kind, r_kind
use mpimod, only: mype
use constants, only: zero,one
use gsi_4dvar, only: ibdate,lsqrtb,jsiga
use jfunc, only: jiter
use lanczos, only : congrad_siga
use state_vectors, only: allocate_state,deallocate_state
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert_set
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert_final
use gsi_bundlemod, only: gsi_bundle
implicit none
! declare local variables
character(len=*),parameter:: myname_ = "getsiga"
type(gsi_bundle)     :: siga                      ! vector to analysis errors
integer(i_kind)      :: nymd                      ! date as in YYYYMMDD
integer(i_kind)      :: nhms                      ! time as in HHMMSS
integer(i_kind)      :: nvecs
integer(i_kind)      :: ier

! consistency checks
if (jiter/=jsiga) return
if (.not.lsqrtb) then
   write(6,*)trim(myname_),': must set lsqrt=.t. to get analysis errors'
   call stop2(331)
end if

nymd = 10000*ibdate(1)+ibdate(2)*100+ibdate(3)
nhms = 10000*ibdate(4)
if(mype==0) write(6,'(2a,i8.8,2x,i6.6)')trim(myname_),': starting to calculate analysis errors at ',&
             nymd, nhms

! allocate memory for working arrays
call allocate_state(siga)

! calculate estimate of analysis errors
call congrad_siga(siga,nvecs,ier)

! write out analysis errors
if(ier==0) then
   call gsi_4dcoupler_putpert_set (nymd,nhms,ier)
   call gsi_4dcoupler_putpert (siga,nymd,nhms,'tlm','siga')
   if(mype==0) write(6,'(2a,i5,a)')trim(myname_),': complete calculating analysis errors using ',&
                                   nvecs, ' eigenvectors'
   call gsi_4dcoupler_putpert_final(ier)
else
   if(mype==0) write(6,'(2a,i6)')trim(myname_),': failed to calculate analysis errors, ier= ', ier
endif

! clean up
call deallocate_state(siga)

return
end subroutine getsiga

subroutine view_cv (xhat,mydate,filename,writecv)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    view_cv
!   prgmmr: todling
!
! abstract:  this allow writing CV to file for visualization
!
! program history log:
!   2011-02-23  todling  - initial code
!                          (not sure we'll keep this here)
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

use kinds, only: i_kind, r_kind
use mpimod, only: mype
use constants, only: zero,one
use gsi_4dvar, only: nsubwin,lsqrtb
use state_vectors, only: allocate_state,deallocate_state
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert_set
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert_final
use gsi_bundlemod, only: gsi_bundle
use control_vectors, only: control_vector,write_cv
use state_vectors, only: allocate_state,deallocate_state,prt_state_norms
use bias_predictors, only: predictors,allocate_preds,deallocate_preds
use bias_predictors, only: write_preds
implicit none
type(control_vector)        :: xhat
integer(i_kind), intent(in) :: mydate(5) ! as in iadate or ibdate, or similar
character(len=*),intent(in) :: filename
logical,         intent(in) :: writecv   ! when .t., simply write out CV directly
! declare local variables
character(len=*),parameter:: myname_ = "view_cv"
integer(i_kind)      :: nymd                      ! date as in YYYYMMDD
integer(i_kind)      :: nhms                      ! time as in HHMMSS
integer(i_kind)      :: ii,ier
type(gsi_bundle) :: mval(nsubwin)
type(predictors) :: sbias

! in case CV not required to be transformed ...
if (writecv) then
   call write_cv(xhat,filename)
   return
else
!  if(mype==0) write(6,*) trim(myname_),': not writing CV to disk for now'
!  return
endif

! otherwise, transform CV to state-space and write out ...
nymd = 10000*mydate(1)+mydate(2)*100+mydate(3)
nhms = 10000*mydate(4)
if(mype==0) write(6,'(2a,i8.8,2x,i6.6)')trim(myname_),': start writing state on ',&
             nymd, nhms

! Allocate local variables
do ii=1,nsubwin
   call allocate_state(mval(ii))
end do
call allocate_preds(sbias)

if (lsqrtb) then
   call control2model(xhat,mval,sbias)
else
   call control2state(xhat,mval,sbias)
endif

! write out analysis errors
call gsi_4dcoupler_putpert_set (nymd,nhms,ier)
do ii=1,nsubwin
   call gsi_4dcoupler_putpert (mval(ii),nymd,nhms,'tlm',filename) ! will need to be smart for nsubwin>1
   call prt_state_norms(mval(ii),'output-state')
enddo
call gsi_4dcoupler_putpert_final (ier)
call write_preds(sbias,'preds_'//trim(filename),mype)

! Allocate local variables
call deallocate_preds(sbias)
do ii=1,nsubwin
   call deallocate_state(mval(ii))
end do

if(mype==0) write(6,'(3a)')trim(myname_),': complete writing state ', trim(filename)

return
end subroutine view_cv

subroutine view_cv_ad (xhat,mydate,filename,readcv)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    view_cv_ad
!   prgmmr: todling
!
! abstract:  this allows reading state/control vector and transforming to CV
!
! program history log:
!   2012-05-22  todling  - based on view_ad, performs opposite of view_cv
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

use kinds, only: i_kind, r_kind
use mpimod, only: mype
use constants, only: zero,one
use gsi_4dvar, only: nsubwin,lsqrtb
use state_vectors, only: allocate_state,deallocate_state
use gsi_4dcouplermod, only: gsi_4dcoupler_getpert
use gsi_bundlemod, only: gsi_bundle
use gsi_bundlemod, only: gsi_bundlegetpointer
use gsi_bundlemod, only: assignment(=)
use control_vectors, only: control_vector,read_cv,assignment(=)
use state_vectors, only: allocate_state,deallocate_state,prt_state_norms
use bias_predictors, only: predictors,allocate_preds,deallocate_preds,assignment(=)
use bias_predictors, only: read_preds
implicit none
type(control_vector)        :: xhat
integer(i_kind), intent(in) :: mydate(5) ! as in iadate or ibdate, or similar
character(len=*),intent(in) :: filename
logical,         intent(in) :: readcv    ! when .t. simply read in CV
! declare local variables
character(len=*),parameter:: myname_ = "view_cv_ad"
integer(i_kind)      :: nymd                      ! date as in YYYYMMDD
integer(i_kind)      :: nhms                      ! time as in HHMMSS
integer(i_kind)      :: ii
type(gsi_bundle) :: mval(nsubwin)
type(predictors) :: sbias

! in case CV not required to be transformed ...
if (readcv) then
   call read_cv(xhat,filename)
   return
else ! for now only
   xhat=zero
   if(mype==0) write(6,*) trim(myname_),': input vector set to zero for now'
   return
endif

! otherwise read state-vector and transform to control-space ...
nymd = 10000*mydate(1)+mydate(2)*100+mydate(3)
nhms = 10000*mydate(4)
if(mype==0) write(6,'(2a,i8.8,2x,i6.6)')trim(myname_),': start reading state on ',&
             nymd, nhms

! Allocate local variables
do ii=1,nsubwin
   call allocate_state(mval(ii))
   mval(ii)=zero
end do
call allocate_preds(sbias)
sbias=zero
xhat=zero

if (lsqrtb) then
   call control2model(xhat,mval,sbias) ! dirty trick
endif
! read in (model/state) vector
do ii=1,nsubwin
   mval(ii)=zero
   call gsi_4dcoupler_getpert (mval(ii),'tlm',filename) ! will need better for nsubwin>1
   call prt_state_norms(mval(ii),'input-state')
enddo
call read_preds(sbias,'preds_'//trim(filename))

! convert to control vector
if (lsqrtb) then
   call control2model_ad(mval,sbias,xhat)
else
   call control2state_ad(mval,sbias,xhat)
endif

! Allocate local variables
call deallocate_preds(sbias)
do ii=1,nsubwin
   call deallocate_state(mval(ii))
end do

if(mype==0) write(6,'(3a)')trim(myname_),': complete reading state ', trim(filename)
return
end subroutine view_cv_ad

subroutine view_st (sval,filename)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    view_st
!   prgmmr: todling
!
! abstract:  this allow writing CV to file for visualization
!
! program history log:
!   2011-02-23  todling  - initial code
!                          (not sure we'll keep this here)
!   2020-02-26  todling   obsbin time now in minute
!   2020-03-03  todling   split set/final from put
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:
!
!$$$ end documentation block

use kinds, only: i_kind, r_kind
use mpimod, only: mype
use constants, only: zero,one
use gsi_4dvar, only: ibdate,nobs_bins,nmn_obsbin
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert_set
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert
use gsi_4dcouplermod, only: gsi_4dcoupler_putpert_final
use gsi_bundlemod, only: gsi_bundle
implicit none
type(gsi_bundle)            :: sval(nobs_bins)
character(len=*),intent(in) :: filename
! declare local variables
character(len=*),parameter:: myname_ = "view_st"
integer(i_kind)      :: nymd                      ! date as in YYYYMMDD
integer(i_kind)      :: nhms                      ! time as in HHMMSS
integer(i_kind)      :: ii,status
integer(i_kind)      :: mydate(5)

integer(i_kind),dimension(8) :: ida,jda
real(r_kind),dimension(5)    :: fha

! initial date/time and writer
mydate = ibdate
nymd = 10000*mydate(1)+mydate(2)*100+mydate(3)
nhms = 10000*mydate(4)+mydate(5)*100

call GSI_4dCoupler_putpert_set(nymd,nhms,status)

! write out analysis errors/increments
do ii=1,nobs_bins
   nymd = 10000*mydate(1)+mydate(2)*100+mydate(3)
   nhms = 10000*mydate(4)+mydate(5)*100
   ! iwrtinc ...

   if(mype==0) write(6,'(2a,i8.8,2x,i6.6)')trim(myname_),': start writing state on ', nymd, nhms
   call gsi_4dcoupler_putpert (sval(ii),nymd,nhms,'tlm',filename)
!  call prt_state_norms(sval(ii),'output-state')

   ! increment mydate ...
   fha(:)=0.0; ida=0; jda=0
   fha(3)=nmn_obsbin! relative time interval in minutes
   ida(1)=mydate(1) ! year
   ida(2)=mydate(2) ! month
   ida(3)=mydate(3) ! day
   ida(4)=0         ! time zone
   ida(5)=mydate(4) ! hour
   ida(6)=mydate(5) ! hour
   ! Move date-time forward by nhr_assimilation hours
   call w3movdat(fha,ida,jda)
   mydate(1)=jda(1)
   mydate(2)=jda(2)
   mydate(3)=jda(3)
   mydate(4)=jda(5)
   mydate(5)=jda(6)
enddo
! finalize writer
call GSI_4dCoupler_putpert_final(status)

if(mype==0) write(6,'(3a)')trim(myname_),': complete writing state ', trim(filename)

return
end subroutine view_st