text:
program gen_be_cov2d3d_contrib
!------------------------------------------------------------------------
!  Purpose: Computes contribution of 2D Field on the balanced part
!           of 3D Field for WRFDA cv_options=6
!
!  Auothor: Syed RH Rizvi (MMM/NESL/NCAR)   Date: 02/01/2010
!           & Monika Krysta, (CTBTO, Vienna, Austria)
!
!  Note: Please acknowledge author/institute in work that uses this code.
!------------------------------------------------------------------------
   use da_control, only : stderr, stdout, filename_len
   use da_tools_serial, only : da_get_unit,da_advance_cymdh

   implicit none

   character*10        :: start_date, end_date       ! Starting and ending dates.
   character*10        :: date, new_date             ! Current date (ccyymmddhh).
   character*80        :: variable
   character*80        :: regcoeff_var
   character*10        :: variable1                  ! Variable name
   character*10        :: variable2                  ! Variable name

   real, allocatable   :: regcoeff_psi_chi(:)        ! psi/chi    regression cooefficient.
   real, allocatable   :: regcoeff_psi_t(:,:,:)      ! psi/t      regression cooefficient.
   real, allocatable   :: regcoeff_psi_ps(:,:)       ! psi/ps     regression cooefficient.
   real, allocatable   :: regcoeff_psi_rh(:,:,:)     ! psi/rh     regression cooefficient.
   real, allocatable   :: regcoeff_chi_u_t(:,:,:)    ! chi_u/t    regression coefficient
   real, allocatable   :: regcoeff_chi_u_ps(:,:)     ! chi_u/ps   regression coefficient
   real, allocatable   :: regcoeff_chi_u_rh(:,:,:)   ! chi_u/rh   regression coefficient
   real, allocatable   :: regcoeff_t_u_rh(:,:,:)     ! t_u/rh     regression coefficient
   real, allocatable   :: regcoeff_ps_u_rh(:,:)      ! ps_u/rh    regression coefficient

   character(len=filename_len)        :: filename                   ! Input filename.
   character*3         :: ce                         ! Member index -> character.
   integer             :: ni, nj, nk, nkdum                 ! Grid dimensions.
   integer             :: i, j, k, member            ! Loop counters.
   integer             :: b                          ! Bin marker.
   integer             :: sdate, cdate, edate        ! Starting, current ending dates.
   integer             :: interval                   ! Period between dates (hours).
   integer             :: ne                         ! Number of ensemble members.
   integer             :: bin_type                   ! Type of bin to average over.
   integer             :: num_bins                   ! Number of bins (3D fields).
   integer             :: num_bins2d                 ! Number of bins (3D fields).
   real                :: lat_min, lat_max           ! Used if bin_type = 2 (degrees).
   real                :: binwidth_lat               ! Used if bin_type = 2 (degrees).
   real                :: hgt_min, hgt_max           ! Used if bin_type = 2 (m).
   real                :: binwidth_hgt               ! Used if bin_type = 2 (m).
   real                :: coeffa, coeffb             ! Accumulating mean coefficients.
   logical             :: first_time                 ! True if first file.

   real, allocatable   :: field1(:,:)              ! Field 1.
   real, allocatable   :: field2(:,:,:)              ! Field 2 (full)
   real, allocatable   :: field2_balanced(:,:,:)   ! Field 2 (unbalanced)
   real*4, allocatable   :: field(:,:)   ! Field 2 (unbalanced)
   integer, allocatable:: bin(:,:,:)                 ! Bin assigned to each 3D point.
   integer, allocatable:: bin2d(:,:)                 ! Bin assigned to each 2D point.
   integer, allocatable:: bin_pts(:)                 ! Number of points in bin (3D fields).
   real, allocatable   :: covar(:)                   ! Covariance between input fields.
   real, allocatable   :: var(:)                     ! Autocovariance of field.

   real, allocatable   :: regcoeff_field1_field2(:,:) ! Regression coefficient

   namelist / gen_be_cov3d_nl / start_date, end_date, interval, &
                                ne, variable1, variable2

   integer :: ounit,iunit,namelist_unit

   stderr = 0
   stdout = 6

!---------------------------------------------------------------------------------------------
   write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
!---------------------------------------------------------------------------------------------

   call da_get_unit(ounit)
   call da_get_unit(iunit)
   call da_get_unit(namelist_unit)

   start_date = '2004030312'
   end_date = '2004033112'
   interval = 24
   ne = 1
   variable1 = 'ps_u'
   variable2 = 'rh'

   open(unit=namelist_unit, file='gen_be_cov_contrib_nl.nl', &
        form='formatted', status='old', action='read')
   read(namelist_unit, gen_be_cov3d_nl)
   close(namelist_unit)

   read(start_date(1:10), fmt='(i10)')sdate
   read(end_date(1:10), fmt='(i10)')edate
   write(6,'(4a)')' Computing covariance for fields ', variable1 , ' and ', variable2
   write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
   write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
   write(6,'(a,i8)')' Number of ensemble members at each time = ', ne

   date = start_date
   cdate = sdate

   member = 1
   write(ce,'(i3.3)')member

!  Read grid dimensions from first field file
   filename = trim(variable2)//'/'//date(1:10)
   filename = trim(filename)//'.'//trim(variable2)//'.e'//ce
   open (iunit, file = filename, form='unformatted')
   read(iunit)ni, nj, nk
   close(iunit)

!  Allocate
   allocate( bin(1:ni,1:nj,1:nk) )
   allocate( bin2d(1:ni,1:nj) )
   allocate( field1(1:ni,1:nj) )
   allocate( field2(1:ni,1:nj,1:nk) )
   allocate( field2_balanced(1:ni,1:nj,1:nk) )
   allocate( field(1:ni,1:nj) )

!  Read bin info:
   filename = 'bin.data'
   open (iunit, file = filename, form='unformatted')
   read(iunit)bin_type
   read(iunit)lat_min, lat_max, binwidth_lat
   read(iunit)hgt_min, hgt_max, binwidth_hgt
   read(iunit)num_bins, num_bins2d
   read(iunit)bin(1:ni,1:nj,1:nk)
   read(iunit)bin2d(1:ni,1:nj)
   close(iunit)

   allocate  (regcoeff_psi_chi(1:num_bins))
   allocate  (regcoeff_psi_t(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_psi_ps(1:nk,1:num_bins2d))
   allocate  (regcoeff_psi_rh(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_chi_u_t(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_chi_u_ps(1:nk,1:num_bins2d))
   allocate  (regcoeff_chi_u_rh(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_t_u_rh(1:nk,1:nk,1:num_bins2d))
   allocate  (regcoeff_ps_u_rh(1:nk,1:num_bins2d))

   allocate( regcoeff_field1_field2(1:nk,1:num_bins2d) )

   allocate( bin_pts(1:num_bins) )
   allocate( covar(1:num_bins) )
   allocate( var(1:num_bins) )
   bin_pts(:) = 0
   covar(:) = 0.0
   var(:) = 0.0

!---------------------------------------------------------------------------------------------
   write(6,'(2a)')' [2] Read fields, and calculate correlations'
!---------------------------------------------------------------------------------------------

!  Read regression coefficients
   filename = 'gen_be_stage2.dat'
   open(iunit, file = filename, form='unformatted',status='old')
   do k = 1, 2
   read (iunit)
   end do

   regcoeff_var='regcoeff_'//trim(variable1)//'_'//trim(variable2)
   do k = 1 , 9
   read (iunit) variable
   select case( trim(adjustl(variable)) )

   case ('regcoeff_psi_chi')
   read (iunit) regcoeff_psi_chi
   case ('regcoeff_psi_t')
   read (iunit) regcoeff_psi_t
   case ('regcoeff_psi_ps')
   read (iunit) regcoeff_psi_ps
   case ('regcoeff_psi_rh')
   read (iunit) regcoeff_psi_rh
   case ('regcoeff_chi_u_t')
   read (iunit) regcoeff_chi_u_t
   case ('regcoeff_chi_u_ps')
   read (iunit) regcoeff_chi_u_ps
   case ('regcoeff_chi_u_rh')
   read (iunit) regcoeff_chi_u_rh
   case ('regcoeff_t_u_rh')
   read (iunit) regcoeff_t_u_rh
   case ('regcoeff_ps_u_rh')
   read (iunit) regcoeff_ps_u_rh
   if(trim(adjustl(variable)) == trim(adjustl(regcoeff_var)) )&
   regcoeff_field1_field2 = regcoeff_ps_u_rh
   exit
   case default;
      write(6,fmt='(a)')'Read problem gen_be_stage2.dat'
      write(6,'(a,a)')' Trying to read regression coefficients:',trim(adjustl(variable))
      stop
   end select
   end do

   close(iunit)
   deallocate (regcoeff_psi_chi)
   deallocate (regcoeff_psi_t)
   deallocate (regcoeff_psi_ps)
   deallocate (regcoeff_psi_rh)
   deallocate (regcoeff_chi_u_t)
   deallocate (regcoeff_chi_u_ps)
   deallocate (regcoeff_chi_u_rh)
   deallocate (regcoeff_t_u_rh)
   deallocate (regcoeff_ps_u_rh)

   do while ( cdate <= edate )

      do member = 1, ne

         write(ce,'(i3.3)')member

!        Read first field:
         filename = trim(variable1)//'/'//date(1:10)
         filename = trim(filename)//'.'//trim(variable1)//'.e'//ce//'.01'
         open (iunit, file = filename, form='unformatted')
         read(iunit)ni, nj, nkdum
         read(iunit)field1
         close(iunit)

!        Read second field:
         filename = trim(variable2)//'/'//date(1:10)
         filename = trim(filename)//'.'//trim(variable2)//'.e'//ce
         open (iunit, file = filename, form='unformatted')
         read(iunit)ni, nj, nk
         read(iunit)field2
         close(iunit)

!        Calculate balanced contribution to field2
         do j = 1, nj
            do i = 1, ni
               b = bin2d(i,j)
               do k = 1, nk
                  field2_balanced(i,j,k) = regcoeff_field1_field2(k,b) * field1(i,j)
               end do
            end do
         end do

!        Calculate covariances:
         do k = 1, nk
            do j = 1, nj
               do i = 1, ni
                  b = bin(i,j,k)
                  bin_pts(b) = bin_pts(b) + 1
                  coeffa = 1.0 / real(bin_pts(b))
                  coeffb = real(bin_pts(b)-1) * coeffa
                  covar(b) = coeffb * covar(b) + coeffa * field2_balanced(i,j,k) * field2(i,j,k)
                  var(b) = coeffb * var(b) + coeffa * field2(i,j,k) * field2(i,j,k)
               end do
            end do
         end do

      end do  ! End loop over ensemble members.

!     Calculate next date:
      call da_advance_cymdh( date, interval, new_date )
      date = new_date
      read(date(1:10), fmt='(i10)')cdate

   end do     ! End loop over times.

   filename = trim(variable1)//'.'//trim(variable2)//'.dat'
   open (ounit, file = filename, status='unknown')
   do k = 1, nk
      do j = 1, nj
         do i = 1, ni
            b = bin(i,j,k) ! Take value from center of i dimension.
            if ( var(b) /= 0.0 ) then
               write(ounit,'(f16.8)')covar(b) / var(b)
            else
               write(ounit,'(f16.8)')0.0
            end if
         end do
      end do
   end do
   close(ounit)

   filename = trim(variable1)//'.'//trim(variable2)//'.bin'
   open (ounit, file = filename, form = 'unformatted')
   do k = 1, nk
      do j = 1, nj
         do i = 1, ni
            b = bin(i,j,k) ! Take value from center of i dimension.
            if ( var(b) /= 0.0 ) then
!               write(ounit,'(f16.8)')covar(b) / var(b)
               field(i,j) = covar(b) / var(b)
            else
!               write(ounit,'(f16.8)')0.0
               field(i,j) = 0.
            end if
         end do
      end do
      write(ounit) field
   end do
   close(ounit)

end program gen_be_cov2d3d_contrib

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
Module ply_reader

    use types
    use triangleclass

    implicit none

    public  :: read_ply
    private :: read_vert, read_faces, make_triangle

    Contains

        subroutine read_ply(filename, tarray)

            use Image

            implicit none

            character(*),                intent(IN)    :: filename
            type(triangle), allocatable, intent(INOUT) :: tarray(:)

            type(vector), allocatable :: varray(:)
            type(ivec),   allocatable :: farray(:,:)

            integer            :: u, io, verts, faces, prop, lines

            print*,'Reading: ',filename
            open(newunit=u,file=filename,iostat=io)
            if(io /= 0)stop "file not found"
            verts = 0
            faces = 0
            close(u)

            call read_header(filename, verts, faces, prop, lines)

            print*,verts,faces, prop
            allocate(varray(verts), farray(faces,1))
            call read_vert(filename, varray, lines)
            call read_faces(filename, farray, lines, verts)

            allocate(tarray(faces))
            call make_triangle(varray, farray, tarray)
            deallocate(varray,farray)

        end subroutine read_ply

        subroutine read_header(filename, verts, faces, property, i)

            implicit none

            character(len=*), intent(IN)  :: filename
            integer,          intent(OUT) :: verts, faces, property

            integer            :: io, u, i
            character(len=256) :: line

            verts = 0
            faces = 0
            property = 0

            open(newunit=u,file=filename,iostat=io)

            if(io /= 0)stop "Can't open file in ply_reader"

            i = 1
            do
                read(u, '(a)',iostat=io) line

                if(io /= 0)exit  !eof
                line = adjustl(line)

                !check if ply file
                if(i == 1 .and. verify(adjustl(line(1:3)), 'ply') /= 0)stop 'not valid ply file'
                if(i == 2 .and. verify(adjustl(line(8:12)), 'ascii') /= 0)stop "Binary .ply files not supported"

                !ignore comments
                if(verify(line(1:7), "comment") == 0)then
                    i = i + 1
                    cycle
                end if

                !get # of faces
                if(verify(line(1:13), "element face") == 0)then
                    line = adjustl(line(13:))
                    read(line,*)faces
                end if

                !get # of verts
                if(verify(line(1:14), "element vertex") == 0)then
                    line = adjustl(line(15:))
                    read(line,*)verts
                end if

                !get # of properties
                if(verify(line(1:8), "property") == 0)then
                    line = adjustl(line(9:))

                    if(verify(line(1:4), "list") == 0)then
                        continue
                    else
                        property = property + 1
                    end if
                end if

                if(verify(line(1:10), "end_header") == 0)exit
                i = i + 1

            end do
            close(u)
        end subroutine read_header

        subroutine make_triangle(varray, farray, tarray)

            implicit none

            type(vector), intent(IN)      :: varray(:)
            type(ivec),   intent(IN)      :: farray(:,:)
            type(triangle), intent(INOUT) :: tarray(:)

            type(vector) :: tmp(3), tmp2(3), tmp3(3)
            integer :: i

            do i = 1, size(farray,1)
                tmp(:) = [varray(farray(i,1)%x), varray(farray(i,1)%y), varray(farray(i,1)%z)]
                tarray(i) = triangle(rgb(255,0,0), tmp, tmp2, tmp3)
            end do

        end subroutine make_triangle

        subroutine read_vert(filename, array, lines)

            implicit none

            character(*), intent(IN)  :: filename
            type(vector), intent(OUT) :: array(:)

            integer            :: i, u, io,lines
            character(len=256) :: line
            real               :: valMax

            open(newunit=u, file=filename, iostat=io)

            read(u, '(a)',iostat=io) line
            if(io /= 0)stop  !eof

            line = adjustl(line)
            do i = 1, lines-1
                read(u, '(a)',iostat=io)line
            end do

            do i = 1, size(array)
                read(u, '(a)',iostat=io) line
                line = adjustl(line)
                read(line,*)array(i)%x, array(i)%y, array(i)%z
            end do
            close(u)

            valMax = max(maxval(abs(array%x)),maxval(abs(array%y)),maxval(abs(array%z)))
            if(valMax > 1.0)then
                array%x = array%x / valMax
                array%y = array%y / valMax
                array%z = array%z / valMax
            end if
        end subroutine read_vert

        subroutine read_faces(filename, array, lines, verts)

            implicit none

            character(*), intent(IN)   :: filename
            type(ivec),   intent(OUT)  :: array(:, :)
            integer :: lines, verts

            integer            :: i, u, io
            character(len=256) :: line

            open(newunit=u, file=filename, iostat=io)

            i = 1

            read(u, '(a)',iostat=io) line
            if(io /= 0)stop  !eof

            line = adjustl(line)
            do i = 1, lines+verts-1
                read(u, '(a)',iostat=io)line
            end do

            do i = 1, size(array)
                read(u, '(a)',iostat=io) line
                line = adjustl(line(2:))
                read(line,*)array(i,1)%x, array(i,1)%y, array(i,1)%z
                array(i,1)%x = array(i,1)%x + 1
                array(i,1)%y = array(i,1)%y + 1
                array(i,1)%z = array(i,1)%z + 1
            end do
            close(u)
        end subroutine read_faces
end module ply_reader

max_stars_repo_name: lewisfish/FortGL

max_issues_repo_name: lewisfish/FortGL

max_forks_repo_name: lewisfish/FortGL


text:
FUNCTION area_circle(r)
IMPLICIT NONE

REAL*8, PARAMETER :: Pi = 3.1415927
REAL*8, INTENT(IN) :: r
REAL*8 :: area_circle

area_circle = Pi * r * r
RETURN

END FUNCTION area_circle

max_stars_repo_name: dvuckovic/ecbuild

max_issues_repo_name: dvuckovic/ecbuild

max_forks_repo_name: dvuckovic/ecbuild


text:
c     W. H. Bell
c     A fortran subroutine and function to demonstrate the
c     principle of wrapping.
c
      SUBROUTINE COMMONS
      IMPLICIT NONE
      INCLUDE 'FORTRAN.INC'

      INTEGER I,J

      DO I=1,3
         DO J=1,2
            PRINT *,'COMMONS::INTARRAY(',I,',',J,')=',INTARRAY(I,J)
         ENDDO
      ENDDO

      DO I=1,2
         DO J=1,3
            PRINT *,'COMMONS::REALARRAY(',I,',',J,')=',REALARRAY(I,J)
         ENDDO
      ENDDO

      PRINT *,'SOMESTRING=',SOMESTRING

      RETURN
      END

      FUNCTION CALL_BACK(A,NAME)
      IMPLICIT NONE
      REAL A, C, CALL_BACK, MULT_A
      CHARACTER*(*) NAME

      PRINT *,'CALL_BACK::NAME=',NAME

      C = MULT_A(A)

      PRINT *,'CALL_BACK::C=',C

      CALL_BACK = C
      RETURN
      END

ext: for

max_stars_repo_name: williamhbell/HepCppIntro

max_issues_repo_name: williamhbell/HepCppIntro

max_forks_repo_name: williamhbell/HepCppIntro


text:
c
!> Check that the potential grid mnew is completely
!! contained in the (coarser) finest grid which stays
!! fixed, at level lbase. projec algo. Will guarantee
!! containment in all finer grids twixt them.
!!
!! If grid not contained in some coarse grid, then
!! bisect in long direction and add the two new grids
!! to the cluster (grid) list and return false.
!! EVENTUALLY this has to work, since flagged pts were
!! checked for proper nesting.
!!
!! For grid mnew (on level **lev**), check if it's properly nested in
!! - fixed base level (**lbase**) grids for current regridding
!! - level **lev**-2 grids
!!
!! If grid mnew is properly nested, return true.
!!
!! Need to check base level grids since they are fixed during
!! currign regridding.
!!
!! TODO: check these
!!
!! Need to check **lev**-2 grids since the **lev**-1 grids to be
!! generated at next step is in the range of level **lev**-2 grids.
!! If level **lev**-2 grids do not properly contain mnew, there is
!! no way the **lev**-1 grids to be generated at next step can
!! contain grid mnew.
!! Grid mnew will be projected to level **lev**-2 grids so level
!! **lev**-1 grids are guaranteed to contain grid mnew after they
!! are re-generated.
!!
!!
c
c ---------------------------------------------------------
c
      logical function nestck2(mnew,lbase,badpts,npts,numptc,icl,
     1                        nclust, nvar,naux)
c
      use amr_module
      implicit double precision (a-h,o-z)
      dimension  badpts(2,npts)
      logical baseCheck, isNested1, isNested2
      logical projecCheck

      integer   numptc(maxcl), zeroBuff
c
c ::::::::::::::::::::::: NESTCK :::::::::::::::::::::::::::::::::::

c nestck - check that the potential grid mnew is completely
c          contained in the (coarser) finest grid which stays
c          fixed, at level lbase. projec algo. will guarantee
c          containment in all finer grids twixt them.
c          if grid not contained in some coarse grid,  then
c          bisect in long direction.
c          EVENTUALLY this has to work, since flagged pts were
c          checked for proper nesting.
c
c input parameter:
c   mnew    - grid descriptor of potential grid
c   lbase   - level which stays fixed during regridding
c   badpts  - only the flagged pts. in this cluster (# icl)
c :::::::::::::::::::::::::::::::::;::::::::::::::::::::::::::::::::
c
       nestck2 = .true.
       levnew  = node(nestlevel,mnew)
       lratiox = intratx(levnew-1)
       lratioy = intraty(levnew-1)
c
!--c      # for CONVEX coarsest grid at level 1, nothing to check
       if (lbase .eq. 1)  then
            isNested1 = .true.

c POTENTIAL BUG FIX
c  need to also check that new grid can be projected to level-2 grids so is
c  properly nested.  might be an accident of grid gen that makes it stick out
c  for now check using same call for both purposes
       else

         zeroBuff = 0   ! dont count buffer zone around grid in checking
         isNested1 = baseCheck(mnew,lbase,node(ndilo,mnew),
     .                         node(ndihi,mnew),node(ndjlo,mnew),
     .                         node(ndjhi,mnew),nvar,naux,zeroBuff)
       endif

c again using new second definition of proper nesting (must have existing grid to project to
c to insure new lev-2 grid generated containing levnew grids
       levToCheck = levnew - 2
       if (levToCheck .le. 1) then
         isNested2 = .true.   ! base grid convex, no L shaped domains
       else if (levToCheck .le. lbase) then
         isNested2 = isNested1
       else
         mbuff = max(nghost,ibuff+1) ! you can use buffer zone in checking, since is only to flag points
         isNested2 = baseCheck(mnew,levToCheck,node(ndilo,mnew),
     .                         node(ndihi,mnew),node(ndjlo,mnew),
     .                         node(ndjhi,mnew),nvar,naux,mbuff)

       endif

       if (isNested1 .and. isNested2) then
           nestck2 = .true.
           go to 99
       endif
c
c  ### use grid indices coarsened by 1 level in checking
c  ### remember to offset by 1 since 1st grid cell is 0,0

c  ### grid not properly nested. bisect in long direction, and return
c  ### two clusters instead of 1.
c
 50    if (npts .gt. 1) go to 55
           write(outunit,101) levnew
           write(*,101)       levnew
 101       format(' nestck2: 1 pt. cluster at level ',i5,' still not',
     1       ' nested',/,'          pt. too close to boundary')
           write(outunit,104) badpts(1,npts),badpts(2,npts)
           write(*,104)       badpts(1,npts),badpts(2,npts)
 104       format(' non-nested flagged pt. at: ',2e15.7)
           call outtre(mstart, .false.,nvar,naux)
           call outmsh(mnew, .false.,nvar,naux)
           stop

 55    if (nclust .lt. maxcl) go to 60
           write(outunit,102) maxcl
           write(*,102)       maxcl
 102       format(' too many clusters: > ',i5,' (from nestck2) ')
           stop

 60   if (nprint) write(outunit,103) icl, npts
 103  format(' bisecting cluster ',i5,' with ',i5,' pts. in nestck2')
      if (rnode(cornxhi,mnew)-rnode(cornxlo,mnew) .gt.
     1    rnode(cornyhi,mnew) - rnode(cornylo,mnew)) then
           rmid  = (rnode(cornxhi,mnew) + rnode(cornxlo,mnew) ) / 2.
           rmid  = (node(ndihi,mnew) + node(ndilo,mnew) + 1 ) / 2.
           rmid  = rmid / lratiox
           idir = 1
       else
           rmid  = (rnode(cornyhi,mnew) + rnode(cornylo,mnew) ) / 2.
           rmid  = (node(ndjhi,mnew) + node(ndjlo,mnew) + 1) / 2.
           rmid  = rmid / lratioy
           idir = 2
       endif
c
       ipt  = 1
       ntop = npts

 90    if (badpts(idir,ipt) .lt. rmid) go to 100
c
c  ### swap with a point in top half not yet tested. keep smaller
c  ### half of rect. in bottom half
c
       temp           = badpts(1,ipt)
       badpts(1,ipt)  = badpts(1,ntop)
       badpts(1,ntop) = temp
       temp           = badpts(2,ipt)
       badpts(2,ipt)  = badpts(2,ntop)
       badpts(2,ntop) = temp
       ntop           = ntop - 1
       if (ipt .le. ntop) go to 90
       go to 110
 100   ipt = ipt +1
       if (ipt .le. ntop) go to 90
c
c  ### ntop points to top of 1st cluster (= no. of points in 1st cluster)
c
 110   numptc(icl)     = npts - ntop
       do 120 i        = icl, nclust
          nmove           = nclust + icl - i
 120      numptc(nmove+1) = numptc(nmove)
       numptc(icl)     = ntop
       nclust          = nclust + 1
       nestck2         = .false.
c
 99    return
       end

ext: f

max_stars_repo_name: navravi/amrclaw


text:
PROGRAM twoexp

! Code converted using TO_F90 by Alan Miller
! Date: 1999-03-15  Time: 11:15:27

!           SAMPLE DRIVER PROGRAM FOR NONLINEAR LEAST SQUARES ROUTINE

!            SOLUTION TO SAMPLE PROBLEM (TO 3 SIGNIFICANT DIGITS):
!               LINEAR PARAMETERS  .375, 1.94, -1.46
!               NONLINEAR PARAMETERS  .0129, .0221

USE separable_leastsq
IMPLICIT NONE
REAL (dp) :: y(50), t(50,1), alf(12), beta(8), a(50,13), w(50), eta
INTEGER   :: i, ierr, im1, iprint, iv, j, l, lpnl, n, nl, nmax, p
INTEGER, PARAMETER :: INPUT = 5, output = 6

INTERFACE
  SUBROUTINE ada (lp1, nl, n, nmax, lpp2, iv, a, inc, t, alf, isel)
    USE separable_leastsq
    IMPLICIT NONE
    INTEGER, INTENT(IN)     :: lp1
    INTEGER, INTENT(IN)     :: nl
    INTEGER, INTENT(IN)     :: n
    INTEGER, INTENT(IN)     :: nmax
    INTEGER, INTENT(IN)     :: lpp2
    INTEGER, INTENT(IN)     :: iv
    REAL (dp), INTENT(OUT)  :: a(:,:)
    INTEGER, INTENT(OUT)    :: inc(:,:)
    REAL (dp), INTENT(IN)   :: t(:,:)
    REAL (dp), INTENT(IN)   :: alf(:)
    INTEGER, INTENT(IN)     :: isel
  END SUBROUTINE ada
END INTERFACE

!              READ AND LIST THE DATA

OPEN (INPUT, FILE='a.dat', STATUS='OLD')
OPEN (output, FILE='b.out', STATUS='NEW')
nmax = 50
iprint = 1
READ (INPUT,*) n, l, nl, p, iv
WRITE (output,110) n, l, nl, p, iv
DO  i=1,n
  w(i) = 1.0_dp
  READ (INPUT,*) t(i,1:iv), y(i)
END DO
DO  i=1,n
  WRITE (output,210) i, t(i,1:iv), y(i)
END DO

READ (INPUT,120) alf(1:nl)
WRITE (output,130) alf(1:nl)
WRITE (output,140)

CALL varpro (l, nl, n, nmax, l+p+2, iv, t, y, w, ada, a, iprint,  &
             alf, beta, ierr)
WRITE (6,30) ierr
30 FORMAT (t16, 'IERR =', i7)

!              PRINT LOWER TRIANGLE OF COVARIANCE MATRIX

IF (ierr <= -4) GO TO 90
lpnl = l + nl
WRITE (output,160)
DO  i=1,lpnl
  WRITE (output,170) a(i,1:i)
  a(i,i) = SQRT(a(i,i))
END DO

!              PRINT STANDARD DEVIATIONS OF PARAMETER ESTIMATES

WRITE (output,180)
WRITE (output,170) (a(j,j),j=1,lpnl)

!              CALCULATE AND PRINT LOWER TRIANG. OF CORRELATION MATRIX

IF (lpnl == 1) GO TO 60
DO  i=2,lpnl
  im1 = i-1
  DO  j=1,im1
    a(i,j) = a(i,j) / (a(i,i)*a(j,j))
  END DO
END DO

60 WRITE (output,190)
DO  i=1,lpnl
  a(i,i) = 1.0_dp
  WRITE (output,170) a(i,1:i)
END DO

!              PRINT RESIDUALS

WRITE (output,200)
DO  i=1,n
  eta = y(i) - a(i,lpnl+1) / SQRT(w(i))
  WRITE (output,210) i, w(i), t(i,1:iv), y(i), eta, a(i,lpnl+1)
END DO

90 STOP

110 FORMAT ('1   NON-LINEAR LEAST SQUARES PROBLEM'//  &
    ' NUMBER OF OBSERVATIONS =',i5, '   NUMBER OF LINEAR PARAMETERS =',i4//  &
    ' NUMBER OF NONLINEAR PARAMETERS =',i4,  &
    ' NUMBER OF NONVANISHING PARTIAL DERIVATIVES =',i4//  &
    ' NUMBER OF INDEPENDENT VARIABLES =',i4// '    I    T(I)            Y(I)'//)
120 FORMAT (4E20.7)
130 FORMAT ('0 INITIAL NONLINEAR PARAMETERS'//(4E20.7))
140 FORMAT ('0',50('*'))
160 FORMAT (/'       COVARIANCE MATRIX'/)
170 FORMAT (8E15.7)
180 FORMAT ('0    STANDARD DEVIATIONS OF PARAMETER ESTIMATES'/)
190 FORMAT (/'       CORRELATION MATRIX'/)
200 FORMAT (/'   I         W(I)            T(I)            Y(I)',  &
    '     PREDICTED Y     WEIGHTED RESIDUAL'//)
210 FORMAT (i5,7E16.7)
END PROGRAM twoexp

SUBROUTINE ada (lp1, nl, n, nmax, lpp2, iv, a, inc, t, alf, isel)

!        EVALUATE THE FUNCTIONS PHI AND THEIR PARTIAL DERIVATIVES
!        D PHI(J)/D ALF(K), AT THE SAMPLE POINTS T(I).
!        ISEL = 1  MEANS COMPUTE BOTH FUNCTIONS AND DERIVATIVES AND
!                  INITIALIZE INC AND CONSTANT PHI'S
!             = 2  MEANS COMPUTE ONLY THE NONCONSTANT FUNCTIONS PHI
!             = 3  MEANS COMPUTE ONLY THE DERIVATIVES.

!        THIS PARTICULAR ROUTINE IS FOR FITTING TWO EXPONENTIALS AND
!        ONE CONSTANT TERM:

!        ETA(C, ALF; T) = C  + C  * EXP(-ALF *T) + C * EXP(-ALF * T)
!                          1    2           1       3          2

!        WHERE C IS THE VECTOR OF LINEAR PARAMETERS TO BE DETERMINED,
!        AND ALF IS THE VECTOR OF NONLINEAR PARAMETERS TO BE DETERMINED.
!        HERE N = 33, L = 3, NL = P = 2, NCFUN = 1, PHI(1) = 1, PHI(2) =
!        EXP(-ALF(1)*T), AND PHI(3) = EXP(-ALF(2)*T).  THIS EXAMPLE AND
!        TEST DATA ARE TAKEN FROM REFERENCE 3 (SEE VARPRO).

!     ..................................................................

USE separable_leastsq
IMPLICIT NONE

INTEGER, INTENT(IN)     :: lp1
INTEGER, INTENT(IN)     :: nl
INTEGER, INTENT(IN)     :: n
INTEGER, INTENT(IN)     :: nmax
INTEGER, INTENT(IN)     :: lpp2
INTEGER, INTENT(IN)     :: iv
REAL (dp), INTENT(OUT)  :: a(:,:)        ! a(nmax,lpp2)
INTEGER, INTENT(OUT)    :: inc(:,:)      ! inc(12,lp1)
REAL (dp), INTENT(IN)   :: t(:,:)        ! t(nmax)
REAL (dp), INTENT(IN)   :: alf(:)        ! alf(nl)
INTEGER, INTENT(IN)     :: isel

! Local variable

INTEGER :: i

SELECT CASE ( isel )
  CASE (1)
    GO TO 10
  CASE (2)
    GO TO 30
  CASE (3)
    GO TO 50
END SELECT

!              ISEL = 1, IN THIS CASE THE INCIDENCE MATRIX INC IS:

!               0  1  0  0     <- ALF(1)
!               0  0  1  0     <- ALF(2)

10 inc = 0
inc(1,2) = 1
inc(2,3) = 1

!              ISEL = 1, SET CONSTANT FUNCTION PHI(1) (REMOVE IF NO
!              CONSTANT FUNCTIONS)

a(1:n,1) = 1.0_dp

!              ISEL = 1 OR 2, COMPUTE NONCONSTANT FUNCTIONS PHI

30 DO  i=1,n
  a(i,2) = EXP(-alf(1)*t(i,1))
  a(i,3) = EXP(-alf(2)*t(i,1))
END DO
!                           COLUMN L+1 = 4 IS LEFT FOR WORKSPACE.
IF (isel == 2) GO TO 70

!              ISEL = 1 OR 3, COMPUTE DERIVATIVES DPHI(I) / D ALF(J)

50 DO  i=1,n
  a(i,5) = -t(i,1)*EXP(-alf(1)*t(i,1))
  a(i,6) = -t(i,1)*EXP(-alf(2)*t(i,1))
END DO

70 RETURN
END SUBROUTINE ada

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


! ******************************************************************************************************************************** !
! atchem_data_netCDF.f90
! ATCHEM netCDF ROUTINES
! ******************************************************************************************************************************** !

MODULE atchem_data_netCDF

  USE gem_netcdf
  USE atchem_lib
  IMPLICIT NONE
  SAVE

CONTAINS

  ! ****************************************************************************************************************************** !
  ! SAVE NETCDF RESTART DATA
  ! ****************************************************************************************************************************** !
  SUBROUTINE sub_data_netCDF_ncrstsave(dum_name,dum_yr,dum_iou)
    ! -------------------------------------------------------- !
    ! DUMMY ARGUMENTS
    ! -------------------------------------------------------- !
    character(LEN=*),INTENT(IN)::dum_name                      !
    REAL,INTENT(in)::dum_yr                                    !
    INTEGER,INTENT(OUT)::dum_iou                               !
    ! -------------------------------------------------------- !
    ! DEFINE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    integer::ia,l
    integer::loc_ntrec,loc_iou
    integer::loc_id_lonm,loc_id_latm,loc_id_lon_e,loc_id_lat_e
    integer,dimension(1:2)::loc_it_1
    integer,dimension(1:3)::loc_it_2
    character(127)::loc_title,loc_timunit
    character(7)::loc_string_year
    real::loc_c0,loc_c1
    real,dimension(0:n_i)::loc_lon_e
    real,dimension(0:n_j)::loc_lat_e
    real,dimension(n_i,n_j)::loc_ij,loc_ij_mask
    ! -------------------------------------------------------- !
    ! INITIALIZE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    loc_c0 = 0.0
    loc_c1 = 1.0
    ! -------------------------------------------------------- !
    ! WRITE TO FILE
    ! -------------------------------------------------------- !
    ! -------------------------------------------------------- ! open file
    call sub_opennew(dum_name,loc_iou)
    ! -------------------------------------------------------- ! start definitions
    call sub_redef(loc_iou)
    ! -------------------------------------------------------- ! set global attributes
    loc_string_year = fun_conv_num_char_n(8,int(dum_yr))
    loc_title = 'ATCHEM restart @ year '//loc_string_year
    call sub_putglobal(loc_iou,dum_name,loc_title,string_ncrunid,loc_timunit)
    ! -------------------------------------------------------- ! define dimensions
    call sub_defdim ('lon', loc_iou, n_i, loc_id_lonm)
    call sub_defdim ('lat', loc_iou, n_j, loc_id_latm)
    call sub_defdim ('lon_edges', loc_iou, n_i+1, loc_id_lon_e)
    call sub_defdim ('lat_edges', loc_iou, n_j+1, loc_id_lat_e)
    ! -------------------------------------------------------- ! define 1d data (t)
    loc_it_1(1) = loc_id_lonm
    call sub_defvar ('lon',loc_iou,1,loc_it_1,loc_c0,loc_c0,'X','D','longitude of the t grid','longitude','degrees_east')
    loc_it_1(1) = loc_id_latm
    call sub_defvar ('lat',loc_iou,1,loc_it_1,loc_c0,loc_c0,'Y','D','latitude of the t grid','latitude','degrees_north')
    loc_it_1(1) = loc_id_lon_e
    call sub_defvar ('lon_edges',loc_iou,1,loc_it_1,loc_c0,loc_c0,' ','D' ,'longitude of t grid edges',' ','degrees')
    loc_it_1(1) = loc_id_lat_e
    call sub_defvar ('lat_edges',loc_iou,1,loc_it_1,loc_c0,loc_c0,' ','D','latitude of t grid edges',' ','degrees')
    loc_it_2(1) = loc_id_lonm
    loc_it_2(2) = loc_id_latm
    ! -------------------------------------------------------- ! define (2D) tracer variables
    DO l=1,n_l_atm
       ia = conv_iselected_ia(l)
       call sub_defvar('atm_'//trim(string_atm(ia)),loc_iou,2,loc_it_2,loc_c0,loc_c0,' ','F', &
            & string_longname_atm(ia),'Atmosphere tracer - '//trim(string_atm(ia)),' ')
    end DO
    ! -------------------------------------------------------- ! end definitions
    call sub_enddef (loc_iou)
    call sub_sync(loc_iou)
    ! -------------------------------------------------------- !
    loc_ntrec = 1
    ! -------------------------------------------------------- ! write 1D variables
    call sub_putvar1d ('lon',loc_iou,n_i,loc_ntrec,n_i,phys_atm(ipa_lon,:,1),loc_c1,loc_c0)
    call edge_maker (1,loc_lon_e,phys_atm(ipa_lon,:,1),phys_atm(ipa_lone,:,1),phys_atm(ipa_dlon,:,1),n_i)
    call sub_putvar1d ('lon_edges',loc_iou,n_i+1,loc_ntrec,n_i+1,loc_lon_e,loc_c1,loc_c0)
    call sub_putvar1d ('lat',loc_iou,n_j,loc_ntrec,n_j,phys_atm(ipa_lat,1,:),loc_c1,loc_c0)
    call edge_maker (1,loc_lat_e, phys_atm(ipa_lat,1,:),phys_atm(ipa_latn,1,:),phys_atm(ipa_dlat,1,:), n_j)
    call sub_putvar1d ('lat_edges',loc_iou,n_j+1,loc_ntrec,n_j+1,loc_lat_e,loc_c1,loc_c0)
    ! -------------------------------------------------------- ! write (2D) tracer variables
    loc_ij_mask(:,:) = 1.0
    loc_ij(:,:)= 0.0
    DO l=1,n_l_atm
       ia = conv_iselected_ia(l)
       loc_ij(:,:) = atm(ia,:,:)
       call sub_putvar2d('atm_'//trim(string_atm(ia)),loc_iou,n_i,n_j,loc_ntrec,loc_ij(:,:),loc_ij_mask(:,:))
    end DO
    ! -------------------------------------------------------- ! close file and return IOU
    call sub_closefile(loc_iou)
    dum_iou = loc_iou
    ! -------------------------------------------------------- !
    ! END
    ! -------------------------------------------------------- !
  END SUBROUTINE sub_data_netCDF_ncrstsave
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! INITIALIZE netCDF
  ! ****************************************************************************************************************************** !
  SUBROUTINE sub_data_netCDF_ncrstinit_OLD(dum_name,dum_iou)
    ! ======================================================== !
    !     initialize netcdf restart file
    !
    !     input:
    !       dum_name  = name of variable to be defined
    !       dum_iou   = unit
    !
    ! loosely based on code by: M. Eby
    ! ======================================================== !
    ! -------------------------------------------------------- !
    ! DUMMY ARGUMENTS
    ! -------------------------------------------------------- !
    character(LEN=*),INTENT(IN)::dum_name
    INTEGER,INTENT(OUT)::dum_iou
    ! -------------------------------------------------------- !
    ! DEFINE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    character(255) :: loc_title, loc_timunit
    real           :: loc_c0, loc_c1
    integer        :: loc_it(6), loc_id_time, loc_id_lonm
    integer        :: loc_id_latm, loc_id_lon_e, loc_id_xu, loc_id_yu
    integer        :: loc_id_lat_e, loc_id_xu_e, loc_id_yu_e
    integer        :: loc_id_misc
    ! -------------------------------------------------------- !
    ! INITIALIZE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    loc_c0 = 0.0
    loc_c1 = 1.0
    ! -------------------------------------------------------- !
    ! WRITE TO FILE
    ! -------------------------------------------------------- !
    ! -------------------------------------------------------- ! open file
    call sub_opennew(dum_name,dum_iou)
    ! -------------------------------------------------------- ! start definitions
    call sub_redef(dum_iou)
    ! -------------------------------------------------------- ! set global attributes
    loc_title = 'Time averaged integrals'
    write (loc_timunit,'(a)') 'Year mid-point'
    call sub_putglobal(dum_iou,dum_name,loc_title,string_ncrunid,loc_timunit)
    ! -------------------------------------------------------- ! define dimensions
    call sub_defdim ('time', dum_iou, const_integer_zero, loc_id_time)
    call sub_defdim ('lon', dum_iou, n_i, loc_id_lonm)
    call sub_defdim ('lat', dum_iou, n_j, loc_id_latm)
    call sub_defdim ('xu', dum_iou, n_i, loc_id_xu)
    call sub_defdim ('yu', dum_iou, n_j, loc_id_yu)
    call sub_defdim ('lon_edges', dum_iou, n_i+1, loc_id_lon_e)
    call sub_defdim ('lat_edges', dum_iou, n_j+1, loc_id_lat_e)
    call sub_defdim ('xu_edges', dum_iou, n_i+1, loc_id_xu_e)
    call sub_defdim ('yu_edges', dum_iou, n_j+1, loc_id_yu_e)
    call sub_defdim ('para', dum_iou, const_integer_one, loc_id_misc)
    ! -------------------------------------------------------- ! define 1d data (t)
    loc_it(1) = loc_id_time
    call sub_defvar ('time', dum_iou, 1, loc_it, loc_c0, loc_c0, 'T', 'D' &
         &, 'Year', 'time', trim(loc_timunit))
    call sub_defvar ('year', dum_iou, 1, loc_it, loc_c0, loc_c0, ' ', 'F','year', ' ',' ')
    ! -------------------------------------------------------- ! define 1d data (x, y)
    loc_it(1) = loc_id_lonm
    call sub_defvar ('lon', dum_iou, 1, loc_it, loc_c0, loc_c0, 'X', 'D' , &
         &'longitude of the t grid', 'longitude', 'degrees_east')
    loc_it(1) = loc_id_latm
    call sub_defvar ('lat', dum_iou, 1, loc_it, loc_c0, loc_c0, 'Y', 'D' , &
         &'latitude of the t grid', 'latitude', 'degrees_north')
    loc_it(1) = loc_id_xu
    call sub_defvar ('xu', dum_iou, 1, loc_it, loc_c0, loc_c0, 'X', 'D' , &
         &'longitude of the u grid', 'longitude', 'degrees_east')
    loc_it(1) = loc_id_yu
    call sub_defvar ('yu', dum_iou, 1, loc_it, loc_c0, loc_c0, 'Y', 'D' , &
         &'latitude of the u grid', 'latitude', 'degrees_north')
    loc_it(1) = loc_id_lon_e
    call sub_defvar ('lon_edges', dum_iou, 1, loc_it, loc_c0, loc_c0, ' ', 'D' , &
         &'longitude of t grid edges', ' ', 'degrees')
    loc_it(1) = loc_id_lat_e
    call sub_defvar ('lat_edges', dum_iou, 1, loc_it, loc_c0, loc_c0, ' ', 'D' , &
         &'latitude of t grid edges', ' ', 'degrees')
    loc_it(1) = loc_id_xu_e
    call sub_defvar ('xu_edges', dum_iou, 1, loc_it, loc_c0, loc_c0, ' ', 'D' , &
         &'longitude of u grid edges', ' ', 'degrees')
    loc_it(1) = loc_id_yu_e
    call sub_defvar ('yu_edges', dum_iou, 1, loc_it, loc_c0, loc_c0, ' ', 'D' , &
         &'latitude of u grid edges', ' ', 'degrees')
    ! -------------------------------------------------------- ! close file
    call sub_enddef (dum_iou)
    call sub_closefile (dum_iou)
    ! -------------------------------------------------------- !
    ! END
    ! -------------------------------------------------------- !
  END SUBROUTINE sub_data_netCDF_ncrstinit_OLD
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CREATE ATCHEM RESTART FILE
  ! ****************************************************************************************************************************** !
  SUBROUTINE sub_data_netCDF_ncrstsave_OLD(dum_yr)
    ! -------------------------------------------------------- !
    ! DUMMY ARGUMENTS
    ! -------------------------------------------------------- !
    REAL,INTENT(in)::dum_yr
    ! -------------------------------------------------------- !
    ! DEFINE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    character(255) :: loc_name
    real           :: loc_c0, loc_c1
    integer        :: loc_i, loc_iou, loc_ntrec
    real,dimension(n_i+1) :: loc_lon_e, loc_xu_e
    real,dimension(n_j+1) :: loc_lat_e, loc_yu_e
    logical :: loc_defined
    ! -------------------------------------------------------- !
    ! INITIALIZE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    loc_c0 = 0.0
    loc_c1 = 1.0
    !
    loc_lon_e(:)       = 0.0
    loc_lat_e(:)       = 0.0
    loc_xu_e(:)        = 0.0
    loc_yu_e(:)        = 0.0
    ! -------------------------------------------------------- ! set local netCDF variables
    loc_name = string_ncrst
    loc_iou = ncrst_iou
    loc_ntrec = ncrst_ntrec
    ! -------------------------------------------------------- ! open file and get latest record number
    loc_defined = .true.
    loc_i = 0
    if (loc_ntrec .eq. 0) then
       loc_defined = .false.
       loc_i = 1
    end if
    call sub_opennext (loc_name, dum_yr, loc_i, loc_ntrec, loc_iou)
    ! -------------------------------------------------------- ! write 1d data: time
    call sub_putvars  ('time', loc_iou, loc_ntrec, dum_yr, loc_c1, loc_c0)
    call sub_putvarIs  ('year', loc_iou, loc_ntrec, nint(dum_yr), loc_c1, loc_c0)
    if(.not. loc_defined) then
       ! ----------------------------------------------------- ! write 1d data: x
       call sub_putvar1d ('lon', loc_iou, n_i, loc_ntrec, n_i,phys_atm(ipa_lon,:,1), loc_c1, loc_c0)
       call edge_maker (1, loc_lon_e, phys_atm(ipa_lon,:,1),phys_atm(ipa_lone,:,1), phys_atm(ipa_dlon,:,1), n_i)
       call sub_putvar1d ('lon_edges', loc_iou, n_i+1, loc_ntrec, n_i+1,loc_lon_e, loc_c1, loc_c0)
       call sub_putvar1d ('xu', loc_iou, n_i, loc_ntrec, n_i,loc_lon_e(1:n_i), loc_c1, loc_c0)
       call edge_maker (2, loc_xu_e, phys_atm(ipa_lon,:,1),phys_atm(ipa_lone,:,1), phys_atm(ipa_dlon,:,1), n_i)
       call sub_putvar1d ('xu_edges', loc_iou, n_i+1, loc_ntrec, n_i+1,loc_xu_e, loc_c1, loc_c0)
       ! ----------------------------------------------------- ! write 1d data: y
       call sub_putvar1d ('lat', loc_iou, n_j, loc_ntrec, n_j,phys_atm(ipa_lat,1,:), loc_c1, loc_c0)
       call edge_maker (1, loc_lat_e, phys_atm(ipa_lat,1,:),phys_atm(ipa_latn,1,:), phys_atm(ipa_dlat,1,:), n_j)
       call sub_putvar1d ('lat_edges', loc_iou, n_j+1, loc_ntrec, n_j+1,loc_lat_e, loc_c1, loc_c0)
       call sub_putvar1d ('yu', loc_iou, n_j, loc_ntrec, n_j,loc_lat_e(1:n_j), loc_c1, loc_c0)
       call edge_maker (2, loc_yu_e, phys_atm(ipa_lat,1,:),phys_atm(ipa_latn,1,:), phys_atm(ipa_dlat,1,:), n_j)
       call sub_putvar1d ('yu_edges', loc_iou, n_j+1, loc_ntrec, n_j+1,loc_yu_e, loc_c1, loc_c0)
    end if
    ! -------------------------------------------------------- !
    ncrst_ntrec = loc_ntrec
    ncrst_iou = loc_iou
    call sub_sync(loc_iou)
    ! -------------------------------------------------------- !
    ! END
    ! -------------------------------------------------------- !
  END SUBROUTINE sub_data_netCDF_ncrstsave_OLD
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! WRITE NETCDF RESTART DATA
  ! ****************************************************************************************************************************** !
  SUBROUTINE sub_data_netCDF_ncrstwrite_OLD()
    ! -------------------------------------------------------- !
    ! DEFINE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    INTEGER::l,ia
    INTEGER::loc_iou,loc_ntrec
    CHARACTER(len=255)::loc_unitsname
    real,DIMENSION(1:n_i,1:n_j)::loc_ij,loc_mask
    ! -------------------------------------------------------- !
    ! INITIALIZE LOCAL VARIABLES
    ! -------------------------------------------------------- !
    loc_iou = ncrst_iou
    loc_ntrec = ncrst_ntrec
    loc_mask(:,:) = 1.0
    loc_unitsname = 'atm'
    ! -------------------------------------------------------- !
    ! WRITE DATA
    ! -------------------------------------------------------- !
    DO l=1,n_l_atm
       ia = conv_iselected_ia(l)
       loc_ij(:,:) = atm(ia,:,:)
       call sub_adddef_netcdf(loc_iou,3,'atm_'//trim(string_atm_tname(l)), &
            & trim(string_atm_tlname(l)),trim(string_atm_unit(l)),atm_mima(l,1),atm_mima(l,2))
       call sub_putvar2d('atm_'//trim(string_atm(ia)),loc_iou,n_i,n_j, &
            & loc_ntrec,loc_ij(:,:),loc_mask(:,:))
    END DO
    ! -------------------------------------------------------- !
    ! END
    ! -------------------------------------------------------- !
  END SUBROUTINE sub_data_netCDF_ncrstwrite_OLD
  ! ****************************************************************************************************************************** !

END MODULE atchem_data_netCDF


text:
!
! basic search routines for range search
! (with and without distance list)
!

subroutine search_d2(node,h,ppos,npnear,nblist,d2list)
  include 'globals.h'
  integer, intent(in):: node
  real, intent(in) :: h,ppos(3)
  integer, intent(inout)::npnear,nblist(*)
  real, intent(inout) :: d2list(*)
  integer :: i,daughter,nstack,cnode,cstack(nsubcell*81)
  real :: dist2,xnode,ynode,znode,dx,dy,dz,hsearch2

  if(periodic.AND.h.GT.rsize/2) call terror('search h>rsize/2')
  hsearch2=h**2
  cstack(1)=node
  nstack=1
  do while(nstack.GT.0.AND.nstack.LE.nsubcell*80)
    daughter=cstack(nstack)
    nstack=nstack-1
    if(daughter.LT.nbods1) then
      if(daughter.GT.0) then
        dx=ppos(1)-pos(daughter,1)
        dy=ppos(2)-pos(daughter,2)
        dz=ppos(3)-pos(daughter,3)
        if(dx.GE.hboxsize.AND.periodic) dx=dx-pboxsize
        if(dx.LT.-hboxsize.AND.periodic) dx=dx+pboxsize
        if(dy.GE.hboxsize.AND.periodic) dy=dy-pboxsize
        if(dy.LT.-hboxsize.AND.periodic) dy=dy+pboxsize
        if(dz.GE.hboxsize.AND.periodic) dz=dz-pboxsize
        if(dz.LT.-hboxsize.AND.periodic) dz=dz+pboxsize
        dist2=dx**2+dy**2+dz**2
        if(dist2.LT.hsearch2) then
          npnear=npnear+1
          nblist(npnear)=daughter
          d2list(npnear)=dist2
        endif
      endif
    else
      xnode=bottom(daughter,1)+0.5*cellsize(daughter)
      ynode=bottom(daughter,2)+0.5*cellsize(daughter)
      znode=bottom(daughter,3)+0.5*cellsize(daughter)
      dx=ppos(1)-xnode
      dy=ppos(2)-ynode
      dz=ppos(3)-znode
      if(dx.GE.hboxsize.AND.periodic) xnode=xnode+pboxsize
      if(dx.LT.-hboxsize.AND.periodic) xnode=xnode-pboxsize
      if(dy.GE.hboxsize.AND.periodic) ynode=ynode+pboxsize
      if(dy.LT.-hboxsize.AND.periodic) ynode=ynode-pboxsize
      if(dz.GE.hboxsize.AND.periodic) znode=znode+pboxsize
      if(dz.LT.-hboxsize.AND.periodic) znode=znode-pboxsize
      dist2=0
      if(ppos(1).LT.xnode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(1)-xnode+0.5*cellsize(daughter))**2
      else
        if(ppos(1).GT.xnode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(1)-xnode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(2).LT.ynode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(2)-ynode+0.5*cellsize(daughter))**2
      else
        if(ppos(2).GT.ynode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(2)-ynode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(3).LT.znode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(3)-znode+0.5*cellsize(daughter))**2
      else
        if(ppos(3).GT.znode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(3)-znode-0.5*cellsize(daughter))**2
        endif
      endif
      if(dist2.LT.hsearch2) then
        cstack(nstack+1)=subp(daughter,8)
        cstack(nstack+2)=subp(daughter,7)
        cstack(nstack+3)=subp(daughter,6)
        cstack(nstack+4)=subp(daughter,5)
        cstack(nstack+5)=subp(daughter,4)
        cstack(nstack+6)=subp(daughter,3)
        cstack(nstack+7)=subp(daughter,2)
        cstack(nstack+8)=subp(daughter,1)
        nstack=nstack+8
      endif
    endif
  enddo
  if(nstack.NE.0) then
!$omp critical
    call terror(" search_d2 error")
!$omp end critical
  endif
end subroutine search_d2

subroutine search(node,h,ppos,npnear,nblist)
  include 'globals.h'
  integer, intent(in):: node
  real, intent(in) :: h,ppos(3)
  integer, intent(inout)::npnear,nblist(*)
  integer :: i,daughter,nstack,cnode,cstack(nsubcell*81)
  real :: dist2,xnode,ynode,znode,dx,dy,dz,hsearch2
  if(periodic.AND.h.GT.rsize/2) call terror('search h>rsize/2')
  hsearch2=h**2
  cstack(1)=node
  nstack=1
  do while(nstack.GT.0.AND.nstack.LE.nsubcell*80)
    daughter=cstack(nstack)
    nstack=nstack-1
    if(daughter.LT.nbods1) then
      if(daughter.GT.0) then
        dx=ppos(1)-pos(daughter,1)
        dy=ppos(2)-pos(daughter,2)
        dz=ppos(3)-pos(daughter,3)
        if(dx.GE.hboxsize.AND.periodic) dx=dx-pboxsize
        if(dx.LT.-hboxsize.AND.periodic) dx=dx+pboxsize
        if(dy.GE.hboxsize.AND.periodic) dy=dy-pboxsize
        if(dy.LT.-hboxsize.AND.periodic) dy=dy+pboxsize
        if(dz.GE.hboxsize.AND.periodic) dz=dz-pboxsize
        if(dz.LT.-hboxsize.AND.periodic) dz=dz+pboxsize
        dist2=dx**2+dy**2+dz**2
        if(dist2.LT.hsearch2) then
          npnear=npnear+1
          nblist(npnear)=daughter
        endif
      endif
    else
      xnode=bottom(daughter,1)+0.5*cellsize(daughter)
      ynode=bottom(daughter,2)+0.5*cellsize(daughter)
      znode=bottom(daughter,3)+0.5*cellsize(daughter)
      dx=ppos(1)-xnode
      dy=ppos(2)-ynode
      dz=ppos(3)-znode
      if(dx.GE.hboxsize.AND.periodic) xnode=xnode+pboxsize
      if(dx.LT.-hboxsize.AND.periodic) xnode=xnode-pboxsize
      if(dy.GE.hboxsize.AND.periodic) ynode=ynode+pboxsize
      if(dy.LT.-hboxsize.AND.periodic) ynode=ynode-pboxsize
      if(dz.GE.hboxsize.AND.periodic) znode=znode+pboxsize
      if(dz.LT.-hboxsize.AND.periodic) znode=znode-pboxsize
      dist2=0
      if(ppos(1).LT.xnode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(1)-xnode+0.5*cellsize(daughter))**2
      else
        if(ppos(1).GT.xnode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(1)-xnode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(2).LT.ynode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(2)-ynode+0.5*cellsize(daughter))**2
      else
        if(ppos(2).GT.ynode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(2)-ynode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(3).LT.znode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(3)-znode+0.5*cellsize(daughter))**2
      else
        if(ppos(3).GT.znode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(3)-znode-0.5*cellsize(daughter))**2
        endif
      endif
      if(dist2.LT.hsearch2) then
        cstack(nstack+1)=subp(daughter,8)
        cstack(nstack+2)=subp(daughter,7)
        cstack(nstack+3)=subp(daughter,6)
        cstack(nstack+4)=subp(daughter,5)
        cstack(nstack+5)=subp(daughter,4)
        cstack(nstack+6)=subp(daughter,3)
        cstack(nstack+7)=subp(daughter,2)
        cstack(nstack+8)=subp(daughter,1)
        nstack=nstack+8
      endif
    endif
  enddo
  if(nstack.NE.0) then
!$omp critical
    call terror(" search error")
!$omp end critical
  endif
end subroutine search

!
! basic search routines for complementary range search
! (with and without distance list)
!

subroutine comsearch_d2(node,h,delta,ppos,npnear,nblist,d2list)
  include 'globals.h'
  integer, intent(in):: node
  real, intent(in) :: h,ppos(3),delta
  integer, intent(inout)::npnear,nblist(*)
  real,intent(inout) ::  d2list(*)
  integer :: i,daughter,nstack,cnode,cstack(nsubcell*81)
  real :: dist2,xnode,ynode,znode,dx,dy,dz,hsearch2

  if(periodic.AND. &
      (h.GT.rsize/2.OR.(2*hsmooth(node)+delta).GT.rsize/2)) &
    call terror('comsearch h>rsize/2')
  hsearch2=h**2
  cstack(1)=node
  nstack=1
  do while(nstack.GT.0.AND.nstack.LE.nsubcell*80)
    daughter=cstack(nstack)
    nstack=nstack-1
    if(daughter.LT.nbods1) then
      if(daughter.GT.0) then
        dx=ppos(1)-pos(daughter,1)
        dy=ppos(2)-pos(daughter,2)
        dz=ppos(3)-pos(daughter,3)
        if(dx.GE.hboxsize.AND.periodic) dx=dx-pboxsize
        if(dx.LT.-hboxsize.AND.periodic) dx=dx+pboxsize
        if(dy.GE.hboxsize.AND.periodic) dy=dy-pboxsize
        if(dy.LT.-hboxsize.AND.periodic) dy=dy+pboxsize
        if(dz.GE.hboxsize.AND.periodic) dz=dz-pboxsize
        if(dz.LT.-hboxsize.AND.periodic) dz=dz+pboxsize
        dist2=dx**2+dy**2+dz**2
        if(dist2.LT.hsearch2.OR. &
                   dist2.LT.(2*hsmooth(daughter)+delta)**2) then
          npnear=npnear+1
          nblist(npnear)=daughter
          d2list(npnear)=dist2
        endif
      endif
    else
      xnode=bottom(daughter,1)+0.5*cellsize(daughter)
      ynode=bottom(daughter,2)+0.5*cellsize(daughter)
      znode=bottom(daughter,3)+0.5*cellsize(daughter)
      dx=ppos(1)-xnode
      dy=ppos(2)-ynode
      dz=ppos(3)-znode
      if(dx.GE.hboxsize.AND.periodic) xnode=xnode+pboxsize
      if(dx.LT.-hboxsize.AND.periodic) xnode=xnode-pboxsize
      if(dy.GE.hboxsize.AND.periodic) ynode=ynode+pboxsize
      if(dy.LT.-hboxsize.AND.periodic) ynode=ynode-pboxsize
      if(dz.GE.hboxsize.AND.periodic) znode=znode+pboxsize
      if(dz.LT.-hboxsize.AND.periodic) znode=znode-pboxsize
      dist2=0
      if(ppos(1).LT.xnode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(1)-xnode+0.5*cellsize(daughter))**2
      else
        if(ppos(1).GT.xnode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(1)-xnode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(2).LT.ynode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(2)-ynode+0.5*cellsize(daughter))**2
      else
        if(ppos(2).GT.ynode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(2)-ynode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(3).LT.znode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(3)-znode+0.5*cellsize(daughter))**2
      else
        if(ppos(3).GT.znode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(3)-znode-0.5*cellsize(daughter))**2
        endif
      endif
      if(dist2.LT.hsearch2.OR. &
         dist2.LT.(2*hsmooth(daughter)+delta)**2) then
        cstack(nstack+1)=subp(daughter,8)
        cstack(nstack+2)=subp(daughter,7)
        cstack(nstack+3)=subp(daughter,6)
        cstack(nstack+4)=subp(daughter,5)
        cstack(nstack+5)=subp(daughter,4)
        cstack(nstack+6)=subp(daughter,3)
        cstack(nstack+7)=subp(daughter,2)
        cstack(nstack+8)=subp(daughter,1)
        nstack=nstack+8
      endif
    endif
  enddo
  if(nstack.NE.0) then
!$omp critical
    call terror(" comsearch_d2 error")
!$omp end critical
  endif
end subroutine comsearch_d2

subroutine comsearch(node,h,delta,ppos,npnear,nblist)
  include 'globals.h'
  integer, intent(in):: node
  real, intent(in) :: h,ppos(3),delta
  integer, intent(inout)::npnear,nblist(*)
  integer :: i,daughter,nstack,cnode,cstack(nsubcell*81)
  real :: dist2,xnode,ynode,znode,dx,dy,dz,hsearch2

  if(periodic.AND. &
      (h.GT.rsize/2.OR.(2*hsmooth(node)+delta).GT.rsize/2)) then
    print*,h,rsize,hsmooth(node),delta
    print*,node,ppos
    call terror('comsearch h>rsize/2')
  endif
  hsearch2=h**2
  cstack(1)=node
  nstack=1
  do while(nstack.GT.0.AND.nstack.LE.nsubcell*80)
    daughter=cstack(nstack)
    nstack=nstack-1
    if(daughter.LT.nbods1) then
      if(daughter.GT.0) then
        dx=ppos(1)-pos(daughter,1)
        dy=ppos(2)-pos(daughter,2)
        dz=ppos(3)-pos(daughter,3)
        if(dx.GE.hboxsize.AND.periodic) dx=dx-pboxsize
        if(dx.LT.-hboxsize.AND.periodic) dx=dx+pboxsize
        if(dy.GE.hboxsize.AND.periodic) dy=dy-pboxsize
        if(dy.LT.-hboxsize.AND.periodic) dy=dy+pboxsize
        if(dz.GE.hboxsize.AND.periodic) dz=dz-pboxsize
        if(dz.LT.-hboxsize.AND.periodic) dz=dz+pboxsize
        dist2=dx**2+dy**2+dz**2
        if(dist2.LT.hsearch2.OR. &
             dist2.LT.(2*hsmooth(daughter)+delta)**2) then
          npnear=npnear+1
          nblist(npnear)=daughter
        endif
      endif
    else

      xnode=bottom(daughter,1)+0.5*cellsize(daughter)
      ynode=bottom(daughter,2)+0.5*cellsize(daughter)
      znode=bottom(daughter,3)+0.5*cellsize(daughter)
      dx=ppos(1)-xnode
      dy=ppos(2)-ynode
      dz=ppos(3)-znode
      if(dx.GE.hboxsize.AND.periodic) xnode=xnode+pboxsize
      if(dx.LT.-hboxsize.AND.periodic) xnode=xnode-pboxsize
      if(dy.GE.hboxsize.AND.periodic) ynode=ynode+pboxsize
      if(dy.LT.-hboxsize.AND.periodic) ynode=ynode-pboxsize
      if(dz.GE.hboxsize.AND.periodic) znode=znode+pboxsize
      if(dz.LT.-hboxsize.AND.periodic) znode=znode-pboxsize
      dist2=0
      if(ppos(1).LT.xnode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(1)-xnode+0.5*cellsize(daughter))**2
      else
        if(ppos(1).GT.xnode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(1)-xnode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(2).LT.ynode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(2)-ynode+0.5*cellsize(daughter))**2
      else
        if(ppos(2).GT.ynode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(2)-ynode-0.5*cellsize(daughter))**2
        endif
      endif
      if(ppos(3).LT.znode-0.5*cellsize(daughter)) then
        dist2=dist2+(ppos(3)-znode+0.5*cellsize(daughter))**2
      else
        if(ppos(3).GT.znode+0.5*cellsize(daughter)) then
          dist2=dist2+(ppos(3)-znode-0.5*cellsize(daughter))**2
        endif
      endif
      if(dist2.LT.hsearch2.OR. &
           dist2.LT.(2*hsmooth(daughter)+delta)**2) then
        cstack(nstack+1)=subp(daughter,8)
        cstack(nstack+2)=subp(daughter,7)
        cstack(nstack+3)=subp(daughter,6)
        cstack(nstack+4)=subp(daughter,5)
        cstack(nstack+5)=subp(daughter,4)
        cstack(nstack+6)=subp(daughter,3)
        cstack(nstack+7)=subp(daughter,2)
        cstack(nstack+8)=subp(daughter,1)
        nstack=nstack+8
      endif
    endif
  enddo
  if(nstack.NE.0) then
!$omp critical
    call terror(" comsearch error")
!$omp end critical
  endif
end subroutine comsearch

!
! particle-index search routines for range search
! (with and without distance list)
!
subroutine psearch(node,p,npnear,nblist)
  include 'globals.h'
  integer, intent(in):: node,p
  real :: h,ppos(3)
  integer, intent(inout)::npnear,nblist(*)
  h=2*hsmooth(p)
  ppos=pos(p,1:3)
  call search(node,h,ppos,npnear,nblist)
end subroutine

subroutine psearch_d2(node,p,npnear,nblist,d2list)
  include 'globals.h'
  integer, intent(in):: node,p
  real :: h,ppos(3)
  integer, intent(inout)::npnear,nblist(*)
  real  :: d2list(*)
  h=2*hsmooth(p)
  ppos=pos(p,1:3)
  call search_d2(node,h,ppos,npnear,nblist,d2list)
end subroutine

!
! particle-index search routines for complementary range search
! (with and without distance list)
!

subroutine pcomsearch(node,p,npnear,nblist)
  include 'globals.h'
  integer, intent(in):: node,p
  real :: h,ppos(3), delta
  integer, intent(inout)::npnear,nblist(*)
  h=2*hsmooth(p)
  ppos=pos(p,1:3)
  call comsearch(node,h,0.,ppos,npnear,nblist)
end subroutine

subroutine pcomsearch_d2(node,p,npnear,nblist,d2list)
  include 'globals.h'
  integer, intent(in):: node,p
  real :: h,ppos(3),delta
  integer, intent(inout)::npnear,nblist(*)
  real :: d2list(*)
  h=2*hsmooth(p)
  ppos=pos(p,1:3)
  delta=0.
  call comsearch_d2(node,h,delta,ppos,npnear,nblist,d2list)
end subroutine

!
! particle-index search routines for conditional (ie reusable) range search
! (with and without distance list)
!

subroutine pcond_srch(node,p,nneigh,nblist)
  include 'globals.h'
  real h,h2,ppos(3)
  integer nneigh,p,node,nblist(*)
  h=hsmooth(p)
  ppos=pos(p,1:3)
  call cond_srch(ppos,h,node,nneigh,nblist)
end subroutine

subroutine pcond_srch_d2(node,p,nneigh,nblist,d2list)
  include 'globals.h'
  real h,h2,ppos(3)
  integer nneigh,p,node,nblist(*)
  real :: d2list(*)
  h=hsmooth(p)
  ppos=pos(p,1:3)
  call cond_srch_d2(ppos,h,node,nneigh,nblist,d2list)
end subroutine

!
! driver search routines for conditional (ie reusable) range search
! (with and without distance list)
!

subroutine cond_srch(ppos,h,node,nneigh,nblist)
  include 'globals.h'
  real h,ppos(3),dist
  integer nneigh,node,nblist(*)
  ncalls=ncalls+1
  if(searchreuse.EQ.0.OR.reuseflag.EQ.0) then
    searchpos=ppos
    searchh=h
    searchn=0
  else
    dist=sqrt(sum((ppos-searchpos)**2))
    if(h+dist/2.LE.searchh.AND.h.GT.searchmagic*searchh) then
      nneigh=searchn
      return
    endif
    searchn=0
    searchpos=ppos
    searchh=h
  endif
  nsearches=nsearches+1
  call search(node,2*searchh,searchpos,searchn,nblist)
  nneigh=searchn
  searchreuse=1
end subroutine

subroutine cond_srch_d2(ppos,h,node,nneigh,nblist,d2list)
  include 'globals.h'
  real h,ppos(3),dist
  integer nneigh,node,nblist(*)
  real :: d2list(*)
  ncalls=ncalls+1
  if(searchreuse.EQ.0.OR.reuseflag.EQ.0) then
    searchpos=ppos
    searchh=h
    searchn=0
  else
    dist=sqrt(sum((ppos-searchpos)**2))
    if(h+dist/2.LE.searchh.AND.h.GT.searchmagic*searchh) then
      nneigh=searchn
      return
    endif
    searchn=0
    searchpos=ppos
    searchh=h
  endif
  nsearches=nsearches+1
  call search_d2(node,2*searchh,searchpos,searchn,nblist,d2list)
  nneigh=searchn
  searchreuse=1
end subroutine

!
! particle-index search routines for conditional (ie reusable)
! and complementary range search
! (with and without distance list)
!

subroutine pcond_comsrch(node,p,nneigh,nblist)
  include 'globals.h'
  real h,h2,ppos(3),delta
  integer nneigh,p,node,nblist(*)
  h=hsmooth(p)
  ppos=pos(p,1:3)
  delta=0.
  call cond_comsrch(ppos,h,delta,node,nneigh,nblist)
end subroutine

subroutine pcond_comsrch_d2(node,p,nneigh,nblist,d2list)
  include 'globals.h'
  real h,h2,ppos(3),delta
  integer nneigh,p,node,nblist(*)
  real :: d2list(*)
  h=hsmooth(p)
  ppos=pos(p,1:3)
  delta=0.
  call cond_comsrch_d2(ppos,h,delta,node,nneigh,nblist,d2list)
end subroutine

!
! driver search routines for conditional (ie reusable)
!  and complementary range search
! (with and without distance list)
!
subroutine cond_comsrch(ppos,h,delta,node,nneigh,nblist)
  include 'globals.h'
  real h,ppos(3),dist,delta
  integer nneigh,node,nblist(*)
  ncalls=ncalls+1
  if(searchreuse.EQ.0.OR.reuseflag.EQ.0) then
    searchpos=ppos
    searchh=h
    searchn=0
    searchdelta=delta
  else
    dist=sqrt(sum((ppos-searchpos)**2))
    if(h+dist/2.LE.searchh.AND. &
        delta+dist.LE.searchdelta.AND. &
         h.GT.searchmagic*searchh) then
      nneigh=searchn
      return
    endif
    searchn=0
    searchpos=ppos
    searchh=h
    searchdelta=delta
  endif
  nsearches=nsearches+1
  call comsearch(node,2*searchh,searchdelta,searchpos,searchn,nblist)
  nneigh=searchn
  searchreuse=1
end subroutine

subroutine cond_comsrch_d2(ppos,h,delta,node,nneigh,nblist,d2list)
  include 'globals.h'
  real h,ppos(3),dist,delta
  integer nneigh,node,nblist(*)
  real d2list(*)
  ncalls=ncalls+1
  if(searchreuse.EQ.0.OR.reuseflag.EQ.0) then
    searchpos=ppos
    searchh=h
    searchn=0
    searchdelta=delta
  else
    dist=sqrt(sum((ppos-searchpos)**2))
    if(h+dist/2.LE.searchh.AND. &
         delta+dist.LE.searchdelta.AND. &
         h.GT.searchmagic*searchh) then
      nneigh=searchn
      return
    endif
    searchn=0
    searchpos=ppos
    searchh=h
    searchdelta=delta
  endif
  nsearches=nsearches+1
  call comsearch_d2(node,2*searchh,searchdelta, &
                       searchpos,searchn,nblist,d2list)
  nneigh=searchn
  searchreuse=1
end subroutine

!
! pre search routines
!  these try to determine a search such that multiple particles
!  can reuse the same search list

subroutine presearch(k,kmax,nbuf,buf,ntodo,todo)
  include 'globals.h'
  integer :: k,kmax,nbuf,buf(nbuf),todo(nbuf),ntodo,ib,p,q,nneigh
  real :: dist2,spos(3),ppos(3),ph,qh,sh,qpos(3)
  ntodo=0
  if(buf(1).EQ.0) then
    p=pactive(k)
    ntodo=ntodo+1
    todo(ntodo)=p
    if(hsmooth(p).GT.0) then
      ppos=pos(p,1:3); ph=hsmooth(p)
      spos=ppos; sh=ph
      do ib=1,MIN(nbuf-1,kmax-k)
        if(buf(ib+1).NE.0) cycle
        q=pactive(k+ib)
        qpos=pos(q,1:3); qh=hsmooth(q)
        dist2=sum((ppos-qpos)**2)
        if(dist2.LT.ph**2) then
          buf(ib+1)=1
          ntodo=ntodo+1; todo(ntodo)=q
          spos=((ntodo-1)*spos)/ntodo+qpos/ntodo
        endif
      enddo
      do ib=1,ntodo
        q=todo(ib)
        qpos=pos(q,1:3); qh=hsmooth(q)
        dist2=sum((spos-qpos)**2)
        sh=MAX(sh,qh+sqrt(dist2)/2)
      enddo
      sh=sh*1.02
      searchreuse=0
      call cond_srch(spos,sh,root,nneigh,srlist)
    endif
  endif
  do ib=1,nbuf-1
    buf(ib)=buf(ib+1)
  enddo
  buf(nbuf)=0
end subroutine

subroutine precomsearch(k,kmax,nbuf,buf,ntodo,todo)
  include 'globals.h'
  integer :: k,kmax,nbuf,buf(nbuf),todo(nbuf),ntodo,ib,p,q,nneigh
  real :: dist2,spos(3),ppos(3),ph,qh,sh,qpos(3),delta
  ntodo=0
  if(buf(1).EQ.0) then
    p=pactive(k)
    ntodo=ntodo+1
    todo(ntodo)=p
    if(hsmooth(p).GT.0) then
      ppos=pos(p,1:3); ph=hsmooth(p)
      spos=ppos; sh=ph; delta=0
      do ib=1,MIN(nbuf-1,kmax-k)
        if(buf(ib+1).NE.0) cycle
        q=pactive(k+ib)
        qpos=pos(q,1:3); qh=hsmooth(q)
        dist2=sum((ppos-qpos)**2)
        if(dist2.LT.ph**2) then
          buf(ib+1)=1
          ntodo=ntodo+1; todo(ntodo)=q
          spos=((ntodo-1)*spos)/ntodo+qpos/ntodo
        endif
      enddo
      do ib=1,ntodo
        q=todo(ib)
        qpos=pos(q,1:3); qh=hsmooth(q)
        dist2=sum((spos-qpos)**2)
        sh=MAX(sh,qh+sqrt(dist2)/2)
        delta=MAX(delta,sqrt(dist2))
      enddo
      searchreuse=0
      call cond_comsrch(spos,sh,delta,root,nneigh,srlist)
    endif
  endif
  do ib=1,nbuf-1
    buf(ib)=buf(ib+1)
  enddo
  buf(nbuf)=0
end subroutine

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
      subroutine real_default(r1,r2,r3)
          implicit none
          real, intent(in) :: r1
          real, intent(inout) :: r2
          real, intent(out) :: r3
          r3 = r1 + r2
      end subroutine
      subroutine real_x_len(r10,r11,r12,r13,r14,r15)
          implicit none
          real*4, intent(in) :: r10
          real*4, intent(inout) :: r11
          real*4, intent(out) :: r12
          real*8, intent(in) :: r13
          real*8, intent(inout) :: r14
          real*8, intent(out) :: r15
          r12 = r10 + r11
          r15 = r13 + r14
      end subroutine
      subroutine real_kind_x(r7,r8,r9,r10,r11,r12)
          implicit none
          real(kind=4), intent(in) :: r7
          real(kind=4), intent(inout) :: r8
          real(kind=4), intent(out) :: r9
          real(kind=8), intent(in) :: r10
          real(kind=8), intent(inout) :: r11
          real(kind=8), intent(out) :: r12
          r9 = r7 + r8
          r12 = r11 + r10
      end subroutine
      ! subroutine real_kind_call(r1,r2,r3,r4,r5,r6)
          ! implicit none
          ! real(kind=kind(0.0)), intent(in) :: r1
          ! real(kind=kind(0.0)), intent(inout) :: r2
          ! real(kind=kind(0.0)), intent(out) :: r3
          ! real(kind=kind(0.0D0)), intent(in) :: r4
          ! real(kind=kind(0.0D0)), intent(inout) :: r5
          ! real(kind=kind(0.0D0)), intent(out) :: r6
          ! r3 = r2 + r1
          ! r6 = r4 + r5
      ! end subroutine
      ! subroutine real_srk_call(r1,r2,r3,r4,r5,r6,r7,r8,r9)
          ! implicit none
          ! real(kind=selected_real_kind(1)), intent(in) :: r1
          ! real(kind=selected_real_kind(1)), intent(inout) :: r2
          ! real(kind=selected_real_kind(1)), intent(out) :: r3
          ! real(kind=selected_real_kind(7)), intent(in) :: r4
          ! real(kind=selected_real_kind(7)), intent(inout) :: r5
          ! real(kind=selected_real_kind(7)), intent(out) :: r6
          ! real(kind=selected_real_kind(14)), intent(in) :: r7
          ! real(kind=selected_real_kind(14)), intent(inout) :: r8
          ! real(kind=selected_real_kind(14)), intent(out) :: r9
          ! r3 = r1 + r2
          ! r6 = r4 + r5
          ! r9 = r7 + r8
      ! end subroutine

max_stars_repo_name: wilsonify/fwrap

max_issues_repo_name: fwrap/fwrap

max_forks_repo_name: fwrap/fwrap


text:
module write_plotfileLM_module

  use ml_layout_module
  use multifab_module
  use fabio_module
  use convert_stag_module
  use convert_rhoc_to_c_module
  use probin_multispecies_module, only: nspecies, plot_stag
  use probin_common_module, only: prob_lo, prob_hi

  implicit none

contains

  subroutine write_plotfileLM(mla,name,rho,rhotot,Temp,umac,pres,istep,dx,time)

    type(ml_layout),    intent(in)    :: mla
    character(len=*),   intent(in)    :: name
    type(multifab),     intent(inout) :: rho(:)
    type(multifab),     intent(in)    :: rhotot(:)
    type(multifab),     intent(in)    :: Temp(:)
    type(multifab),     intent(in)    :: umac(:,:)
    type(multifab),     intent(in)    :: pres(:)
    integer,            intent(in)    :: istep
    real(kind=dp_t),    intent(in)    :: dx(:,:),time

    ! local variables
    character(len=20), allocatable  :: plot_names(:)
    character(len=20)               :: plot_names_stagx(1)
    character(len=20)               :: plot_names_stagy(1)
    character(len=20)               :: plot_names_stagz(1)
    character(len=20)               :: plotfile_name
    character(len=20)               :: plotfile_namex
    character(len=20)               :: plotfile_namey
    character(len=20)               :: plotfile_namez
    integer                         :: i,dm,n,nlevs

    ! multifab of size nlevs
    type(multifab), allocatable     :: plotdata(:)
    type(multifab), allocatable     :: plotdata_stag(:,:)

    type(multifab) :: conc(mla%nlevel)

    type(bl_prof_timer), save :: bpt

    call build(bpt,"write_plotfileLM")

    nlevs = mla%nlevel
    dm = mla%dim

    ! rho + species (rho) + nspeces (conc) + Temp + dm (averaged umac) + dm (shifted umac) + pres
    allocate(plot_names(2*nspecies+2*dm+3))
    allocate(plotdata(nlevs))
    allocate(plotdata_stag(nlevs,dm))

    plot_names(1) = "rho"
    do n=1,nspecies
       write(plot_names(n+1),'(a,i0)') "rho", n
    enddo
    do n=1,nspecies
       write(plot_names(nspecies+n+1),'(a,i0)') "c", n
    enddo
    plot_names(2*nspecies+2) = "Temp"
    plot_names(2*nspecies+3) = "averaged_velx"
    plot_names(2*nspecies+4) = "averaged_vely"
    if (dm > 2) plot_names(2*nspecies+5) = "averaged_velz"
    plot_names(2*nspecies+dm+3) = "shifted_velx"
    plot_names(2*nspecies+dm+4) = "shifted_vely"
    if (dm > 2) plot_names(2*nspecies+dm+5) = "shifted_velz"
    plot_names(2*nspecies+2*dm+3) = "pres"

    plot_names_stagx(1) = "velx"
    plot_names_stagy(1) = "vely"
    plot_names_stagz(1) = "velz"

    ! compute concentrations
    do n=1,nlevs
       call multifab_build(conc(n),mla%la(n),nspecies,0)
    end do
    call convert_rhoc_to_c(mla,rho,rhotot,conc,.true.)

    ! build plotdata for 2*nspecies+2*dm+3 and 0 ghost cells
    do n=1,nlevs
       call multifab_build(plotdata(n),mla%la(n),2*nspecies+2*dm+3,0)
       do i=1,dm
          call multifab_build_edge(plotdata_stag(n,i), mla%la(n), 1, 0, i)
       end do
    enddo

    ! copy rhotot, rho, conc, and Temp into plotdata
    do n = 1,nlevs
       call multifab_copy_c(plotdata(n),1           ,rhotot(n),1,       1,0)
       call multifab_copy_c(plotdata(n),2           ,rho(n)   ,1,nspecies,0)
       call multifab_copy_c(plotdata(n),nspecies+2  ,conc(n)  ,1,nspecies,0)
       call multifab_copy_c(plotdata(n),2*nspecies+2,Temp(n)  ,1,1       ,0)
    enddo

    ! vel averaged
    do i=1,dm
       call average_face_to_cc(mla,umac(:,i),1,plotdata,2*nspecies+2+i,1)
    end do

    ! vel shifted
    do i=1,dm
       call shift_face_to_cc(mla,umac(:,i),1,plotdata,2*nspecies+dm+2+i,1)
    end do

    ! pressure
    do n = 1,nlevs
       call multifab_copy_c(plotdata(n),2*nspecies+2*dm+3,pres(n),1,1,0)
    enddo

    ! copy staggered velocity and momentum into plotdata_stag
    do n=1,nlevs
       do i=1,dm
          call multifab_copy_c(plotdata_stag(n,i),1,umac(n,i),1,1)
       end do
    end do

    ! define the name of the plotfile that will be written
    write(unit=plotfile_name,fmt='(a,i8.8)') name, istep
    if ( parallel_IOProcessor() ) then
      write(*,'(2A)') "Saving PLOT FILEs to directory ", trim(plotfile_name)
      write(*,*)
    end if
    write(unit=plotfile_namex,fmt='("stagx",i8.8)') istep
    write(unit=plotfile_namey,fmt='("stagy",i8.8)') istep
    if (dm > 2) then
       write(unit=plotfile_namez,fmt='("stagz",i8.8)') istep
    end if

    ! write the plotfile
    call fabio_ml_multifab_write_d(plotdata, mla%mba%rr(:,1), trim(plotfile_name), &
                                   plot_names, mla%mba%pd(1), prob_lo, prob_hi, &
                                   time, dx(1,:))

    if (plot_stag) then
       ! staggered plotfiles
       ! the data appears "shifted" in amrvis, but you can post process these and still keep
       ! all the boundary data
       call fabio_ml_multifab_write_d(plotdata_stag(:,1), mla%mba%rr(:,1), trim(plotfile_namex), &
                                      plot_names_stagx, mla%mba%pd(1), prob_lo, prob_hi, &
                                      time, dx(1,:))

       call fabio_ml_multifab_write_d(plotdata_stag(:,2), mla%mba%rr(:,1), trim(plotfile_namey), &
                                      plot_names_stagy, mla%mba%pd(1), prob_lo, prob_hi, &
                                      time, dx(1,:))

       if (dm > 2) then
          call fabio_ml_multifab_write_d(plotdata_stag(:,3), mla%mba%rr(:,1), trim(plotfile_namez), &
                                         plot_names_stagz, mla%mba%pd(1), prob_lo, prob_hi, &
                                         time, dx(1,:))
       end if
    end if

    ! make sure to destroy the multifab or you'll leak memory
    do n=1,nlevs
       call multifab_destroy(conc(n))
       call multifab_destroy(plotdata(n))
      do i=1,dm
         call multifab_destroy(plotdata_stag(n,i))
      end do
    enddo

    deallocate(plotdata)
    deallocate(plotdata_stag)
    deallocate(plot_names)

    call destroy(bpt)

  end subroutine write_plotfileLM

end module write_plotfileLM_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
 SUBROUTINE texplit( )
 !********************************************************************
 !
 ! *** time stepping routine
 !
 !********************************************************************
 USE ctrl_db, ONLY: tdtime, ttime, ndoft, nheat, npoin, therm, neqt
 USE curv_db
 !      USE kinc_db

 USE ift_db, ONLY : npret,lctmp,nprev,prtmp,fxtem
 USE npo_db, ONLY : tempe,dtemp,tresi,heata,heati,iftmp,trate,tmass
 USE heat_db, ONLY : ncsur
 IMPLICIT NONE
! local variables
 INTEGER (kind=4), PARAMETER :: nn = 1000000, nn1 = 1000001
 INTEGER (kind=4) :: ieq,i,n1,nt0,j
 REAL (kind=8) :: dt2,facts(nheat),factv(npret)
 REAL (kind=8), SAVE :: tdtold
! Functions
 REAL (kind=8):: functs

 !   compute prescribed temperatures
 DO i=1,npret
   factv(i) = functs (lctmp(i),ttime)*prtmp(nprev+1,i)
 END DO

 IF(npret > 0)THEN
   DO i=1,nprev
     prtmp(i,npret+1)= DOT_PRODUCT(prtmp(i,1:npret),factv(1:npret))
   END DO
 END IF

 dtemp = tempe                      !keep original values

 IF( therm )THEN                    !for heat transfer problems

   IF(ttime == 0d0) tdtold = 0d0    !first time
   dt2 = (tdtime+tdtold)/2d0        !average time increment

   IF( ncsur > 0 ) CALL tcsurf( )    !compute convection and radiation terms

   trate = -tresi                   !initializes residual

   !   compute and add external equivalent heats
   DO i=1,nheat
     facts(i) = functs (heata(i),ttime)*heati(neqt+1,i)
   END DO
   IF(nheat > 0)THEN
     n1  = nheat+1
     DO ieq=1,neqt
       heati(ieq,n1) = DOT_PRODUCT(heati(ieq,1:nheat),facts(1:nheat))
       trate(ieq) = trate(ieq) + heati(ieq,n1)
     END DO
   END IF
   !   integrate equations
   DO ieq = 1,neqt
     trate(ieq) =  trate(ieq)/tmass(ieq)*dt2      !increment in fact
   END DO

   !  update temperatures

   DO i=1,npoin             !loop over each node
     DO j=1,ndoft
       nt0 = iftmp(j,i)       !DOF
       SELECT CASE (nt0)
       CASE (1: )             ! active DOF
         tempe(j,i) = tempe(j,i) + trate(nt0)
       CASE (0 )
         !nothing
       CASE (-nn:-1)          ! fixed DOF
         tempe(j,i) = fxtem(-nt0)
       CASE (:-nn1)           ! prescribed DOF
         tempe(j,i) = prtmp(-nt0-nn,npret+1)
       END SELECT
     END DO
   END DO

   tdtold = tdtime

 ELSE  ! itemp only

   !  update temperatures

   DO i=1,npoin             !loop over each node
     DO j=1,ndoft             !loop over each node
       nt0 = iftmp(j,i)       !DOF
       SELECT CASE (nt0)
       CASE (0 )
         !nothing
       CASE (-nn:-1)
         tempe(j,i) = fxtem(-nt0)
       CASE (:-nn1)
         tempe(j,i) = prtmp(-nt0-nn,npret+1)
       END SELECT
     END DO
   END DO

 END IF

 dtemp = tempe - dtemp  !incremental temperature

 RETURN
 END SUBROUTINE texplit

max_stars_repo_name: jerebenitez/IFE-simpact-openfoam

max_issues_repo_name: jerebenitez/IFE-simpact-openfoam

max_forks_repo_name: jerebenitez/IFE-simpact-openfoam


text:
    SUBROUTINE SGEES_F95( A, WR, WI, VS, SELECT, SDIM, INFO )
!
!  -- LAPACK95 interface driver routine (version 3.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September, 2000
!
!  .. USE STATEMENTS ..
   USE LA_PRECISION, ONLY: WP => SP
   USE LA_AUXMOD, ONLY: ERINFO, LSAME
   USE F77_LAPACK, ONLY: GEES_F77 => LA_GEES
!  USE LA_EXTERNAL, ONLY: SELECT
!  .. IMPLICIT STATEMENT ..
   IMPLICIT NONE
   INTERFACE
      LOGICAL FUNCTION SELECT(WR, WI)
         USE LA_PRECISION, ONLY: WP => SP
         REAL(WP), INTENT(IN) :: WR, WI
      END FUNCTION SELECT
   END INTERFACE
!  .. SCALAR ARGUMENTS ..
   INTEGER, INTENT(OUT), OPTIONAL :: INFO, SDIM
!  .. ARRAY ARGUMENTS ..
   REAL(WP), INTENT(INOUT) :: A(:,:)
   REAL(WP), INTENT(OUT) :: WR(:), WI(:)
   REAL(WP), INTENT(OUT), OPTIONAL, TARGET :: VS(:,:)
!  .. EXTERNAL ARGUMENTS ..
   OPTIONAL :: SELECT
!----------------------------------------------------------------------
!
! Purpose
! =======
!
!      LA_GEES computes for a real/complex square matrix A, the
! eigenvalues, the real-Schur/complex-Schur form T , and, optionally, the
! matrix of Schur vectors Z, where Z is orthogonal/unitary. This gives the
! Schur factorization
!                        A = Z*T*Z^H.
! Optionally, it also orders the eigenvalues on the diagonal of the Schur
! form so that selected eigenvalues are at the top left. The leading
! columns of Z then form an orthonormal basis for the invariant subspace
! corresponding to the selected eigenvalues.
!      A real matrix is in real-Schur form if it is block upper triangular
! with 1 by 1 and 2 by 2 blocks along the main diagonal. 2 by 2 blocks are
! standardized in the form
!                       [ a  b ]
! 		        [ c  a ]
! where b*c < 0. The eigenvalues of such a block are a +/- Sqrt(b*c).
! A complex matrix is in complex-Schur form if it is upper triangular.
!
! =========
!
!        SUBROUTINE LA_GEES( A, <w>, VS=vs, SELECT=select, &
!                                   SDIM=sdim, INFO=info )
!           <type>(<wp>), INTENT(INOUT) :: A(:,:)
!           <type>(<wp>), INTENT(OUT) :: <w(:)>
!           <type>(<wp>), INTENT(OUT), OPTIONAL :: VS(:,:)
!           INTERFACE
!               LOGICAL FUNCTION SELECT(<w(j)>)
!                  <type>(<wp>), INTENT(IN) :: <w(j)>
!               END FUNCTION SELECT
!           END INTERFACE
!           OPTIONAL :: SELECT
!           INTEGER, INTENT(OUT), OPTIONAL :: SDIM, INFO
!        where
!           <type> ::= REAL | COMPLEX
!           <wp>   ::= KIND(1.0) | KIND(1.0D0)
!           <w>    ::= WR, WI | W
!           <w(:)> ::= WR(:), WI(:) | W(:)
!           <w(j)> ::= WR(j) , WI(j) |  W(j)
!
! Arguments
! =========
!
! A        (input/output) REAL or COMPLEX square array, shape (:,:).
!          On entry, the matrix A.
!          On exit, the Schur form T.
! <w>      (output) REAL or COMPLEX array, shape (:) with size(w) =
!          size(A,1).
!          The computed eigenvalues in the order in which they appear on
!          the diagonal of the Schur form T.
!          <w(:)> ::= WR(:), WI(:) | W(:),
!          where
!          WR(:), WI(:) are of REAL type (for the real and imaginary
! 	   parts) and W(:) is of COMPLEX type.
!          Note: If A is real, then a complex-conjugate pair appear
!          consecutively, with the eigenvalue having the positive
! 	   imaginary part appearing first.
! VS       Optional (output) REAL or COMPLEX square array, shape (:,:)
!          with size(VS,1) = size(A,1).
!          The matrix Z of Schur vectors.
! SELECT   Optional (input) LOGICAL FUNCTION.
!          LOGICAL FUNCTION SELECT( <w(j)> )
!            <type>(<wp>), INTENT(IN) :: <w(j)>
!          where
!            <type> ::= REAL | COMPLEX
!            <wp>   ::= KIND(1.0) | KIND(1.0D0)
!            <w(j)> ::= WR(j) , WI(j) | W(j)
!          1. SELECT must be declared as EXTERNAL or as an explicit
!             interface in the calling (sub)program.
!          2. SELECT is called by LA_GEES for every computed eigenvalue
! 	    w(j) (but only once for a complex conjugate pair when A is
! 	    real). It is used to select the eigenvalues that will be
!             ordered to the top left of the Schur form. The eigenvalue
! 	    w(j) is selected if SELECT(w(j)) has the value .TRUE.
!          3. A selected complex eigenvalue may no longer satisfy
! 	    SELECT(w(j)) = .TRUE. after ordering, since ordering may
! 	    change the value of complex eigenvalues (especially if the
! 	    eigenvalue is ill-conditioned). In this case INFO is set to
! 	    size(A,1) + 2 (see INFO below).
!          Note: Select must be present if SDIM is desired.
! SDIM     Optional (output) INTEGER.
!          The number of eigenvalues (after sorting) for which
! 	   SELECT=.TRUE. (If A is real, complex conjugate pairs for which
! 	   SELECT=.TRUE. for either eigenvalue count as 2).
! INFO     Optional (output) INTEGER.
!          = 0: successful exit.
!          < 0: if INFO = -i, the i-th argument had an illegal value.
!          > 0: if INFO = i, and i is
!              <= n: the QR algorithm failed to compute all the
!                    eigenvalues; elements 1:ilo-1 and i+1:n of w contain
!                    those eigenvalues which have converged. VS contains
!                    the matrix which reduces A to its partially converged
!                    Schur form.
!              = n+1: the eigenvalues could not be reordered because some
!                    eigenvalues were not sufficiently separated (the
! 		   problem is very ill-conditioned).
!              = n+2: after reordering, some leading complex eigenvalues
! 	           in the Schur form no longer satisfy SELECT = .TRUE.
!                    This can be caused by ordinary roundoff or underflow
!                    due to scaling.
!              n is the order of A.
!          If INFO is not present and an error occurs, then the program is
!          terminated with an error message.
!-------------------------------------------------------------------------
!  .. LOCAL PARAMETERS ..
   CHARACTER(LEN=7), PARAMETER :: SRNAME = 'LA_GEES'
!  .. LOCAL SCALARS ..
   CHARACTER(LEN=1) :: LJOBVS, LSORT
   INTEGER, SAVE :: LWORK = 0
   INTEGER :: N, LINFO, LD, ISTAT, ISTAT1, S1VS, S2VS, LSDIM
!  .. LOCAL ARRAYS ..
   LOGICAL, TARGET :: LLBWORK(1)
   LOGICAL, POINTER :: BWORK(:)
   REAL(WP), TARGET :: LLVS(1,1)
   REAL(WP), POINTER :: WORK(:)
!  .. INTRINSIC FUNCTIONS ..
   INTRINSIC MAX, PRESENT, SIZE
!  .. EXECUTABLE STATEMENTS ..
   LSDIM = 0
   LINFO = 0; ISTAT = 0; N = SIZE(A,1); LD = MAX(1,N)
   IF( PRESENT(VS) )THEN; S1VS = SIZE(VS,1); S2VS = SIZE(VS,2); LJOBVS = 'V'
   ELSE; S1VS = 1; S2VS = 1; LJOBVS = 'N'; END IF
   IF( PRESENT(SDIM) .OR. PRESENT(SELECT) )THEN; LSORT = 'S'; ELSE; LSORT = 'N'; END IF
!  .. TEST THE ARGUMENTS
   IF( N < 0 .OR. SIZE(A,2) /= N )THEN; LINFO = -1
   ELSE IF( SIZE( WR ) /= N )THEN; LINFO = -2
   ELSE IF( SIZE( WI ) /= N )THEN; LINFO = -3
   ELSE IF( PRESENT(VS) .AND. ( S1VS /= N .OR. S2VS /= N ) )THEN; LINFO = -4
   ELSE IF( PRESENT(SDIM) .AND. .NOT.PRESENT(SELECT) )THEN; LINFO = -5
   ELSE IF( N > 0 )THEN
      IF( ISTAT == 0 ) THEN
         LWORK = MAX( 1, 3*N, LWORK); ALLOCATE(WORK(LWORK), STAT=ISTAT)
         IF( ISTAT /= 0 )THEN; DEALLOCATE(WORK,STAT=ISTAT1)
            LWORK = MAX( 1, 3*N ); ALLOCATE(WORK(LWORK), STAT=ISTAT)
            IF( ISTAT == 0) CALL ERINFO( -200, SRNAME, LINFO )
         END IF
      END IF
      IF( ISTAT == 0 ) THEN
        IF( PRESENT(VS) )THEN
          IF( PRESENT(SDIM) )THEN
	     ALLOCATE(BWORK(N),STAT=ISTAT)
             CALL GEES_F77( LJOBVS, LSORT, SELECT, N, A, LD, LSDIM, WR, WI, &
                        VS, S1VS, WORK, LWORK, BWORK, LINFO )
	    ELSE
	     CALL GEES_F77( LJOBVS, LSORT, SELECT, N, A, LD, LSDIM, WR, WI, &
                        VS, S1VS, WORK, LWORK, LLBWORK, LINFO )
	    ENDIF
	  ELSE
	   IF( PRESENT(SDIM) )THEN
	      ALLOCATE(BWORK(N),STAT=ISTAT)
	      CALL GEES_F77( LJOBVS, LSORT, SELECT, N, A, LD, LSDIM, WR, WI, &
                        LLVS, S1VS, WORK, LWORK, BWORK, LINFO )
	   ELSE
	      CALL GEES_F77( LJOBVS, LSORT, SELECT, N, A, LD, LSDIM, WR, WI, &
                        LLVS, S1VS, WORK, LWORK, LLBWORK, LINFO )
	  ENDIF
	 ENDIF
         IF( LINFO == 0 ) LWORK = INT(WORK(1)+1)
      ELSE; LINFO = -100; ENDIF
!      IF( PRESENT(SDIM) ) SDIM = LSDIM
      IF( LSAME(LSORT,'S') ) DEALLOCATE(BWORK, STAT=ISTAT1)
      DEALLOCATE(WORK, STAT=ISTAT1)
   ENDIF
   IF( PRESENT(SDIM) ) SDIM = LSDIM
   CALL ERINFO(LINFO,SRNAME,INFO,ISTAT)
END SUBROUTINE SGEES_F95

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


text:
program test
  integer :: a
  integer :: b

  a = 6 + 3 + 4 + 4 +  &
    3 + 6
  b = 7 + 6 + 6 + 5   &
    + 5 + 6
  print *, a + b
end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
c---- ANOMALOUS COUPLINGS AVAILABLE FOR SPIN 0 HIGGS
c---- g1: SM, g2: Higher-order scalar, g3: Form factor dependence on g2, g4: Pseudoscalar
c---- prime: Dipole ansatz product for q1,2**2 dependence
c---- prime 2: |q1**2| + |q2**2|
c---- prime 3: |q1**2| - |q2**2|
c---- prime 4: (q1 + q2)**2
c---- prime 5: q1**4 + q2**4
c---- prime 6: q1**4 - q2**4
c---- prime 7: |q1**2| * |q2**2|
c---- NOTE: Please add new future couplings in the same order for both declaration and the common block

      integer AllowAnomalousCouplings
      integer distinguish_HWWcouplings
      integer AnomalCouplPR
      integer AnomalCouplDK
      integer channeltoggle_stu ! 0, 1, 2 for s, t+u and s+t+u
      integer vvhvvtoggle_vbfvh ! 0, 1, 2 for VBF, VH and VB+VH

      integer cz_q1sq,cz_q2sq,cz_q12sq
      integer cw_q1sq,cw_q2sq,cw_q12sq
      integer c2z_q1sq,c2z_q2sq,c2z_q12sq
      integer c2w_q1sq,c2w_q2sq,c2w_q12sq

      double precision mb_4gen,mt_4gen

      double precision LambdaBSM,Lambda_Q
      double precision Lambda_zgs1
      double precision Lambda_z1,Lambda_z2,Lambda_z3,Lambda_z4
      double precision Lambda_z11,Lambda_z21,Lambda_z31,Lambda_z41
      double precision Lambda_z12,Lambda_z22,Lambda_z32,Lambda_z42
      double precision Lambda_z10,Lambda_z20,Lambda_z30,Lambda_z40
      double precision Lambda_w1,Lambda_w2,Lambda_w3,Lambda_w4
      double precision Lambda_w11,Lambda_w21,Lambda_w31,Lambda_w41
      double precision Lambda_w12,Lambda_w22,Lambda_w32,Lambda_w42
      double precision Lambda_w10,Lambda_w20,Lambda_w30,Lambda_w40

      double precision h2mass,h2width

      double precision Lambda2BSM,Lambda2_Q
      double precision Lambda2_zgs1
      double precision Lambda2_z1,Lambda2_z2,Lambda2_z3,Lambda2_z4
      double precision Lambda2_z11,Lambda2_z21,Lambda2_z31,Lambda2_z41
      double precision Lambda2_z12,Lambda2_z22,Lambda2_z32,Lambda2_z42
      double precision Lambda2_z10,Lambda2_z20,Lambda2_z30,Lambda2_z40
      double precision Lambda2_w1,Lambda2_w2,Lambda2_w3,Lambda2_w4
      double precision Lambda2_w11,Lambda2_w21,Lambda2_w31,Lambda2_w41
      double precision Lambda2_w12,Lambda2_w22,Lambda2_w32,Lambda2_w42
      double precision Lambda2_w10,Lambda2_w20,Lambda2_w30,Lambda2_w40

      double complex kappa_top,kappa_tilde_top
      double complex kappa_bot,kappa_tilde_bot
      double complex ghg2,ghg3,ghg4
      double complex kappa_4gen_top,kappa_tilde_4gen_top
      double complex kappa_4gen_bot,kappa_tilde_4gen_bot
      double complex ghg2_4gen,ghg3_4gen,ghg4_4gen

      double complex ghz1,ghz2,ghz3,ghz4
      double complex ghz1_prime,ghz2_prime,ghz3_prime,ghz4_prime
      double complex ghz1_prime2,ghz2_prime2,ghz3_prime2,ghz4_prime2
      double complex ghz1_prime3,ghz2_prime3,ghz3_prime3,ghz4_prime3
      double complex ghz1_prime4,ghz2_prime4,ghz3_prime4,ghz4_prime4
      double complex ghz1_prime5,ghz2_prime5,ghz3_prime5,ghz4_prime5
      double complex ghz1_prime6,ghz2_prime6,ghz3_prime6,ghz4_prime6
      double complex ghz1_prime7,ghz2_prime7,ghz3_prime7,ghz4_prime7

      double complex ghzgs1_prime2,ghzgs2,ghzgs3,ghzgs4
      double complex ghgsgs2,ghgsgs3,ghgsgs4

      double complex ghw1,ghw2,ghw3,ghw4
      double complex ghw1_prime,ghw2_prime,ghw3_prime,ghw4_prime
      double complex ghw1_prime2,ghw2_prime2,ghw3_prime2,ghw4_prime2
      double complex ghw1_prime3,ghw2_prime3,ghw3_prime3,ghw4_prime3
      double complex ghw1_prime4,ghw2_prime4,ghw3_prime4,ghw4_prime4
      double complex ghw1_prime5,ghw2_prime5,ghw3_prime5,ghw4_prime5
      double complex ghw1_prime6,ghw2_prime6,ghw3_prime6,ghw4_prime6
      double complex ghw1_prime7,ghw2_prime7,ghw3_prime7,ghw4_prime7

      double complex kappa2_top,kappa2_tilde_top
      double complex kappa2_bot,kappa2_tilde_bot
      double complex gh2g2,gh2g3,gh2g4
      double complex kappa2_4gen_top,kappa2_tilde_4gen_top
      double complex kappa2_4gen_bot,kappa2_tilde_4gen_bot
      double complex gh2g2_4gen,gh2g3_4gen,gh2g4_4gen

      double complex gh2z1,gh2z2,gh2z3,gh2z4
      double complex gh2z1_prime,gh2z2_prime,gh2z3_prime,gh2z4_prime
      double complex gh2z1_prime2,gh2z2_prime2,gh2z3_prime2,gh2z4_prime2
      double complex gh2z1_prime3,gh2z2_prime3,gh2z3_prime3,gh2z4_prime3
      double complex gh2z1_prime4,gh2z2_prime4,gh2z3_prime4,gh2z4_prime4
      double complex gh2z1_prime5,gh2z2_prime5,gh2z3_prime5,gh2z4_prime5
      double complex gh2z1_prime6,gh2z2_prime6,gh2z3_prime6,gh2z4_prime6
      double complex gh2z1_prime7,gh2z2_prime7,gh2z3_prime7,gh2z4_prime7

      double complex gh2zgs1_prime2,gh2zgs2,gh2zgs3,gh2zgs4
      double complex gh2gsgs2,gh2gsgs3,gh2gsgs4

      double complex gh2w1,gh2w2,gh2w3,gh2w4
      double complex gh2w1_prime,gh2w2_prime,gh2w3_prime,gh2w4_prime
      double complex gh2w1_prime2,gh2w2_prime2,gh2w3_prime2,gh2w4_prime2
      double complex gh2w1_prime3,gh2w2_prime3,gh2w3_prime3,gh2w4_prime3
      double complex gh2w1_prime4,gh2w2_prime4,gh2w3_prime4,gh2w4_prime4
      double complex gh2w1_prime5,gh2w2_prime5,gh2w3_prime5,gh2w4_prime5
      double complex gh2w1_prime6,gh2w2_prime6,gh2w3_prime6,gh2w4_prime6
      double complex gh2w1_prime7,gh2w2_prime7,gh2w3_prime7,gh2w4_prime7

!     anomalous couplings for triple and quartic gauge boson couplings: VVV, VVVV
!     notation wrt paper: dV_N=d^{NWW}*d_1^N, dP_N=d^{NWW}*d_2^N, dM_N=d^{NWW}*d_3^N, dFour_N=d^{NWW}*d_4^N, for N=Z,\gamma
!                         dZZWpWm = d^{ZZWW}*sw^2/cw^2, dZAWpWm = d^{Z\gammaWW}*sw/cw, dAAWpWm = d^{\gamma\gammaWW}
      double complex dV_A,dP_A,dM_A,dFour_A
      double complex dV_Z,dP_Z,dM_Z,dFour_Z
      double complex dZZWpWm,dZAWpWm,dAAWpWm

      common/spinzerohiggs_anomcoupl/
     & AllowAnomalousCouplings,
     & distinguish_HWWcouplings,
     & AnomalCouplPR,AnomalCouplDK,
     & channeltoggle_stu,vvhvvtoggle_vbfvh,
     & cz_q1sq,cz_q2sq,cz_q12sq,
     & cw_q1sq,cw_q2sq,cw_q12sq,
     & c2z_q1sq,c2z_q2sq,c2z_q12sq,
     & c2w_q1sq,c2w_q2sq,c2w_q12sq,

     & mb_4gen,mt_4gen,

     & LambdaBSM,Lambda_Q,
     & Lambda_zgs1,
     & Lambda_z1,Lambda_z2,Lambda_z3,Lambda_z4,
     & Lambda_z11,Lambda_z21,Lambda_z31,Lambda_z41,
     & Lambda_z12,Lambda_z22,Lambda_z32,Lambda_z42,
     & Lambda_z10,Lambda_z20,Lambda_z30,Lambda_z40,
     & Lambda_w1,Lambda_w2,Lambda_w3,Lambda_w4,
     & Lambda_w11,Lambda_w21,Lambda_w31,Lambda_w41,
     & Lambda_w12,Lambda_w22,Lambda_w32,Lambda_w42,
     & Lambda_w10,Lambda_w20,Lambda_w30,Lambda_w40,

     & h2mass,h2width,

     & Lambda2BSM,Lambda2_Q,
     & Lambda2_zgs1,
     & Lambda2_z1,Lambda2_z2,Lambda2_z3,Lambda2_z4,
     & Lambda2_z11,Lambda2_z21,Lambda2_z31,Lambda2_z41,
     & Lambda2_z12,Lambda2_z22,Lambda2_z32,Lambda2_z42,
     & Lambda2_z10,Lambda2_z20,Lambda2_z30,Lambda2_z40,
     & Lambda2_w1,Lambda2_w2,Lambda2_w3,Lambda2_w4,
     & Lambda2_w11,Lambda2_w21,Lambda2_w31,Lambda2_w41,
     & Lambda2_w12,Lambda2_w22,Lambda2_w32,Lambda2_w42,
     & Lambda2_w10,Lambda2_w20,Lambda2_w30,Lambda2_w40,

     & kappa_top,kappa_tilde_top,
     & kappa_bot,kappa_tilde_bot,
     & ghg2,ghg3,ghg4,
     & kappa_4gen_top,kappa_tilde_4gen_top,
     & kappa_4gen_bot,kappa_tilde_4gen_bot,
     & ghg2_4gen,ghg3_4gen,ghg4_4gen,

     & ghz1,ghz2,ghz3,ghz4,
     & ghz1_prime,ghz2_prime,ghz3_prime,ghz4_prime,
     & ghz1_prime2,ghz2_prime2,ghz3_prime2,ghz4_prime2,
     & ghz1_prime3,ghz2_prime3,ghz3_prime3,ghz4_prime3,
     & ghz1_prime4,ghz2_prime4,ghz3_prime4,ghz4_prime4,
     & ghz1_prime5,ghz2_prime5,ghz3_prime5,ghz4_prime5,
     & ghz1_prime6,ghz2_prime6,ghz3_prime6,ghz4_prime6,
     & ghz1_prime7,ghz2_prime7,ghz3_prime7,ghz4_prime7,

     & ghzgs1_prime2,ghzgs2,ghzgs3,ghzgs4,
     & ghgsgs2,ghgsgs3,ghgsgs4,

     & ghw1,ghw2,ghw3,ghw4,
     & ghw1_prime,ghw2_prime,ghw3_prime,ghw4_prime,
     & ghw1_prime2,ghw2_prime2,ghw3_prime2,ghw4_prime2,
     & ghw1_prime3,ghw2_prime3,ghw3_prime3,ghw4_prime3,
     & ghw1_prime4,ghw2_prime4,ghw3_prime4,ghw4_prime4,
     & ghw1_prime5,ghw2_prime5,ghw3_prime5,ghw4_prime5,
     & ghw1_prime6,ghw2_prime6,ghw3_prime6,ghw4_prime6,
     & ghw1_prime7,ghw2_prime7,ghw3_prime7,ghw4_prime7,

     & kappa2_top,kappa2_tilde_top,
     & kappa2_bot,kappa2_tilde_bot,
     & gh2g2,gh2g3,gh2g4,
     & kappa2_4gen_top,kappa2_tilde_4gen_top,
     & kappa2_4gen_bot,kappa2_tilde_4gen_bot,
     & gh2g2_4gen,gh2g3_4gen,gh2g4_4gen,

     & gh2z1,gh2z2,gh2z3,gh2z4,
     & gh2z1_prime,gh2z2_prime,gh2z3_prime,gh2z4_prime,
     & gh2z1_prime2,gh2z2_prime2,gh2z3_prime2,gh2z4_prime2,
     & gh2z1_prime3,gh2z2_prime3,gh2z3_prime3,gh2z4_prime3,
     & gh2z1_prime4,gh2z2_prime4,gh2z3_prime4,gh2z4_prime4,
     & gh2z1_prime5,gh2z2_prime5,gh2z3_prime5,gh2z4_prime5,
     & gh2z1_prime6,gh2z2_prime6,gh2z3_prime6,gh2z4_prime6,
     & gh2z1_prime7,gh2z2_prime7,gh2z3_prime7,gh2z4_prime7,

     & gh2zgs1_prime2,gh2zgs2,gh2zgs3,gh2zgs4,
     & gh2gsgs2,gh2gsgs3,gh2gsgs4,

     & gh2w1,gh2w2,gh2w3,gh2w4,
     & gh2w1_prime,gh2w2_prime,gh2w3_prime,gh2w4_prime,
     & gh2w1_prime2,gh2w2_prime2,gh2w3_prime2,gh2w4_prime2,
     & gh2w1_prime3,gh2w2_prime3,gh2w3_prime3,gh2w4_prime3,
     & gh2w1_prime4,gh2w2_prime4,gh2w3_prime4,gh2w4_prime4,
     & gh2w1_prime5,gh2w2_prime5,gh2w3_prime5,gh2w4_prime5,
     & gh2w1_prime6,gh2w2_prime6,gh2w3_prime6,gh2w4_prime6,
     & gh2w1_prime7,gh2w2_prime7,gh2w3_prime7,gh2w4_prime7,

     & dV_A,dP_A,dM_A,dFour_A,
     & dV_Z,dP_Z,dM_Z,dFour_Z,
     & dZZWpWm,dZAWpWm,dAAWpWm

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
        function pct_rejected(n_good_qc,n_bad_qc)

        n_total = n_good_qc + n_bad_qc

        if(n_total .gt. 0)then
            pct_rejected = (float(n_bad_qc) / float(n_total)) * 100.
        else
            pct_rejected = 0.
        endif

        return
        end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense