text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine MACNUT     ( pmsa   , fl     , ipoint , increm, noseg , &
                              noflux , iexpnt , iknmrk , noq1  , noq2  , &
                              noq3   , noq4   )
!
!*******************************************************************************
!
      implicit none
!
!     type    name         i/o description
!
      real(4) pmsa(*)     !i/o process manager system array, window of routine to process library
      real(4) fl(*)       ! o  array of fluxes made by this process in mass/volume/time
      integer ipoint( 43) ! i  array of pointers in pmsa to get and store the data
      integer increm( 43) ! i  increments in ipoint for segment loop, 0=constant, 1=spatially varying
      integer noseg       ! i  number of computational elements in the whole model schematisation
      integer noflux      ! i  number of fluxes, increment in the fl array
      integer iexpnt(4,*) ! i  from, to, from-1 and to+1 segment numbers of the exchange surfaces
      integer iknmrk(*)   ! i  active-inactive, surface-water-bottom, see manual for use
      integer noq1        ! i  nr of exchanges in 1st direction (the horizontal dir if irregular mesh)
      integer noq2        ! i  nr of exchanges in 2nd direction, noq1+noq2 gives hor. dir. reg. grid
      integer noq3        ! i  nr of exchanges in 3rd direction, vertical direction, pos. downward
      integer noq4        ! i  nr of exchanges in the bottom (bottom layers, specialist use only)
      integer ipnt( 43)   !    local work array for the pointering
      integer iseg        !    local loop counter for computational element loop
!
!*******************************************************************************
!
!     type    name         i/o description                                        unit
!
      real(4) depth       ! i  depth of segment                                   (m)
      real(4) totaldepth  ! i  total depth water column                           (m)
      real(4) locseddept  ! i  sediment layer depth to bottom of segment          (m)
      real(4) po4s12      ! i  concentration of PO4 in S12 bottom (pores)         (g/m3)
      real(4) nh4s12      ! i  concentration of NH4 in S12 bottom (pores)         (g/m3)
      integer ibotseg     ! i  bottom segment number                              (-)
      real(4) frbmlay     ! i  fraction in layer of SM biomass in column          (-)
      real(4) rootdesm01  ! i  rooting depth sm01                                 (m)
      real(4) poros       ! i  volumetric porosity                                (-)
      real(4) nh4         ! i  ammonium (nh4)                                     (gn/m3)
      real(4) no3         ! i  nitrate (no3)                                      (gn/m3)
      real(4) po4         ! i  ortho-phosphate (po4)                              (gp/m3)
      real(4) disco2      ! i  concentration of dissolved carbon dioxide          (g/m3)
      real(4) dish2co3    ! i  concentration of dissolved true h2co3              (gc/m3)
      real(4) dishco3     ! i  concentration of dissolved hco3(-)                 (gc/m3)
      real(4) prfnh4sm01  ! i  ammonium preferency over nitrate sm01              (-)
      real(4) kmdinsm01w  ! i  half-saturation value n sm01 in water              (gn/m3)
      real(4) kmpsm01w    ! i  half-saturation value p sm01 in water              (gp/m3)
      real(4) kmco2sm01   ! i  half-saturation value co2 + h2c03 sm01             (gc/m3)
      real(4) kmhco3sm01  ! i  half-saturation value hc03 sm01                    (gc/m3)
      real(4) kmdinsm01b  ! i  half-saturation value n sm01 in bottom             (gn/m3)
      real(4) kmpsm01b    ! i  half-saturation value p sm01 in bottom             (gp/m3)
      real(4) limnsm01w   ! o  nitrogen limitation function sm01 <0-1> water      (-)
      real(4) limpsm01w   ! o  phosphorus limitation function sm01 <0-1> water    (-)
      real(4) lco2sm01    ! o  co2+h2c03 limitation function sm01 <0-1>           (-)
      real(4) limnsm01b   ! o  nitrogen limitation function sm01 <0-1> bottom     (-)
      real(4) limpsm01b   ! o  phosphorus limitation function sm01 <0-1>bottom    (-)
      real(4) limnutsm01  ! o  nutrient limitation function sm01 <0-1>            (-)
      real(4) frootuptn   ! o  fraction root uptake nitrogen sm01                 (-)
      real(4) frootuptp   ! o  fraction root uptake phosphorus                    (-)
      real(4) cdinsm01w   ! o  average water concentration din for sm01           (m)
      real(4) cpo4sm01w   ! o  average water concentration po4 for sm01           (m)
      real(4) cco2sm01    ! o  average concentration co2+h2co3 for sm01           (m)
      real(4) chco3sm01   ! o  average concentration hco3 for sm01                (m)
      real(4) cdinsm01b   ! o  average sediment concentration din for sm01        (m)
      real(4) cpo4sm01b   ! o  average sediment concentration po4 for sm01        (m)

      ! local

      integer ikmrk1      !    first attribute
      integer ikmrk2      !    second attribute
      real(4) din         ! l  dissolved inorganic nitrogen, corrected for preference
      real(4) z1          ! l  z1
      real(4) fr_avg      ! l  fr_avg
      real(4) hroot       ! l  hroot
      real(4) limn        ! l  limn
      real(4) limp        ! l  limp

      ! zero the average concentrations for all segments

      ipnt  = ipoint
      do iseg = 1 , noseg
         pmsa(ipnt(38)) = 0.0
         pmsa(ipnt(39)) = 0.0
         pmsa(ipnt(40)) = 0.0
         pmsa(ipnt(41)) = 0.0
         pmsa(ipnt(42)) = 0.0
         pmsa(ipnt(43)) = 0.0
         ipnt  = ipnt + increm
      enddo

      ! first loop average concentration over height and rooting depth

      ipnt  = ipoint
      do iseg = 1 , noseg

         depth      = pmsa( ipnt(  1) )
         totaldepth = pmsa( ipnt(  2) )
         locseddept = pmsa( ipnt(  4) )
         nh4s12     = pmsa( ipnt(  5) )
         po4s12     = pmsa( ipnt(  6) )
         ibotseg    = nint(pmsa( ipnt(  7) ))
         FrBmLay    = pmsa( ipnt(  8) )
         rootdesm01 = pmsa( ipnt(  9) )
         poros      = pmsa( ipnt( 10) )
         nh4        = max(pmsa( ipnt( 11) ),0.0)
         no3        = max(pmsa( ipnt( 12) ),0.0)
         po4        = max(pmsa( ipnt( 13) ),0.0)
         disco2     = max(pmsa( ipnt( 14) ),0.0)
         dish2co3   = max(pmsa( ipnt( 15) ),0.0)
         dishco3    = max(pmsa( ipnt( 16) ),0.0)
         prfnh4sm01 = pmsa( ipnt( 17) )

         ! convert co2 to carbon and add h2co3, calculated din with preference
         ! adjust for porosity

         din        = (nh4 + no3/prfnh4sm01)/poros
         po4        =  po4/poros
         disco2     = (disco2*12./44. + dish2co3)/poros
         dishco3    =  dishco3/poros

         call dhkmrk(1,iknmrk(iseg),ikmrk1)
         if (ikmrk1.eq.1) then

            ! active water segment

            fr_avg = FrBmLay

            pmsa(botidx(38)) = pmsa(botidx(38)) + din*fr_avg
            pmsa(botidx(39)) = pmsa(botidx(39)) + po4*fr_avg
            pmsa(botidx(40)) = pmsa(botidx(40)) + disco2*fr_avg
            pmsa(botidx(41)) = pmsa(botidx(41)) + dishco3*fr_avg

            ! S12 sediment concentration

            call dhkmrk(2,iknmrk(iseg),ikmrk2)
            if ((ikmrk2.eq.0).or.(ikmrk2.eq.3)) then
              if (nh4s12.gt.0.0) pmsa(botidx(42)) = nh4s12
              if (po4s12.gt.0.0) pmsa(botidx(43)) = po4s12
            endif

         elseif (ikmrk1.eq.3) then

            ! sediment bed segment, distribution of roots in bed

            hroot = min(rootdesm01,totaldepth)
            z1 = locseddept - depth

            if (hroot .gt. locseddept) then
               ! completely in segment:
               fr_avg = min(1.0,depth/hroot)
            elseif (hroot .gt. z1 ) then
               ! partialy in segment:
               fr_avg = (hroot-z1)/hroot
            else
               ! not in segment:
               fr_avg = 0.0
            endif

            pmsa(botidx(42)) = pmsa(botidx(42)) + din*fr_avg
            pmsa(botidx(43)) = pmsa(botidx(43)) + po4*fr_avg

         endif

         ipnt = ipnt + increm

      enddo

      ! second loop limiting factors and water / sediment uptake ratio

      ipnt  = ipoint
      do iseg = 1 , noseg

         call dhkmrk(1,iknmrk(iseg),ikmrk1)
         if (ikmrk1.eq.1) then
            call dhkmrk(2,iknmrk(iseg),ikmrk2)
            if ((ikmrk2.eq.0).or.(ikmrk2.eq.3)) then

               kmdinsm01w = pmsa( ipnt( 18) )
               kmpsm01w   = pmsa( ipnt( 19) )
               kmco2sm01  = pmsa( ipnt( 20) )
               kmhco3sm01 = pmsa( ipnt( 21) )
               kmdinsm01b = pmsa( ipnt( 22) )
               kmpsm01b   = pmsa( ipnt( 23) )
               cdinsm01w  = pmsa( ipnt( 24) )
               cpo4sm01w  = pmsa( ipnt( 25) )
               cco2sm01   = pmsa( ipnt( 26) )
               chco3sm01  = pmsa( ipnt( 27) )
               cdinsm01b  = pmsa( ipnt( 28) )
               cpo4sm01b  = pmsa( ipnt( 29) )

               ! n limitation

               if ( kmdinsm01b .lt. 1e-20 ) then
                  ! only in water
                  limnsm01w  = cdinsm01w/(cdinsm01w+kmdinsm01w)
                  limnsm01b  = -1.
               else
                  if ( kmdinsm01w .lt. 1e-20 ) then
                     ! only in bottom
                     limnsm01b  = cdinsm01b/(cdinsm01b+kmdinsm01b)
                     limnsm01w  = -1.
                  else
                     ! both
                     limnsm01w  = cdinsm01w/(cdinsm01w+kmdinsm01w)
                     limnsm01b  = cdinsm01b/(cdinsm01b+kmdinsm01b)
                  endif
               endif
               limn = max(limnsm01w,limnsm01b)
               if ( cdinsm01w .gt. 1e-10 ) then
                  if ( cdinsm01b .gt. 1e-10 ) then
                     frootuptn = .998/(1.+2.66*(cdinsm01b/cdinsm01w)**-0.83)
                  else
                     frootuptn = 0.0
                  endif
               else
                  frootuptn = 1.0
               endif

               ! p limitation

               if ( kmpsm01b .lt. 1e-20 ) then
                  ! only in water
                  limpsm01w  = cpo4sm01w/(cpo4sm01w+kmpsm01w)
                  limpsm01b  = -1.
               else
                  if ( kmpsm01w .lt. 1e-20 ) then
                     ! only in bottom
                     limpsm01b  = cpo4sm01b/(cpo4sm01b+kmpsm01b)
                     limpsm01w  = -1.
                  else
                     ! both
                     limpsm01w  = cpo4sm01w/(cpo4sm01w+kmpsm01w)
                     limpsm01b  = cpo4sm01b/(cpo4sm01b+kmpsm01b)
                  endif
               endif
               limp = max(limpsm01w,limpsm01b)
               if ( cpo4sm01w .gt. 1e-10 ) then
                  if ( cpo4sm01b .gt. 1e-10 ) then
                     frootuptp = .998/(1.+2.66*(cpo4sm01b/cpo4sm01w)**-0.83)
                  else
                     frootuptp = 0.0
                  endif
               else
                  frootuptp = 1.0
               endif

               ! c limitation

               if ( kmco2sm01 .lt. 1e-20 ) then
                  if ( kmhco3sm01 .lt. 1e-20 ) then
                     ! no c limitation calculation
                     lco2sm01 = 1.0
                  else
                     ! only hco3 limitation calculatation
                     lco2sm01 = chco3sm01/(chco3sm01+kmhco3sm01)
                  endif
               else
                  if ( kmhco3sm01 .lt. 1e-20 ) then
                     ! only co2 limitation calculatation
                     lco2sm01 = cco2sm01/(cco2sm01+kmco2sm01)
                  else
                     ! both co2 and hco3 limitation calculatation, take the max
                     lco2sm01 = max(cco2sm01/(cco2sm01+kmco2sm01),chco3sm01/(chco3sm01+kmhco3sm01))
                  endif
               endif

               ! overall limitation is minimum of n, p and c

               limnutsm01 = min(limn,limp,lco2sm01)

               pmsa( ipnt( 30)   ) = limnsm01w
               pmsa( ipnt( 31)   ) = limpsm01w
               pmsa( ipnt( 32)   ) = lco2sm01
               pmsa( ipnt( 33)   ) = limnsm01b
               pmsa( ipnt( 34)   ) = limpsm01b
               pmsa( ipnt( 35)   ) = limnutsm01
               pmsa( ipnt( 36)   ) = frootuptn
               pmsa( ipnt( 37)   ) = frootuptp

            endif
         endif
         ipnt        = ipnt        + increm

      enddo

      return
      contains

      ! Auxiliary function to compute the index for the associated bottom segment
      integer function botidx(number)
         integer, intent(in) :: number

         botidx = ipoint(number)+(ibotseg-1)*increm(number)
      end function botidx

      end subroutine

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      subroutine backtracking(dim,x,m,lambda,rho,equatn,linear,f,d,gtd,
     +alpha,fp,xp,evalaldim,setpdim,btinfo,inform)

C     SCALAR ARGUMENTS
      integer btinfo,dim,inform,m
      double precision alpha,f,fp,gtd

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision d(dim),lambda(m),rho(m),x(dim),xp(*)

C     SUBROUTINE ARGUMENTS
      external evalaldim,setpdim

C     Backtracking with quadratic interpolation.
C
C     btinfo:
C
C     0: Armijo satisfied
C     2: Unbounded objective function?
C     3: Too small backtracking step. Wrong gradient?

      include "dim.par"
      include "machconst.inc"
      include "algconst.par"
      include "counters.inc"
      include "outtyp.inc"

C     LOCAL SCALARS
      logical samep
      integer i,interp
      double precision atmp

      interp = 0

 2010 continue

C     Test Armijo condition

      if ( fp .le. f + alpha * gamma * gtd ) then

C         Finish backtracking with the current point

          btinfo = 0

          if ( iprintinn .ge. 6 ) then
C              write(*, 900)
              write(10,900)
          end if

          return

      end if

C     Test f going to -inf

      if ( fp .le. fmin ) then

C         Finish backtracking with the current point

          btinfo = 2

          if ( iprintinn .ge. 6 ) then
C              write(*, 920)
              write(10,920)
          end if

          return

      end if

C     Test if we obtained a functional value similar to the current one
C     associated to a very small step

      samep = .true.
      do i = 1,dim
          if ( xp(i) .gt. x(i) + macheps23 * max(1.0d0,abs(x(i))) .or.
     +         xp(i) .lt. x(i) - macheps23 * max(1.0d0,abs(x(i))) ) then
              samep = .false.
          end if
      end do

      if ( samep .and. fp .le. f + macheps23 * abs(f) ) then

C         Finish backtracking with the current point

          btinfo = 3

          if ( iprintinn .ge. 6 ) then
C              write(*, 930)
              write(10,930)
          end if

          return

      end if

C     Compute new step

      interp = interp + 1

      atmp = ( - gtd * alpha ** 2 ) /
     +       ( 2.0d0 * ( fp - f - alpha * gtd ) )

      if ( atmp .ge. sigma1 * alpha .and.
     +     atmp .le. sigma2 * alpha ) then
          alpha = atmp
      else
          alpha = alpha / etaint
      end if

C     Compute new trial point

      do i = 1,dim
          xp(i) = x(i) + alpha * d(i)
      end do

      call setpdim(dim,xp)

      call evalaldim(dim,xp,m,lambda,rho,equatn,linear,fp,inform)
      if ( inform .lt. 0 ) return

C     Print information of this iteration

      if ( iprintinn .ge. 6 ) then
C          write(*, 110) alpha,fp,fcnt
          write(10,110) alpha,fp,fcnt
      end if

      go to 2010

C     NON-EXECUTABLE STATEMENTS

 110  format(  5X,'Alpha = ',1P,D7.1,' F = ',1P,D24.16,' FE = ',I7)
 900  format(  5X,'Flag of backtracking: Armijo condition holds.')
 920  format(  5X,'Flag of backtracking: Unbounded objective function?')
 930  format(  5X,'Flag of backtracking: Too small backtracking step.')

      end

ext: f


text:
: alias ( xt -- ) create , does> @ execute ;
: noop ( -- ) ;
: defer ( <name> -- ) ['] noop alias ;
: is ( xt <name> -- ) ' >body ! ;

max_stars_repo_name: hanshuebner/rekonstrukt

max_issues_repo_name: hanshuebner/rekonstrukt

max_forks_repo_name: hanshuebner/rekonstrukt


text:
program main

  ! These modules are all in BoxLib/Src/F_BaseLib
  use boxlib
  use multifab_module
  use bl_IO_module

  ! These modules are all in this directory
  use init_data_module
  use write_plotfile_module
  use advance_module

  implicit none

  ! stuff you can set with the inputs file (otherwise use default values below)
  integer :: dim, nsteps, plot_int, n_cell, max_grid_size, verbose

  ! dummy indices using for reading in inputs file
  integer :: un, farg, narg
  logical :: need_inputs_file, found_inputs_file
  character(len=128) :: inputs_file_name

  integer, allocatable :: lo(:), hi(:)
  integer :: istep

  real(dp_t), allocatable :: prob_lo(:), prob_hi(:)
  real(dp_t) :: dx, dt, time, start_time, run_time, run_time_IOproc

  logical, allocatable :: is_periodic(:)

  type(box)      :: bx
  type(boxarray) :: ba
  type(layout)   :: la
  type(multifab) :: data

  namelist /probin/ dim, nsteps, plot_int, n_cell, max_grid_size, verbose

  ! if running in parallel, this will print out the number of MPI
  ! processes and OpenMP threads
  call boxlib_initialize()

  ! parallel_wtime() returns the number of wallclock-time seconds since
  ! the program began
  start_time = parallel_wtime()

  ! default values - will get overwritten by the inputs file
  dim           = 2
  nsteps        = 10
  plot_int      = 1
  n_cell        = 64
  max_grid_size = 32
  verbose       = 0

  ! read inputs file and overwrite any default values
  narg = command_argument_count()
  need_inputs_file = .true.
  farg = 1
  if ( need_inputs_file .AND. narg >= 1 ) then
     call get_command_argument(farg, value = inputs_file_name)
     inquire(file = inputs_file_name, exist = found_inputs_file )
     if ( found_inputs_file ) then
        farg = farg + 1
        un = unit_new()
        open(unit=un, file = inputs_file_name, status = 'old', action = 'read')
        read(unit=un, nml = probin)
        close(unit=un)
        need_inputs_file = .false.
     end if
  end if

  ! now that we have dim, we can allocate these
  allocate(lo(dim),hi(dim))
  allocate(is_periodic(dim))
  allocate(prob_lo(dim),prob_hi(dim))

  ! physical problem is a box on (-1,-1) to (1,1), periodic on all sides
  prob_lo(:) = -1.d0
  prob_hi(:) =  1.d0
  is_periodic(:) = .true.

  ! create a box from (0,0) to (n_cell-1,n_cell-1)
  lo(:) = 0
  hi(:) = n_cell-1
  bx = make_box(lo,hi)

  ! the grid spacing is the same in each direction
  dx = (prob_hi(1)-prob_lo(1)) / n_cell

  ! initialize the boxarray to be one single box
  call boxarray_build_bx(ba,bx)

  ! overwrite the boxarray to respect max_grid_size
  call boxarray_maxsize(ba,max_grid_size)

  ! build the layout, la
  ! the third argument is the problem domain, which in this case is bx
  call layout_build_ba(la,ba,bx,pmask=is_periodic)

  call destroy(ba)

  ! build multifab with 2 components and 6 ghost cells
  call multifab_build(data,la,2,6)

  ! initialze data
  call init_data(data,dx,prob_lo)

  istep = 0
  time = 0.d0

  dt = 0.1d0*dx

  ! write out plotfile 0
  call write_plotfile(la,data,istep,dx,time,prob_lo,prob_hi)

  do istep=1,nsteps

     if (verbose .ge. 1) then
        ! we only want one processor to write to screen
        if ( parallel_IOProcessor() ) then
           print*,'Advancing time step',istep
        end if
     end if

     ! advance the data
     call advance(data,dx,dt)

     time = time + dt

     if (mod(istep,plot_int) .eq. 0 .or. istep .eq. nsteps) then
        ! write out plotfile
        call write_plotfile(la,data,istep,dx,time,prob_lo,prob_hi)
     end if

  end do

  ! make sure to destroy the multifab or you'll leak memory
  call destroy(data)
  call destroy(la)

  deallocate(lo,hi,is_periodic,prob_lo,prob_hi)

  ! deallocate temporary boxarrays and communication mappings
  call layout_flush_copyassoc_cache()

  ! check for memory that should have been deallocated
  if ( parallel_IOProcessor() ) then
     print*, 'MEMORY STATS AT END OF PROGRAM'
     print*, ' '
  end if
  call print(multifab_mem_stats(),    "    multifab")
  call print(fab_mem_stats(),         "         fab")
  call print(boxarray_mem_stats(),    "    boxarray")
  call print(layout_mem_stats(),      "      layout")
  call print(boxassoc_mem_stats(),    "    boxassoc")
  call print(fgassoc_mem_stats(),     "     fgassoc")
  call print(syncassoc_mem_stats(),   "   syncassoc")
  call print(copyassoc_mem_stats(),   "   copyassoc")
  call print(fluxassoc_mem_stats(),   "   fluxassoc")

  ! parallel_wtime() returns the number of wallclock-time seconds since
  ! the program began
  run_time = parallel_wtime() - start_time

  ! collect run_time from each processor and store the maximum
  call parallel_reduce(run_time_IOproc, run_time, MPI_MAX, &
                       proc = parallel_IOProcessorNode())

  if ( parallel_IOProcessor() ) then
     print*,"Run time (s) =",run_time_IOproc
  end if

  call boxlib_finalize()

end program main

max_stars_repo_name: memmett/BoxLib

max_issues_repo_name: memmett/BoxLib

max_forks_repo_name: memmett/BoxLib


text:
      PROGRAM XTEST
C***********************************************************************
C*                                                                     *
C*  FORTRAN CODE WRITTEN FOR INCLUSION IN IBM RESEARCH REPORT RC20525, *
C*  'FORTRAN ROUTINES FOR USE WITH THE METHOD OF L-MOMENTS, VERSION 3' *
C*                                                                     *
C*  J. R. M. HOSKING                                                   *
C*  IBM RESEARCH DIVISION                                              *
C*  T. J. WATSON RESEARCH CENTER                                       *
C*  YORKTOWN HEIGHTS                                                   *
C*  NEW YORK 10598, U.S.A.                                             *
C*                                                                     *
C*  VERSION 3     AUGUST 1996                                          *
C*                                                                     *
C*  VERSION 3.03  JUNE 2000                                            *
C*  * CHARACTER variable declarations changed to conform with          *
C*    Fortran-77 standard                                              *
C*                                                                     *
C***********************************************************************
C
C  EXAMPLE PROGRAM TO ILLUSTRATE THE USE OF ROUTINE REGTST
C
C  THE ROUTINE READS IN THE SAMPLE L-MOMENTS FOR DIFFERENT SITES IN A
C  REGION AND CALLS ROUTINE REGTST TO COMPUTE DISCORDANCY, HETEROGENEITY
C  AND GOODNESS-OF-FIT STATISTICS.
C
C  PARAMETERS OF PROGRAM:
C  MAXNS  - MAXIMUM NUMBER OF SITES
C  SSEED  - SEED FOR RANDOM-NUMBER GENERATOR
C  NSIM   - NSIM PARAMETER OF ROUTINE REGTST
C  KPRINT - OUTPUT FLAG, KPRINT PARAMETER OF ROUTINE REGTST
C  INFILE - STREAM NUMBER TO WHICH INPUT FILE IS ATTACHED
C  KOUT   - STREAM NUMBER TO WHICH OUTPUT FILE IS ATTACHED
C  NPROB  - NUMBER OF QUANTILES TO BE FOUND
C  A      - ) PARAMETERS OF
C  B      - ) PLOTTING POSITION
C
C  ARRAYS TO BE INITIALIZED IN DATA STATEMENTS:
C  PROB   - PROBABILITIES FOR WHICH QUANTILES ARE TO BE FOUND
C
C  ROUTINES USED: REGTST AND ROUTINES CALLED BY REGTST
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (SSEED=619145091D0,NSIM=500,MAXNS=200)
      PARAMETER (KPRINT=1,INFILE=7,KOUT=6,NPROB=10)
C
      CHARACTER*12 NAMES(MAXNS)
      CHARACTER*60 REGNAM
      DOUBLE PRECISION D(MAXNS),H(3),PARA(5,6),PROB(NPROB),RMOM(5),
     *  VBAR(3),VOBS(3),VSD(3),XMOM(5,MAXNS),Z(5)
      INTEGER LEN(MAXNS)
C
      DATA A/0D0/,B/0D0/
      DATA PROB/0.01D0,0.02D0,0.05D0,0.1D0,0.2D0,
     *  0.5D0,0.90D0,0.95D0,0.99D0,0.999D0/
C
C         READ IN THE AT-SITE L-MOMENTS.
C         DATA FILE MAY CONTAIN ANY NUMBER OF REGIONAL DATA STRUCTURES.
C         A 'REGIONAL DATA STRUCTURE' CONSISTS OF:
C         1. ONE RECORD CONTAINING:
C            (COLUMNS  1- 4) NUMBER OF SITES IN REGION;
C            (COLUMNS  5-56) IDENTIFYING LABEL FOR THE REGION.
C         2. FOR EACH SITE, ONE RECORD CONTAINING:
C            (COLUMNS  1-12) AN IDENTIFYING LABEL FOR THE SITE;
C            (COLUMNS 13-16) THE RECORD LENGTH AT THE SITE;
C            (COLUMNS 17-56) SAMPLE L-MOMENTS L-1, L-CV, T-3, T-4, T-5.
C
    1 READ(INFILE,5000,END=900)NSITES,REGNAM
      WRITE(KOUT,6000)REGNAM,NSITES
      DO 10 ISITE=1,NSITES
      READ(INFILE,5010)NAMES(ISITE),LEN(ISITE),(XMOM(I,ISITE),I=1,5)
   10 CONTINUE
C
C         CALCULATE TEST STATISTICS
C
      SEED=SSEED
      CALL REGTST(NSITES,NAMES,LEN,XMOM,A,B,SEED,NSIM,NPROB,PROB,
     *  KPRINT,KOUT,RMOM,D,VOBS,VBAR,VSD,H,Z,PARA)
      GOTO 1
C
  900 CONTINUE
      WRITE(KOUT,6010)
      STOP
C
 5000 FORMAT(I4,A52)
 5010 FORMAT(A12,I4,5F8.0)
 6000 FORMAT(///1X,A52,I8,' SITES')
 6010 FORMAT(///' ALL DATA PROCESSED')
      END

ext: for

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
subroutine size(a,b)
        real*16 :: a,b
        end subroutine
program main
   REAL*16 :: x(2)
   call SIZE(x(1),x(2))
end program

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


text:
program main

  use boxlib
  use parallel
  use multifab_module
  use bl_IO_module
  use layout_module
  use init_data_module
  use write_plotfile_module
  use advance_module
  use bl_prof_module

  implicit none
  !
  ! We only support 3-D.
  !
  integer, parameter :: DM = 3
  !
  ! We need four grow cells.
  !
  integer, parameter :: NG  = 4
  !
  ! We have five components (no species).
  !
  integer, parameter :: NC  = 5

  integer            :: nsteps, plot_int, n_cell, max_grid_size
  integer            :: un, farg, narg
  logical            :: need_inputs_file, found_inputs_file
  character(len=128) :: inputs_file_name
  integer            :: i, lo(DM), hi(DM), istep
  double precision   :: prob_lo(DM), prob_hi(DM), cfl, eta, alam
  double precision   :: dx(DM), dt, time, start_time, end_time
  logical            :: is_periodic(DM)
  type(box)          :: bx
  type(boxarray)     :: ba
  type(layout)       :: la
  type(multifab)     :: U

  type(bl_prof_timer), save :: bpt, bpt_init_data, bpt_sleep

  !
  ! What's settable via an inputs file.
  !
  namelist /probin/ nsteps, plot_int, n_cell, max_grid_size, cfl, eta, alam

  call boxlib_initialize()
  call bl_prof_initialize(on = .true.)

  call build(bpt, "bpt_main")

  start_time = parallel_wtime()

  call build(bpt_sleep, "bpt_sleep")
  call Sleep(5)
  call destroy(bpt_sleep)

  !
  ! Namelist default values -- overwritable via inputs file.
  !
  nsteps        = 100
  plot_int      = 10
  n_cell        = 32
  max_grid_size = 32
  cfl           = 0.5d0
  eta           = 1.8d-4 ! Diffusion coefficient.
  alam          = 1.5d2  ! Diffusion coefficient.
  !
  ! Read inputs file and overwrite any default values.
  !
  narg = command_argument_count()
  need_inputs_file = .true.
  farg = 1
  if ( need_inputs_file .AND. narg >= 1 ) then
     call get_command_argument(farg, value = inputs_file_name)
     inquire(file = inputs_file_name, exist = found_inputs_file )
     if ( found_inputs_file ) then
        farg = farg + 1
        un = unit_new()
        open(unit=un, file = inputs_file_name, status = 'old', action = 'read')
        read(unit=un, nml = probin)
        close(unit=un)
        need_inputs_file = .false.
     end if
  end if

  if ( parallel_IOProcessor() ) then
     write(6,probin)
  end if
  !
  ! Physical problem is a box on (-1,-1) to (1,1), periodic on all sides.
  !
  prob_lo     = -0.1d0
  prob_hi     =  0.1d0
  is_periodic = .true.
  !
  ! Create a box from (0,0) to (n_cell-1,n_cell-1).
  !
  lo = 0
  hi = n_cell-1
  bx = make_box(lo,hi)

  do i = 1,DM
     dx(i) = (prob_hi(i)-prob_lo(i)) / n_cell
  end do

  call boxarray_build_bx(ba,bx)

  call boxarray_maxsize(ba,max_grid_size)

  call layout_build_ba(la,ba,boxarray_bbox(ba),pmask=is_periodic)

  call destroy(ba)

  call multifab_build(U,la,NC,NG)

  call build(bpt_init_data, "bpt_init_data")
  call init_data(U,dx,prob_lo,prob_hi)
  call destroy(bpt_init_data)

  istep = 0
  time  = 0.d0

  if (plot_int > 0) then
     call write_plotfile(U,istep,dx,time,prob_lo,prob_hi)
  end if

  do istep=1,nsteps

     if (parallel_IOProcessor()) then
        print*,'Advancing time step',istep,'time = ',time
     end if

     call advance(U,dt,dx,cfl,eta,alam)

     time = time + dt

     if (plot_int > 0) then
        if (mod(istep,plot_int) .eq. 0 .or. istep .eq. nsteps) then
           call write_plotfile(U,istep,dx,time,prob_lo,prob_hi)
        end if
     end if

  end do

  call destroy(U)
  call destroy(la)

  end_time = parallel_wtime()

  if ( parallel_IOProcessor() ) then
     print*,"Run time (s) =",end_time-start_time
  end if

  call destroy(bpt)
  call bl_prof_glean("bl_prof_report")
  call bl_prof_finalize()

  call boxlib_finalize()

end program main


text:
      SUBROUTINE PK_EQUATOR(KFILDO,IPACK,ND5,IS3,NS3,IPKOPT,L3264B,
     1                      LOCN,IPOS,IER,*)
C
C        MARCH   2000   LAWRENCE   GSC/TDL   ORIGINAL CODING
C        JANUARY 2001   GLAHN      COMMENTS; IER = 303 CHANGED TO 304;
C                                  ADDED CHECK FOR SIZE OF IS3( )
C        JANUARY 2001   GLAHN/LAWRENCE CHANGED IS3(45) TO ISC(47)
C                                  IN CALL TO PKBG
C
C        PURPOSE
C            PACKS TEMPLATE 3.110, AN EQUATORIAL AZIMUTHAL
C            EQUIDISTANT PROJECTION TEMPLATE, IN SECTION 3 OF A GRIB2.
C            IT IS THE RESPONSIBILITY OF THE CALLING ROUTINE TO PACK
C            THE FIRST 13 OCTETS IN SECTION 3.
C
C        DATA SET USE
C           KFILDO - UNIT NUMBER FOR OUTPUT (PRINT) FILE. (OUTPUT)
C
C        VARIABLES
C              KFILDO = UNIT NUMBER FOR OUTPUT (PRINT) FILE. (INPUT)
C            IPACK(J) = THE ARRAY THAT HOLDS THE ACTUAL PACKED MESSAGE
C                       (J=1,ND5). (INPUT/OUTPUT)
C                 ND5 = THE SIZE OF THE ARRAY IPACK( ). (INPUT)
C              IS3(J) = CONTAINS THE EQUATORIAL AZIMUTHAL EQUIDISTANT
C                       INFORMATION (IN ELEMENTS 15 THROUGH 57)
C                       THAT WILL BE PACKED INTO IPACK( ) (J=1,NS3).
C                       (INPUT/OUTPUT)
C                 NS3 = SIZE OF IS3( ). (INPUT)
C              IPKOPT = PACKING INDICATOR:
C                       0 = ERROR, DON'T PACK
C                       1 = PACK IA( ), SIMPLE
C                       2 = PACK IA( ) AND IB( ), SIMPLE
C                       3 = PACK COMPLEX OR SPATIAL DIFFERENCING
C                       4 = PACK COMPLEX.
C                       (INPUT)
C              L3264B = THE INTEGER WORD LENGTH IN BITS OF THE MACHINE
C                       BEING USED. VALUES OF 32 AND 64 ARE
C                       ACCOMMODATED. (INPUT)
C                LOCN = THE WORD POSITION TO PLACE THE NEXT VALUE.
C                       (INPUT/OUTPUT)
C                IPOS = THE BIT POSITION IN LOCN TO START PLACING
C                       THE NEXT VALUE. (INPUT/OUTPUT)
C                 IER = RETURN STATUS CODE. (OUTPUT)
C                         0 = GOOD RETURN.
C                       1-4 = ERROR CODES GENERATED BY PKBG. SEE THE
C                             DOCUMENTATION IN THE PKBG ROUTINE.
C                       302 = IS3( ) HAS NOT BEEN DIMENSIONED LARGE
C                             ENOUGH TO CONTAIN THE ENTIRE TEMPLATE.
C                       304 = IS3(13) DOES NOT INDICATE POLAR
C                             STEREOGRAPHIC MAP PROJECTION;
C                             INCORRECTLY CALLED SUBROUTINE.
C                       310 = UNRECOGNIZED OR UNSUPPORTED SHAPE OF
C                             EARTH CODE IN IS3(15) RETURNED BY
C                             SUBROUTINE EARTH.
C                   * = ALTERNATE RETURN WHEN IER .NE. 0.
C
C             LOCAL VARIABLES
C               ISIGN = A FLAG INDICATING WHETHER A VALUE TO BE
C                       PACKED IS POSITIVE OR NEGATIVE.
C                   N = L3264B = THE INTEGER WORD LENGTH IN BITS OF
C                       THE MACHINE BEING USED. VALUES OF 32 AND
C                       64 ARE ACCOMMODATED.
C             MINSIZE = THE SMALLEST ALLOWABLE DIMENSION FOR IS3( ).
C
C        NON SYSTEM SUBROUTINES CALLED
C           EARTH, PKBG
C
      PARAMETER(MINSIZE=57)
C
      DIMENSION IPACK(ND5),IS3(NS3)
C
      N=L3264B
      IER=0
C
C        CHECK TO MAKE SURE THAT THE USER WANTS TO PROCESS
C        A EQUATORIAL AZIMUTHAL EQUIDISTANT PROJECTION.  CHECK FOR
C        NS3 GE 13 ASSURES A VALUE FOR IS3(13).
C
      IF(NS3.GE.13)THEN
C
         IF (IS3(13).NE.110)THEN
D           WRITE(KFILDO,10)IS3(13)
D10         FORMAT(/' MAP PROJECTION CODE ',I4,' INDICATED BY IS3(13)'/
D    1              ' IS NOT EQUATORIAL AZIMUTHAL. PLEASE REFER'/
D    2              ' TO THE GRIB2 DOCUMENTATION TO DETERMINE THE'/
D    3              ' CORRECT MAP PROJECTION PACKER TO CALL.'/)
            IER=304
            GOTO 900
         ENDIF
C
      ELSE
         IER=302
         GO TO 900
      ENDIF
C
C        CHECK THE DIMENSIONS OF IS3( ).
C
      IF(NS3.LT.MINSIZE)THEN
         IER=302
         GO TO 900
      ENDIF
C
C        PACK THE SHAPE OF THE EARTH.
      CALL EARTH(KFILDO,IPACK,ND5,IS3,NS3,N,LOCN,IPOS,IER)
      IF(IER.NE.0)GO TO 900
C
C        PACK THE NUMBER OF POINTS ALONG THE X AND Y AXES
C        (Nx, Ny)
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(31),32,N,IER,*900)
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(35),32,N,IER,*900)
C
C        PACK THE LATITUDE & LONGITUDE OF THE TANGENCY POINT
      ISIGN=0
      IF(IS3(39).LT.0)ISIGN=1
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,ISIGN,1,N,IER,*900)
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,ABS(IS3(39)),
     1          31,N,IER,*900)
      ISIGN=0
      IF(IS3(43).LT.0)ISIGN=1
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,ISIGN,1,N,IER,*900)
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,ABS(IS3(43)),
     1          31,N,IER,*900)
C
C        PACK THE RESOLUTION AND COMPONENT FLAG.
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(47),8,N,IER,*900)
C
C        PACK Dx - X-DIRECTION GRID LENGTH IN UNITS OF 10-5 M
C        AS MEASURED AT THE POINT OF THE AXIS.
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(48),32,N,IER,*900)
C
C        PACK Dy - Y-DIRECTION GRID LENGTH IN UNITS OF 10-5 M
C        AS MEASURED AT THE POINT OF THE AXIS.
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(52),32,N,IER,*900)
C
C        PACK THE PROJECTION CENTER FLAG.
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(56),8,N,IER,*900)
C
C        PACK THE SCANNING MODE. CHECK TO SEE IF THE DATA POINTS
C        HAVE BEEN ORDERED BOUSTROPHEDONICALLY. IF SO, SET BIT 4
C        TO 1 IN THE SCANNING MODE OCTET.
C
      IF((IPKOPT.EQ.3).OR.(IPKOPT.EQ.4))THEN
         IS3(57)=IOR(IS3(57),16)
      ELSE
         IS3(57)=IAND(IS3(57),239)
      ENDIF
C
      CALL PKBG(KFILDO,IPACK,ND5,LOCN,IPOS,IS3(57),8,N,IER,*900)
C
C       ERROR RETURN SECTION
 900  IF(IER.NE.0)RETURN 1
C
      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! $UWHPSC/codes/fortran/sub3.f90

module sub3module

contains

subroutine fsub(x,f)
  ! compute f(x) = x**2 for all elements of the array x.
  implicit none
  real(kind=8), dimension(:), intent(in) :: x
  real(kind=8), dimension(size(x)), intent(out) :: f
  f = x**2
end subroutine fsub

end module sub3module

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

program sub3
    use sub3module
    implicit none
    real(kind=8), dimension(3) :: y,z

    y = (/2., 3., 4./)
    call fsub(y,z)
    print *, "z = ",z
end program sub3

max_stars_repo_name: philipwangdk/HPC

max_issues_repo_name: philipwangdk/HPC

max_forks_repo_name: philipwangdk/HPC


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine getdesiinfo(set,istartset,iendset,ialset,nset,
     &  mi,nactdof,ndesi,nodedesi,ntie,tieset,itmp,nmpc,nodempc,
     &  ipompc)
!
!     reading the input deck: *DESIGNVARIABLES
!
      implicit none
!
      character*81 setname
      character*81 set(*)
      character*81 tieset(3,*)
!
      integer mi(*),istartset(*),iendset(*),ialset(*),ndesi,
     &  node,nodedesi(*),nset,ntie,i,j,k,l,nmpc,nodempc(3,*),
     &  nactdof(0:mi(2),*),itmp(*),ntmp,index,id,ipompc(*)
!
      setname(1:1)=' '
      ndesi=0
!
!     Search for the set name
!
      do i=1,ntie
         if(tieset(1,i)(81:81).eq.'D') then
            setname=tieset(2,i)
         endif
      enddo
!
!     Check for the existence of the name
!
      if(setname(1:1).eq.' ') then
        write(*,*) '*ERROR in getdesiinfo: name of node set '
        write(*,*) '  has not yet been defined. '
        call exit(201)
      endif
!
!     catalogue all nodes (dependent and independent) which
!     belong to MPC's
!
      ntmp=0
      do i=1,nmpc
         index=ipompc(i)
         do
            if(index.eq.0) exit
            node=nodempc(1,index)
            call nident(itmp,node,ntmp,id)
            if(id.gt.0) then
               if(itmp(id).eq.node) then
                  index=nodempc(3,index)
                  cycle
               endif
            endif
            ntmp=ntmp+1
            do j=ntmp,id+2,-1
               itmp(j)=itmp(j-1)
            enddo
            itmp(id+1)=node
            index=nodempc(3,index)
         enddo
      enddo
!
!     Search the name of the node set in "set(i)" and
!     assign the nodes of the set to the appropriate variables
!
      do i=1,nset
         if(setname.eq.set(i)) then
            loop1: do j=istartset(i),iendset(i)
               if(ialset(j).gt.0) then
                  node=ialset(j)
!
!                 check for SPC-constraints: if a node is constrained in all dofs
!                 it is removed from the design node set
!
                  do l=1,3
c                     if(nactdof(l,node).le.0) cycle loop1
                     if(nactdof(l,node).gt.0) exit
                     if(l.eq.3) then
                        write(*,*) '*WARNING in getdesiinfo:'
                        write(*,*) '         node ',node,' has no'
                        write(*,*) '         active dofs and'
                        write(*,*) '         is removed from the set'
                        write(*,*) '         of design variables'
                        cycle loop1
                     endif
                  enddo
!
!                 check for MPC-constraints
!
                  call nident(itmp,node,ntmp,id)
                  if(id.gt.0) then
                     if(itmp(id).eq.node) then
                        write(*,*) '*WARNING in getdesiinfo:'
                        write(*,*) '         node ',node,' is subject'
                        write(*,*) '         to MPC-constraints and'
                        write(*,*) '         is removed from the set'
                        write(*,*) '         of design variables'
                        cycle loop1
                     endif
                  endif
!
                  ndesi=ndesi+1
                  nodedesi(ndesi)=node
               else
                  k=ialset(j-2)
                  loop2: do
                     k=k-ialset(j)
                     if(k.ge.ialset(j-1)) exit
!
!                 check for SPC-constraints: if a node is constrained in all dofs
!                 it is removed from the design node set
!
                     do l=1,3
                        if(nactdof(l,k).gt.0) exit
                        if(l.eq.3) then
                           write(*,*) '*WARNING in getdesiinfo:'
                           write(*,*) '         node ',k,' has no'
                           write(*,*) '         active dofs and'
                           write(*,*) '         is removed from the set'
                           write(*,*) '         of design variables'
                           cycle loop2
                        endif
                     enddo
!
!                 check for MPC-constraints
!
                     call nident(itmp,k,ntmp,id)
                     if(id.gt.0) then
                        if(itmp(id).eq.k) then
                           write(*,*) '*WARNING in getdesiinfo:'
                           write(*,*) '         node ',k,' is subject'
                           write(*,*) '         to MPC-constraints and'
                           write(*,*) '         is removed from the set'
                           write(*,*) '         of design variables'
                           cycle loop2
                        endif
                     endif
!
                     ndesi=ndesi+1
                     nodedesi(ndesi)=k
                  enddo loop2
               endif
            enddo loop1
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
c--- logical variable for determining whether or not to use complex pole
c--- scheme of Passarino et al. (HTO)
      logical CPscheme
      common/CPscheme/CPscheme

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


!=================================================
! SUBROUTINE
!	MNEWT(x1, y1, z1, x2, y2, z2, x3, y3, z3, di, da2)
!		intent(in) :: x1, y1, z1, x2, y2, z2, x3, di, da2
!		intent(out) :: y3, z3
!
!=================================================
module mod_mnewt

	implicit none

contains
!==========
subroutine MNEWT(x1,y1,z1, x2,y2,z2, x3,y3,z3, di,da2)

	use mod_intdep

	implicit none

	real(8), intent(in) :: x1, y1, z1, x2, y2, z2, x3, di, da2
	real(8), intent(out) :: y3, z3

	real(8) :: ct, cx, cy, cz, c2
	real(8) :: dx, dl, dr, dz, du
	real(8) :: arr, arl, ari

	ct = 1.d-6
	cx = x2 - x1
	cy = y2 - y1
	cz = z2 - z1
	c2 = cx*cx + cy*cy + cz*cz

	dx = x3 - x1
	dl = di
	dr = 0.d0
	arr = 0.d0

	do
		y3 = y1 + dl
		CALL INTDEP(x3,y3,z3)

		dz = z3 - z1
		CALL AREA2(cx, cy, cz, dx, dl, dz, c2, arl)

		if(dabs(arl-da2) < ct) return
		if(arl < da2) then
			dr = dl
			arr = arl
			dl = dl * da2 / arl
			cycle
		end if

		do
			du = dr + (dl - dr) * (da2 - arr) / (arl - arr)
			y3 = y1 + du
			CALL INTDEP(x3,y3,z3)

			dz = z3 - z1
			CALL AREA2(cx, cy, cz, dx, du, dz, c2, ari)

			if(dabs(ari-da2) < ct) return

			if(ari > da2) then
				arl = ari
				dl = du
			else if(ari < da2) then
				arr = ari
				dr = du
			end if

		end do

	end do

end subroutine MNEWT

end module mod_mnewt


text:
subroutine griddims( gdp )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: griddims.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/data/src/general/griddims.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Initialises grid related dimensions
!              This routine is introduced with the implementation
!              of DD boundaries
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    use dfparall
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer      , pointer :: nmax
    integer      , pointer :: mmax
    integer      , pointer :: nlb
    integer      , pointer :: nub
    integer      , pointer :: mlb
    integer      , pointer :: mub
    integer      , pointer :: nmlb
    integer      , pointer :: nmub
    integer      , pointer :: ddbound
    integer      , pointer :: jstart
    integer      , pointer :: nmmaxj
    integer      , pointer :: nmmax
    integer      , pointer :: numdomains
    integer      , pointer :: nummappers
!
! Global variables
!
!   NONE
!
! Local variables
!
    integer                :: istat
!
!! executable statements -------------------------------------------------------
!
    nmax         => gdp%d%nmax
    mmax         => gdp%d%mmax
    nlb          => gdp%d%nlb
    nub          => gdp%d%nub
    mlb          => gdp%d%mlb
    mub          => gdp%d%mub
    nmlb         => gdp%d%nmlb
    nmub         => gdp%d%nmub
    ddbound      => gdp%d%ddbound
    jstart       => gdp%d%jstart
    nmmaxj       => gdp%d%nmmaxj
    nmmax        => gdp%d%nmmax
    numdomains   => gdp%gdprognm%numdomains
    nummappers   => gdp%gdprognm%nummappers
    !
    ! basic dimensions:
    ! nmmax  = mmax   *   nmax
    ! nmmaxj = nmmax  + 2*nmax
    ! jstart = 1      - 2*nmax
    !
    !
    ! Probably, ddbound can be 0 when numdomains=1, nummappers>=1
    ! For safety, ddbound is always 1 when nummappers>=1
    !
    ! ddbound is expected to be able to be 0 when parll
    ! Unfortunately, it must be 1
    !
    if (nummappers>=1 .or. parll) then
       gdp%d%ddbound = 1
    else
       gdp%d%ddbound = 0
    endif
    gdp%d%nlb  = 1 - gdp%d%ddbound
    gdp%d%nub  = nmax + gdp%d%ddbound
    gdp%d%mlb  = -1 - gdp%d%ddbound
    gdp%d%mub  = mmax + 2 + gdp%d%ddbound
    nmmax      = (mmax + 2*gdp%d%ddbound)*(nmax + 2*gdp%d%ddbound)
    jstart     = 1 - 2*(nmax + 2*gdp%d%ddbound)
    nmmaxj     = nmmax + 2*(nmax + 2*gdp%d%ddbound)
    gdp%d%nmlb = jstart
    gdp%d%nmub = nmmaxj
    !
    gdp%griddim%mlb    = gdp%d%mlb
    gdp%griddim%mub    = gdp%d%mub
    gdp%griddim%mmax   = mmax
    !
    gdp%griddim%nlb    = gdp%d%nlb
    gdp%griddim%nub    = gdp%d%nub
    gdp%griddim%nmax   = nmax
    !
    gdp%griddim%nmlb   = gdp%d%nmlb
    gdp%griddim%nmub   = gdp%d%nmub
    gdp%griddim%nmmax  = nmmax
    !
    gdp%griddim%mfg    = gdp%gdparall%mfg
    gdp%griddim%mlg    = gdp%gdparall%mlg
    gdp%griddim%mmaxgl = gdp%gdparall%mmaxgl
    !
    gdp%griddim%nfg    = gdp%gdparall%nfg
    gdp%griddim%nlg    = gdp%gdparall%nlg
    gdp%griddim%nmaxgl = gdp%gdparall%nmaxgl
    !
    gdp%griddim%aggrtable => null()
    !
    allocate(gdp%griddim%celltype(gdp%d%nmlb:gdp%d%nmub), stat=istat)
    gdp%griddim%celltype(:) = 1
    !
    allocate(gdp%griddim%nmbnd(0,2), stat=istat)
end subroutine griddims


text:
!====================================================================
!-----------------------------------------------------------------------
!       NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-----------------------------------------------------------------------
!BOI
! !MODULE: m_temperature -- Temperature conversion routines
!
! !DESCRIPTION:
!     This module contains routines that convert temperature values.
!
! !REFERENCES:
!     Rogers, R. R., 1979: A Short Course in Cloud Physics, 2nd. Ed.,
!          Pergamon Press, Oxford. 235pp.
!
! !INTERFACE:
!
      module      m_temperature
      implicit    NONE
      private	! except

      public ::           ! Routines to convert from air temp (AT) to ...
     .   ATtoPT, AT2PT,   ! ... potential temperature (PT)
     .   ATtoVT, AT2VT    ! ... virtual temperature (VT)

      public ::           ! Routines to convert from ...
     .   PTtoAT, PT2AT,   ! ... potential temperature (PT)
     .   VTtoAT, VT2AT    ! ... virtual temperature (VT)
      public ::           ! Equation code for determining the VT
     .   VT_Exact,        ! ... using the exact equations
     .   VT_NCEP          ! ... NCEPS's approximation
!
      interface ATtoVT
         module procedure ! Convert from air temperature and ..
     .      ATtoVT_vp,    ! ... vapor pressure or
     .      ATtoVT_mr     ! ... mixing ratio to virtual temperature
      end interface

      interface AT2VT
         module procedure ! Convert from air temperature and ..
     .      AT2VT_vp,     ! ... vapor pressure or
     .      AT2VT_mr      ! ... mixing ratio to virtual temperature
      end interface
!
      interface VTtoAT
         module procedure ! Convert from virtual temperature and ..
     .      VTtoAT_vp,    ! ... vapor pressure or
     .      VTtoAT_mr     ! ... mixing ratio to air temperature
      end interface
!
      interface VT2AT
         module procedure ! Convert from virtual temperature and ..
     .      VT2AT_vp,     ! ... vapor pressure or
     .      VT2AT_mr      ! ... mixing ratio to air temperature
      end interface
!
!     Note: Routine names of the form xxtoxx and xx2xx implies a
!           subroutine and function, respectively.
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder  Original code
!EOI
!-----------------------------------------------------------------
      real, parameter ::
     .   Ra     =  287.0,   ! The individual gas constant for dry air and ...
     .   Rv     =  461.0,   ! ... water vapor (j/(kg-deg K, from Rogers,
     .   RaRv   =  Ra / Rv  !   1979)
      real, parameter ::
     .   Cp     =  1004.0,  ! Specific heat at constant pressure
     .   Cv     =   717.0,  ! ... at constant volume for dry air
     .                      ! (j/(kg-deg K, from Rogers, 1979)
     .   k      =  ( Cp - Cv ) / Cp,
     .   PS_Def =  1000.0   ! Default value for the final pressure for
                            !   the potential temperature (mb).

      integer, parameter ::
     .   VT_Exact =  0,       ! Exact equation
     .   VT_NCEP  =  1        ! ... NCEP's approximation

      contains
!......................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: ATtoPT () --- Air temperature to potential temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from air temperature
!     to potential temperature
!
! !INTERFACE:
      pure subroutine ATtoPT ( AT, P, PT, PS )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,      intent (in)  ::
     .   AT (:), ! air temperature (deg K)
     .   P  (:)  ! air pressure (mb or hPa)
      real,      intent (in), optional ::
     .   PS      ! reference pressure (mb or hPa).  Default: PS = 1000.0 mb
! !OUTPUT PARAMETERS:
      real,      intent (out) ::
     .   PT (:)  ! potential temperature (deg K)
!
!     Note: The minimum size of the arrays AT and PT determines the number
!           of values to convert.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!
! EOP
!-------------------------------------------------------------------------

      real    ::  T, PP, PSS
      integer :: iT, NT, NP

      PSS  = PS_Def
      if ( present ( PS )) PSS = PS

      NP   = size ( P )
      NT   = min  ( size ( AT ), size ( PT ))
      if ( NP .eq. 0 ) NT = 0
      do iT = 1, NT
         T  = AT ( iT )
         PP = PSS / P  ( min ( iT, NP ))
         PT ( iT ) = T * PP ** k
      end do

      return
      end subroutine ATtoPT

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

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: PTtoAT () --- Potential temperature to air temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from potential
!     temperature to air temperature
!
! !INTERFACE:
      pure subroutine PTtoAT ( PT, P, AT, PS )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,      intent (in)  ::
     .   PT (:), ! potential temperature (deg K)
     .   P  (:)  ! air pressure (mb or hPa)
      real,      intent (in), optional ::
     .   PS      ! reference pressure (mb or hPa).  Default: PS = 1000.0 mb
! !OUTPUT PARAMETERS:
      real,      intent (out) ::
     .   AT (:)  ! air temperature (deg K)
!
!     Note: The minimum size of the arrays PT and AT determines the number
!           of values to convert.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!
! EOP
!-------------------------------------------------------------------------

      real    ::  T, PP, PSS
      integer :: iT, NT, NP

      PSS  = PS_Def
      if ( present ( PS )) PSS = PS

      NP   = size ( P )
      NT   = min  ( size ( AT ), size ( PT ))
      if ( NP .eq. 0 ) NT = 0
      do iT = 1, NT
         T  = PT ( iT )
         PP = PSS / P  ( min ( iT, NP ))
         AT ( iT ) = T * PP ** (-k)
      end do

      return
      end subroutine PTtoAT

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

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: AT2PT () --- Air temperature to potential temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from air
!     temperature to potential temperature
!
! !INTERFACE:
      pure function AT2PT ( AT, P, PS )
      implicit NONE
!
! !INPUT PARAMETERS:
      real, intent (in)           :: AT    ! air temperature (deg K)
      real, intent (in)           :: P     ! air pressure (mb or hPa)
      real, intent (in), optional :: PS    ! reference pressure (mb or hPa).
                                           !   Default: PS = 1000.0 mb
! !OUTPUT PARAMETERS:
      real                        :: AT2PT ! potential temperature (deg K)
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!
! EOP
!-------------------------------------------------------------------------

      real :: PT (1)

      call ATtoPT ( (/AT/), (/P/), PT, PS )
      AT2PT = PT (1)

      return
      end function AT2PT

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

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: PT2AT () --- Potential temperature to air temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from potential
!     temperature to air temperature
!
! !INTERFACE:
      pure function PT2AT ( PT, P, PS )
      implicit NONE
!
! !INPUT PARAMETERS:
      real, intent (in)           :: PT    ! potential temperature (deg K)
      real, intent (in)           :: P     ! air pressure (mb or hPa)
      real, intent (in), optional :: PS    ! reference pressure (mb or hPa).
                                           !   Default: PS = 1000.0 mb
! !OUTPUT PARAMETERS:
      real                        :: PT2AT ! air temperature (deg K)
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!
! EOP
!-------------------------------------------------------------------------

      real :: AT (1)

      call PTtoAT ( (/PT/), (/P/), AT, PS )
      PT2AT = AT (1)

      return
      end function PT2AT
!......................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: ATtoVT_vp () --- Air temperature and vapor pressure to virtual temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from air temperature
!     and vapor pressure to virtual temperature.
!
! !INTERFACE:
      pure subroutine ATtoVT_vp ( AT, P, VP, VT, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,      intent (in)  ::
     .   AT (:), ! air temperature (deg K)
     .   P  (:), ! air pressure   (mb or hPa)
     .   VP (:)  ! vapor pressure (mb or hPa)
      integer,   intent (in), optional ::
     .   Eqtn    ! equation code (default: 0)
!
! !OUTPUT PARAMETERS:
      real,      intent (out) ::
     .   VT (:)  ! virtual temperature (deg K)
!
!     Note: The minimum size of the arrays AT, VP and VT determines the
!           number of values to convert.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real    ::  T, W,  PP
      integer :: iT, NT, NP, Code_

      Code_ = VT_Exact
      if ( present ( Eqtn )) Code_ = Eqtn

      NP   = size ( P )
      NT   = min  ( size ( AT ), size ( VP ), size ( VT ))
      if ( NP .eq. 0 ) NT = 0

      if ( Code_ .eq. VT_NCEP ) then
         do iT = 1, NT
            T  = AT ( iT )
            PP = P  ( min ( iT, NP ))
            W  = RaRv * VP ( iT ) / ( PP - VP ( iT ))
            VT ( iT ) = T * ( 1.0 + 0.61 * W )
         end do

      else
         do iT = 1, NT
            T  = AT ( iT )
            PP = P  ( min ( iT, NP ))
            W  = RaRv * VP ( iT ) / ( PP - VP ( iT ))
            VT ( iT ) = T * ( 1.0 + W / RaRv )
     .                    / ( 1.0 + W )
         end do

      end if

      return
      end subroutine ATtoVT_vp
!......................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: ATtoVT_mr () --- Air temperature and mixing ratio to virtual temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from air temperature
!     and mixing ratio to virtual temperature
!
! !INTERFACE:
      pure subroutine ATtoVT_mr ( AT, MR, VT, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,      intent (in)  ::
     .   AT (:), ! air temperature (deg K)
     .   MR (:)  ! mixing ratio (g/kg)
      integer,   intent (in), optional ::
     .   Eqtn    ! equation code (default: 0)
! !OUTPUT PARAMETERS:
      real,      intent (out) ::
     .   VT (:)  ! virtual temperature (deg K)
!
!     Note: The minimum size of the arrays AT, MR and VT determines the
!           number of values to convert.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real    ::  T, MRR
      integer :: iT, NT, Code_

      Code_ = VT_Exact
      if ( present ( Eqtn )) Code_ = Eqtn

      NT = min  ( size ( AT ), size ( MR ), size ( VT ))

      if ( Code_ .eq. VT_NCEP ) then
         do iT = 1, NT
            T   = AT ( iT )
            MRR = 0.001 * MR ( iT )
            VT ( iT ) = T * ( 1.0 + 0.61 * MRR )
         end do

      else
         do iT = 1, NT
            T   = AT ( iT )
            MRR = 0.001 * MR ( iT )
            VT ( iT ) = T * ( 1.0 + MRR / RaRv )
     .                    / ( 1.0 + MRR )
         end do

      end if

      return
      end subroutine ATtoVT_mr
!......................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: VTtoAT_vp () --- Virtual temperature and vapor pressure to airtemperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from virtual
!     temperature and vapor pressure to air temperature.
!
! !INTERFACE:
      pure subroutine VTtoAT_vp ( VT, P, VP, AT, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,      intent (in)  ::
     .   VT (:), ! virtual temperature (deg K)
     .   P  (:), ! air pressure   (mb or hPa)
     .   VP (:)  ! vapor pressure (mb or hPa)
      integer,   intent (in), optional ::
     .   Eqtn    ! equation code (default: 0)
! !OUTPUT PARAMETERS:
      real,      intent (out) ::
     .   AT (:)  ! air temperature (deg K)
!
!     Note: The minimum size of the arrays VT, VP and AT determines the
!           number of values to convert.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real    ::  T, W,  PP
      integer :: iT, NT, NP, Code_

      Code_ = VT_Exact
      if ( present ( Eqtn )) Code_ = Eqtn

      NP   = size ( P )
      NT   = min  ( size ( VT ), size ( VP ), size ( AT ))
      if ( NP .eq. 0 ) NT = 0

      if ( Code_ .eq. VT_NCEP ) then
         do iT = 1, NT
            T  = VT ( iT )
            PP = P  ( min ( iT, NP ))
            W  = RaRv * VP ( iT ) / ( PP - VP ( iT ))
            AT ( iT ) = T / ( 1.0 + 0.61 * W )

         end do

      else
         do iT = 1, NT
            T  = VT ( iT )
            PP = P  ( min ( iT, NP ))
            W  = RaRv * VP ( iT ) / ( PP - VP ( iT ))
            AT ( iT ) = T * ( 1.0 + W )
     .                    / ( 1.0 + W / RaRv )
         end do

      end if

      return
      end subroutine VTtoAT_vp
!......................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: VTtoAT_mr () --- Virtual temperature and mixing ratio to air temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from virtural
!     temperature and mixing ratio to air temperature
!
! !INTERFACE:
      pure subroutine VTtoAT_mr ( VT, MR, AT, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,      intent (in)  ::
     .   VT (:), ! virtual temperature (deg K)
     .   MR (:)  ! mixing ratio (g/kg)
      integer,   intent (in), optional ::
     .   Eqtn    ! equation code (default: 0)
! !OUTPUT PARAMETERS:
      real,      intent (out) ::
     .   AT (:)  ! air temperature (deg K)
!
!     Note: The minimum size of the arrays AT, MR and VT determines the
!           number of values to convert.
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real    ::  T, MRR
      integer :: iT, NT, Code_

      Code_ = VT_Exact
      if ( present ( Eqtn )) Code_ = Eqtn

      NT   = min  ( size ( VT ), size ( MR ), size ( AT ))

      if ( Code_ .eq. VT_NCEP ) then
         do iT = 1, NT
            T   = VT ( iT )
            MRR = 0.001 * MR ( iT )
            AT ( iT ) = T / ( 1.0 + 0.61 * MRR )
         end do
      else
         do iT = 1, NT
            T   = VT ( iT )
            MRR = 0.001 * MR ( iT )
            AT ( iT ) = T * ( 1.0 + MRR )
     .                    / ( 1.0 + MRR / RaRv )
         end do

      end if

      return
      end subroutine VTtoAT_mr
!......................................................................

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: AT2VT_vp () --- Air temperature and vapor pressure to virtual temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from air
!     temperature and vapor pressure to virtual temperature
!
! !INTERFACE:
      pure function AT2VT_vp ( AT, P, VP, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,    intent (in) :: AT       ! air temperature (deg K)
      real,    intent (in) :: P        ! air pressure (mb or hPa)
      real,    intent (in) :: VP       ! vapor pressure (mb or hPa)
      integer, intent (in),   optional ::
     .                        Eqtn     ! equation code (default: 0)
!
! !OUTPUT PARAMETERS:
      real                 :: AT2VT_vp ! virtual temperature (deg K)
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real :: VT (1)

      call ATtoVT_vp ( (/AT/), (/P/), (/VP/), VT, Eqtn )
      AT2VT_vp = VT (1)

      return
      end function AT2VT_vp

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

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: AT2VT_mr () --- Air temperature and mixing ratio to virtual temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from air
!     temperature and mixing ratio to virtual temperature
!
! !INTERFACE:
      pure function AT2VT_mr ( AT, MR, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,    intent (in) :: AT       ! air temperature (deg K)
      real,    intent (in) :: MR       ! vapor pressure  (mb or hPa)
      integer, intent (in),   optional ::
     .                        Eqtn     ! equation code (default: 0)
!
! !OUTPUT PARAMETERS:
      real                 :: AT2VT_mr ! virtual temperature (deg K)
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real :: VT (1)

      call ATtoVT_mr ( (/AT/), (/MR/), VT, Eqtn )
      AT2VT_mr = VT (1)

      return
      end function AT2VT_mr

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

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: VT2AT_vp () --- Virtual temperature and vapor pressure to air temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from virtual
!     temperature and vapor pressure to air temperature
!
! !INTERFACE:
      pure function VT2AT_vp ( VT, P, VP, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,    intent (in) :: VT       ! virtual temperature (deg K)
      real,    intent (in) :: P        ! air pressure (mb or hPa)
      real,    intent (in) :: VP       ! vapor pressure (mb or hPa)
      integer, intent (in),   optional ::
     .                        Eqtn     ! equation code (default: 0)
!
! !OUTPUT PARAMETERS:
      real                 :: VT2AT_vp ! air temperature (deg K)
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real :: AT (1)

      call VTtoAT_vp ( (/VT/), (/P/), (/VP/), AT, Eqtn )
      VT2AT_vp = AT (1)

      return
      end function VT2AT_vp

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

!-------------------------------------------------------------------------
!         Nasa/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
! BOP
!
! !ROUTINE: VT2AT_mr () --- Virtual temperature and mixing ratio to air temperature
!
! !DESCRIPTION:
!     This routine converts a set of temperature values from virtual
!     temperature and mixing ratio to air temperature
!
! !INTERFACE:
      pure function VT2AT_mr ( VT, MR, Eqtn )
      implicit NONE
!
! !INPUT PARAMETERS:
      real,    intent (in) :: VT       ! virtual temperature (deg K)
      real,    intent (in) :: MR       ! vapor pressure  (mb or hPa)
      integer, intent (in),   optional ::
     .                        Eqtn     ! equation code (default: 0)
!
! !OUTPUT PARAMETERS:
      real                 :: VT2AT_mr ! air temperature (deg K)
!
! !SEE ALSO:
!
! !REVISION HISTORY:
!     23Jul2002  C. Redder   Original code.
!     20Apr2004  C. Redder   Added the optional argument, Eqtn.
!
! EOP
!-------------------------------------------------------------------------

      real :: AT (1)

      call VTtoAT_mr ( (/VT/), (/MR/), AT, Eqtn )
      VT2AT_mr = AT (1)

      return
      end function VT2AT_mr
!......................................................................
      end module m_temperature
!====================================================================

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


text:
      double complex function qlLi2omx2(v,w,x,y)
      implicit none
      include 'qlconstants.f'
c     expression for dilog(1-(v-i*ep)*(w-i*ep)/(x-i*ep)/(y-i*ep))
C     for real v,w,x and y
      double precision v,w,x,y,omarg,arg,ddilog
      double complex qllnrat,lnarg,lnomarg,prod
      arg=(v*w)/(x*y)
      omarg=1d0-arg
      if (arg .le. 1d0) then
         if (arg .eq. 0d0 .or. arg .eq. 1d0) then
            prod=0d0
         else
            lnarg=qllnrat(v,x)+qllnrat(w,y)
            lnomarg=dcmplx(log(omarg))
            prod=lnarg*lnomarg
         endif
         qlLi2omx2=dcmplx(pisqo6-ddilog(arg))-prod
      elseif (arg .gt. 1d0) then
         arg=(x*y)/(v*w)
         lnarg=-qllnrat(v,x)-qllnrat(w,y)
         lnomarg=dcmplx(log(1d0-arg))
         qlLi2omx2=-dcmplx(pisqo6-ddilog(arg))+lnarg*lnomarg
     .   -0.5d0*lnarg**2
      endif
      return
      end


! We need to test how many test cases are allowed.
module passUnderCaseLimit_M

#include "t123/TestFile.f90i"

implicit none

contains

TEST( Case, 1 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 2 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 3 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 4 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 5 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 6 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 7 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 8 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 9 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 10 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 11 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 12 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 13 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 14 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 15 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 16 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 17 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 18 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 19 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 20 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 21 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 22 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 23 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 24 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 25 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 26 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 27 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 28 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 29 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 30 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 31 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 32 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 33 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 34 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 35 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 36 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 37 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 38 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 39 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 40 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 41 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 42 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 43 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 44 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 45 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 46 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 47 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 48 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 49 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 50 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 51 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 52 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 53 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 54 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 55 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 56 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 57 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 58 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 59 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 60 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 61 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 62 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 63 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 64 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 65 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 66 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 67 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 68 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 69 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 70 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 71 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 72 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 73 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 74 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 75 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 76 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 77 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 78 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 79 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 80 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 81 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 82 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 83 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 84 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 85 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 86 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 87 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 88 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 89 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 90 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 91 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 92 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 93 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 94 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 95 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 96 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 97 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 98 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 99 )
    EXPECT_EQ(1,1)
END TEST
TEST( Case, 100 )
    EXPECT_EQ(1,1)
END TEST

end module

program main
use passUnderCaseLimit_M
!
RUN_ALL_TESTS()
!
end program


text:
c
c
c
c     ===============================================================
      subroutine claw1ez(maxmx,meqn,mwaves,mbc,maux,mwork,mthlim,
     &                   q,work,aux)
c     ===============================================================
c
c     An easy-to-use clawpack driver routine for simple applications
c     Documentation is available at
c                 http://www.amath.washington.edu/~claw/doc.html
c
c     Author: Randall J. LeVeque
c     Version of August, 2002 --  CLAWPACK Version 4.1
c
      implicit double precision (a-h,o-z)
      external bc1,rp1,src1,b4step1

      dimension    q(1-mbc:maxmx+mbc, meqn)
      dimension  aux(1-mbc:maxmx+mbc, maux)
      dimension work(mwork)
      dimension mthlim(mwaves)
c
      dimension method(7),dtv(5),cflv(4),nv(2),mthbc(2)
      dimension tout(100)
      logical outt0
      character*12 fname
c
      open(10,file='fort.info',status='unknown',form='formatted')

c
c     # Read the input in standard form from claw.data:
c     # For a description of input parameters see the documentation at
c                 http://www.amath.washington.edu/~claw

c     ## New in 4.4:   Input file name changed from claw1ez.data
c     ## Open file and skip over leading lines with # comments:
      fname = 'claw.data'
      call opendatafile(55,fname)
c
c     Number of space dimensions:  ## New in 4.4
      read(55,*) ndim
c
c     ## The remainder is unchanged from 4.3:

c     Number of grid cells:
      read(55,*) mx

c     i/o variables
      read(55,*) nout
      read(55,*) outstyle
      if (outstyle.eq.1) then
          read(55,*) tfinal
          nstepout = 1
        elseif (outstyle.eq.2) then
          read(55,*) (tout(i), i=1,nout)
          nstepout = 1
        elseif (outstyle.eq.3) then
          read(55,*) nstepout, nstop
	  nout = nstop
        endif

c     timestepping variables
      read(55,*) dtv(1)
      read(55,*) dtv(2)
      read(55,*) cflv(1)
      read(55,*) cflv(2)
      read(55,*) nv(1)
c

c     # input parameters for clawpack routines
      read(55,*) method(1)
      read(55,*) method(2)
      read(55,*) method(3)
      read(55,*) method(4)
      read(55,*) method(5)
      read(55,*) method(6)
      read(55,*) method(7)

      read(55,*) meqn1
      read(55,*) mwaves1
      read(55,*) (mthlim(mw), mw=1,mwaves1)

c     # physical domain:
      read(55,*) t0
      read(55,*) xlower
      read(55,*) xupper
c
c     # boundary conditions:
      read(55,*) mbc1
      read(55,*) mthbc(1)
      read(55,*) mthbc(2)

      if (method(7) .ne. maux) then
         write(6,*) '*** ERROR ***  maux set wrong in input or driver'
         stop
         endif

      if (meqn1 .ne. meqn) then
         write(6,*) '*** ERROR ***  meqn set wrong in input or driver'
         stop
         endif
      if (mwaves1 .ne. mwaves) then
         write(6,*) '*** ERROR ***  mwaves set wrong in input or driver'
         stop
         endif
      if (mbc1 .ne. mbc) then
         write(6,*) '*** ERROR ***  mbc set wrong in input or driver'
         stop
         endif
c
      if ((mthbc(1).eq.2 .and. mthbc(2).ne.2) .or.
     &    (mthbc(2).eq.2 .and. mthbc(1).ne.2)) then
         write(6,*) '*** ERROR ***  periodic boundary conditions'
         write(6,*) ' require mthbc(1) and mthbc(2) BOTH be set to 2'
         stop
         endif

c
c     # check that enough storage has been allocated:
c

      mwork1 = (maxmx + 2*mbc) * (2 + 4*meqn + mwaves + meqn*mwaves)
c
      if (mx.gt.maxmx .or. mwork.lt.mwork1) then
c        # insufficient storage
         maxmx1 = max0(mx,maxmx)

         mwork1 = (maxmx1 + 2*mbc) * (2 + 4*meqn + mwaves + meqn*mwaves)

         write(6,*) ' '
         write(6,*) '*** ERROR *** Insufficient storage allocated'
         write(6,*) 'Recompile after increasing values in driver.f:'
         write(6,611) maxmx1
         write(6,613) mwork1
 611     format(/,'parameter (maxmx = ',i5,')')
 613     format('parameter (mwork = ',i7,')',/)
         stop
         endif

c
c
      write(6,*) 'running...'
      write(6,*) ' '
c
c     # grid spacing
      dx = (xupper - xlower) / float(mx)
c

c     # time increments between outputing solution:
      if (outstyle .eq. 1) then
         dtout = (tfinal - t0)/float(nout)
         endif
c
c
c     # call user's routine setprob to set any specific parameters
c     # or other initialization required.
c
      call setprob
c
c     # set aux array:
c
      if (maux .gt. 0)  then
         call setaux(maxmx,mbc,mx,xlower,dx,maux,aux)
         endif
c
c     # set initial conditions:
c
      call qinit(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux)
c
      outt0 = .true.
      if (outt0) then
c        # output initial data
         call out1(maxmx,meqn,mbc,mx,xlower,dx,q,t0,0,aux,maux)
         write(6,601) 0, t0
         endif

c
c     ----------
c     Main loop:
c     ----------
c
      tend = t0
      do 100 n=1,nout
         tstart = tend
         if (outstyle .eq. 1)  tend = tstart + dtout
         if (outstyle .eq. 2)  tend = tout(n)
         if (outstyle .eq. 3)  tend = tstart - 1.d0  !# single-step mode
c
         call claw1(maxmx,meqn,mwaves,mbc,mx,
     &           q,aux,xlower,dx,tstart,tend,dtv,cflv,nv,method,mthlim,
     &           mthbc,work,mwork,info,bc1,rp1,src1,b4step1)
c
c        # check to see if an error occured:
         if (info .ne. 0) then
            write(6,*) '*** ERROR in claw1 ***  info =',info
            if (info.eq.1) then
               write(6,*) '***   either mx > maxmx or mbc < 2'
               endif
            if (info.eq.2) then
               write(6,*) '***   dt does not divide (tend - tstart)'
               write(6,*) '***   and dt is fixed since method(1)=0'
               endif
            if (info.eq.3) then
               write(6,*) '***   method(1)=1 and cflv(2) > cflv(1)'
               endif
            if (info.eq.4) then
               write(6,*) '***   mwork is too small'
               endif
            if (info.eq.11) then
               write(6,*) '***   Too many times steps, n > nv(1)'
               endif
            if (info.eq.12) then
               write(6,*)
     &          '***   The Courant number is greater than cflv(1)'
               write(6,*) '***   and dt is fixed since method(1)=0'
               endif

            go to 999
            endif
c
         dtv(1) = dtv(5)  !# use final dt as starting value on next call
c
c        # output solution at this time
c        ------------------------------
c
c        # if outstyle=1 or 2, then nstepout=1 and we output every time
c        # we reach this point, since claw1 was called for the entire time
c        # increment between outputs.
c
c        # if outstyle=3 then we only output if we have taken nstepout
c        # time steps since the last output.

c        # iframe is the frame number used to form file names in out1
         iframe = n/nstepout
	 if (iframe*nstepout .eq. n) then
            call out1(maxmx,meqn,mbc,mx,xlower,dx,q,tend,iframe,
     &                aux,maux)
            write(6,601) iframe,tend
            write(10,1010) tend,info,dtv(3),dtv(4),dtv(5),
     &           cflv(3),cflv(4),nv(2)
	    endif
c
c        # formats for writing out information about this call to claw:
c
  601    format('CLAW1EZ: Frame ',i4,
     &           ' output files done at time t =',
     &           d12.4,/)
c
 1010    format('tend =',d15.4,/,
     &       'info =',i5,/,'smallest dt =',d15.4,/,'largest dt =',
     &       d15.4,/,'last dt =',d15.4,/,'largest cfl =',
     &         d15.4,/,'last cfl =',d15.4,/,'steps taken =',i4,/)
c
  100    continue
c
  999 continue
c
      return
      end

max_issues_repo_name: che-wenchao/D-Claw

max_forks_repo_name: che-wenchao/D-Claw


text:
      INTEGER FUNCTION test()
        IMPLICIT NONE
        INCLUDE "acc_testsuite.fh"
        INTEGER :: x !Iterators
        REAL(8),DIMENSION(LOOPCOUNT):: a, b, c !Data
        INTEGER :: errors = 0

        !Initilization
        CALL RANDOM_SEED()
        CALL RANDOM_NUMBER(a)
        CALL RANDOM_NUMBER(b)
        c = 0

        !$acc enter data copyin(a(1:LOOPCOUNT), b(1:LOOPCOUNT), &
           !$acc c(1:LOOPCOUNT)) async(1)
        !$acc kernels wait(1)
          !$acc loop
          DO x = 1, LOOPCOUNT
            c(x) = a(x) + b(x)
          END DO
        !$acc end kernels
        !$acc exit data copyout(c(1:LOOPCOUNT)) &
           !$acc delete(a(1:LOOPCOUNT), b(1:LOOPCOUNT))

        DO x = 1, LOOPCOUNT
          IF (abs(c(x) - (a(x) + b(x))) .gt. PRECISION) THEN
            errors = errors + 1
          END IF
        END DO
        test = errors
      END

      PROGRAM test_kernels_async_main
      IMPLICIT NONE
      INTEGER :: failed, success !Number of failed/succeeded tests
      INTEGER :: num_tests,crosschecked, crossfailed, j
      INTEGER :: temp,temp1
      INCLUDE "acc_testsuite.fh"
      INTEGER test

      CHARACTER*50:: logfilename !Pointer to logfile
      INTEGER :: result

      num_tests = 0
      crosschecked = 0
      crossfailed = 0
      result = 1
      failed = 0

      !Open a new logfile or overwrite the existing one.
      logfilename = "OpenACC_testsuite.log"
!      WRITE (*,*) "Enter logFilename:"
!      READ  (*,*) logfilename

      OPEN (1, FILE = logfilename)

      WRITE (*,*) "######## OpenACC Validation Suite V 2.5 ######"
      WRITE (*,*) "## Repetitions:", N
      WRITE (*,*) "## Loop Count :", LOOPCOUNT
      WRITE (*,*) "##############################################"
      WRITE (*,*)

      WRITE (*,*) "--------------------------------------------------"
      WRITE (*,*) "Test of kernels_wait"
      WRITE (*,*) "--------------------------------------------------"

      crossfailed=0
      result=1
      WRITE (1,*) "--------------------------------------------------"
      WRITE (1,*) "Test of kernels_wait"
      WRITE (1,*) "--------------------------------------------------"
      WRITE (1,*)
      WRITE (1,*) "(Crosstests should fail)"
      WRITE (1,*)

      DO j = 1, N
        temp =  test()
        IF (temp .EQ. 0) THEN
          WRITE (1,*)  j, ". test successfull."
          success = success + 1
        ELSE
          WRITE (1,*) "Error: ",j, ". test failed."
          failed = failed + 1
        ENDIF
      END DO

      IF (failed .EQ. 0) THEN
        WRITE (1,*) "Directive worked without errors."
        WRITE (*,*) "Directive worked without errors."
        result = 0
        WRITE (*,*) "Result:",result
      ELSE
        WRITE (1,*) "Directive failed the test ", failed, " times."
        WRITE (*,*) "Directive failed the test ", failed, " times."
        result = failed * 100 / N
        WRITE (*,*) "Result:",result
      ENDIF
      CALL EXIT (result)
      END PROGRAM

max_stars_repo_name: clementval/OpenACCV-V

max_issues_repo_name: clementval/OpenACCV-V

max_forks_repo_name: clementval/OpenACCV-V


      subroutine sd_hydsed_init

      use input_file_module
      use sd_channel_module
      use channel_velocity_module
      use maximum_data_module
      use hydrograph_module

      implicit none

      real :: kh
      integer :: idb                  !             |
      integer :: i                    !none         |counter
      real :: aa                      !none         |area/area=1 (used to calculate velocity with
                                      !             |Manning"s equation)
      real :: a                       !m^2          |cross-sectional area of channel
      real :: b                       !m            |bottom width of channel
      real :: d                       !m            |depth of flow
      real :: p                       !m            |wetting perimeter
      real :: chside                  !none         |change in horizontal distance per unit
                                      !             |change in vertical distance on channel side
                                      !             |slopes; always set to 2 (slope=1/2)
      real :: fps                     !none         |change in horizontal distance per unit
                                      !             |change in vertical distance on floodplain side
                                      !             |slopes; always set to 4 (slope=1/4)
      integer :: max                  !             |
      real :: rh                      !m            |hydraulic radius
      real :: qman                    !m^3/s or m/s |flow rate or flow velocity
      real :: tt1                     !km s/m       |time coefficient for specified depth
      real :: tt2                     !km s/m       |time coefficient for bankfull depth
      real :: qq1                     !m^3/s        |flow rate for a specified depth

      do i = 1, sp_ob%chandeg
        icmd = sp_ob1%chandeg + i - 1
        idb = ob(icmd)%props

        sd_ch(i)%order = sd_chd(idb)%order
        sd_ch(i)%chw = sd_chd(idb)%chw
        sd_ch(i)%chd = sd_chd(idb)%chd
        sd_ch(i)%chs = sd_chd(idb)%chs
        sd_ch(i)%chl = sd_chd(idb)%chl
        sd_ch(i)%cherod = sd_chd(idb)%cherod
        sd_ch(i)%cov = sd_chd(idb)%cov
        sd_ch(i)%shear_bnk = sd_chd(idb)%shear_bnk
        sd_ch(i)%hc_len = sd_chd(idb)%hc_ini
        sd_ch(i)%hc_hgt = sd_chd(idb)%hc_hgt

        !! compute headcut parameters
        kh = sd_chd(idb)%hc_kh
        if (kh > 1.e-6) then
          sd_ch(i)%hc_co = .37 * (17.83 + 16.56 * kh - 15. * sd_chd(idb)%hc_cov)
          sd_ch(i)%hc_co = amax1 (0., sd_ch(i)%hc_co)
        else
          sd_ch(i)%hc_co = 0.
        end if

        !! compute travel time coefficients
        aa = 1.
        b = 0.
        d = 0.
        chside = sd_chd(idb)%chss
        fps = 4.
        b = sd_ch(i)%chw - 2. * sd_ch(i)%chd * chside

        !! check IF bottom width (b) is < 0
        if (b <= 0.) then
            b = .5 * sd_ch(i)%chw
            b = Max(0., b)
            chside = (sd_ch(i)%chw - b) / (2. * sd_ch(i)%chd)
        end if
        sd_ch_vel(i)%wid_btm = b
        sd_ch_vel(i)%dep_bf = sd_ch(i)%chd

        !! compute flow and travel time at bankfull depth
        p = b + 2. * sd_ch(i)%chd * Sqrt(chside * chside + 1.)
        a = b * sd_ch(i)%chd + chside * sd_ch(i)%chd * sd_ch(i)%chd
        rh = a / p
        sd_ch_vel(i)%area = a
        sd_ch_vel(i)%vel_bf = Qman(a, rh, sd_chd(idb)%chn, sd_ch(i)%chs)
        sd_ch_vel(i)%velav_bf = Qman(aa, rh, sd_chd(idb)%chn, sd_ch(i)%chs)
        sd_ch_vel(i)%celerity_bf = sd_ch_vel(i)%velav_bf * 5. / 3.
        sd_ch_vel(i)%st_dis = sd_ch(i)%chl / sd_ch_vel(i)%celerity_bf / 3.6
        tt2 = sd_ch(i)%chl * a / sd_ch_vel(i)%vel_bf

        !! compute flow and travel time at 1.2 bankfull depth
        d = 1.2 * sd_ch(i)%chd
        a = a + (sd_ch(i)%chw * sd_ch(i)%chd + fps * (sd_ch(i)%chd - sd_ch(i)%chd) ** 2)
        p = p + 4. * sd_ch(i)%chw + (0.4 * sd_ch(i)%chd * Sqrt(fps * fps + 1.))
        rh = a / p
        qq1 = Qman(a, rh, sd_chd(idb)%chn, sd_ch(i)%chs)
        tt1 = sd_ch(i)%chl * a / qq1

        !! compute flow and travel time at 0.1 bankfull depth
        d = 0.1 * sd_ch(i)%chd
        p = b + 2. * sd_ch(i)%chd * Sqrt(chside * chside + 1.)
        a = b * sd_ch(i)%chd + chside * sd_ch(i)%chd * sd_ch(i)%chd
        rh = a / p
        qq1 = Qman(a, rh, sd_chd(idb)%chn, sd_ch(i)%chs)
        tt1 = sd_ch(i)%chl * a / qq1
        sd_ch_vel(i)%vel_1bf = Qman(aa, rh, sd_chd(idb)%chn, sd_ch(i)%chs)
        sd_ch_vel(i)%celerity_1bf = sd_ch_vel(i)%vel_1bf * 5. / 3.
        sd_ch_vel(i)%stor_dis_1bf = sd_ch(i)%chl / sd_ch_vel(i)%celerity_1bf / 3.6
        end do

      return
      end subroutine sd_hydsed_init