text:
      subroutine putdat(flnm, artype,time,
     &                  lsteric,icegln,trcout,
     &                  iexpt,iversn,yrflag,kkout, thbase)
      use mod_mom6  ! HYCOM mom6 array interface
      use mod_za    ! HYCOM array I/O interface
c
      character        flnm*(*)
      double precision time(3)
      real             thbase
      integer          artype,iexpt,iversn,yrflag,kkout
      logical          lsteric,icegln,trcout
c
c --- write model fields.
c --- HYCOM 2.0 array I/O archive file.
c
      character*8 ctype
c
      real      coord,xmin,xmax
      integer   i,j,k,ktr,jversn,l,nop
      data nop/24/
c
c
c ---   output is in "*.[ab]"
c
      l = len_trim(flnm)
      if     (flnm(l-1:l).eq.'.a' .or. flnm(l-1:l).eq.'.b') then
        open (unit=nop,file=flnm(1:l-2)//'.b',form='formatted',
     &          status='new',action='write')
        call zaiopf(flnm(1:l-2)//'.a','new', nop)
      else
        open (unit=nop,file=flnm(1:l)//'.b',form='formatted',
     &          status='new',action='write')
        call zaiopf(flnm(1:l)//'.a','new', nop)
      endif
c
c --- header.
c
      jversn = max(iversn,20)
      if     (artype.eq.1) then
        write(nop,116) ctitle,jversn,iexpt,yrflag,idm,jdm
        write( lp,  *)
        write( lp,116) ctitle,jversn,iexpt,yrflag,idm,jdm
 116    format (a80/a80/a80/a80/
     &   i5,4x,'''iversn'' = hycom version number x10'/
     &   i5,4x,'''iexpt '' = experiment number x10'/
     &   i5,4x,'''yrflag'' = days in year flag'/
     &   i5,4x,'''idm   '' = longitudinal array size'/
     &   i5,4x,'''jdm   '' = latitudinal  array size'/
     &   'field       time step  model day',
     &   '  k  dens        min              max')
      elseif (artype.eq.2) then
        write(nop,118) ctitle,jversn,iexpt,yrflag,idm,jdm
        write( lp,  *)
        write( lp,118) ctitle,jversn,iexpt,yrflag,idm,jdm
 118    format (a80/a80/a80/a80/
     &   i5,4x,'''iversn'' = hycom version number x10'/
     &   i5,4x,'''iexpt '' = experiment number x10'/
     &   i5,4x,'''yrflag'' = days in year flag'/
     &   i5,4x,'''idm   '' = longitudinal array size'/
     &   i5,4x,'''jdm   '' = latitudinal  array size'/
     &   'field         no. recs  mean day',
     &   '  k  dens        min              max')
      else
        write( lp,"(/ a /)")
     &    'error in putdat - only artpe==1 and artype==2 allowed'
        stop
      endif

c
c --- surface fields.
c
      coord=0.0
c
      call zaiowr(montg,ip,.true.,
     &            xmin,xmax, nop, .false.)
      if     (sigver.gt.0) then
c ---   identify the equation of state on the first record
        write (nop,117) 'montg1  ',nstep,time(1),sigver,thbase,xmin,xmax
        call flush(nop)
        write ( lp,117) 'montg1  ',nstep,time(1),sigver,thbase,xmin,xmax
        call flush( lp)
      else
        write (nop,117) 'montg1  ',nstep,time(1),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'montg1  ',nstep,time(1),0,coord,xmin,xmax
        call flush( lp)
      endif !sigver
      call zaiowr(srfht,ip,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'srfhgt  ',nstep,time(2),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'srfhgt  ',nstep,time(2),0,coord,xmin,xmax
      call flush( lp)
c
      if(lsteric) then
        call zaiowr(steric,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'steric  ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'steric  ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
      endif
c
      call zaiowr(surflx,ip,.true., xmin,xmax, nop, .false.)
      write (nop,117) 'surflx  ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'surflx  ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
      call zaiowr(salflx,ip,.true., xmin,xmax, nop, .false.)
      write (nop,117) 'salflx  ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'salflx  ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
c
      call zaiowr(dpbl,ip,.true., xmin,xmax, nop, .false.)
      write (nop,117) 'bl_dpth ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'bl_dpth ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
      call zaiowr(dpmixl,ip,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'mix_dpth',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'mix_dpth',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
      if     (sigver.eq.0) then
        call zaiowr(tmix,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'tmix    ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'tmix    ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
        call zaiowr(smix,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'smix    ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'smix    ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
        call zaiowr(thmix,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'thmix   ',nstep,time(3),0,thbase,xmin,xmax
        call flush(nop)
        write ( lp,117) 'thmix   ',nstep,time(3),0,thbase,xmin,xmax
        call flush( lp)
        call zaiowr(umix,iu,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'umix    ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'umix    ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
        call zaiowr(vmix,iv,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'vmix    ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'vmix    ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
      endif !sigver==0
      if(artype.gt.1) then
        call zaiowr(kemix, ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'kemix   ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'kemix   ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
      endif
      if(icegln) then
        call zaiowr(covice,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'covice  ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'covice  ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
        call zaiowr(thkice,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'thkice  ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'thkice  ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
        call zaiowr(temice,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'temice  ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'temice  ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
      endif
c
c --- depth averaged fields (no mask, in case this is a subregion archive).
c
      call zaiowr(ubaro,iu,.false.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
      call zaiowr(vbaro,iv,.false.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
      if(artype.gt.1) then
        call zaiowr(kebaro,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'kebtrop ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'kebtrop ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
      endif
c
c --- layer loop.
c
      do 75 k=1,kkout
      coord=theta(k)
      call zaiowr(u(1,1,k),iu,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'u-vel.  ',nstep,time(3),k,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'u-vel.  ',nstep,time(3),k,coord,xmin,xmax
      call flush( lp)
      call zaiowr(v(1,1,k),iv,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'v-vel.  ',nstep,time(3),k,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'v-vel.  ',nstep,time(3),k,coord,xmin,xmax
      call flush( lp)
      if(artype.gt.1) then
        call zaiowr(ke(1,1,k),ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'k.e.    ',nstep,time(3),k,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'k.e.    ',nstep,time(3),k,coord,xmin,xmax
        call flush( lp)
      endif
      call zaiowr(dp(1,1,k),  ip,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'thknss  ',nstep,time(3),k,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'thknss  ',nstep,time(3),k,coord,xmin,xmax
      call flush( lp)
      call zaiowr(temp(1,1,k),ip,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'temp    ',nstep,time(3),k,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'temp    ',nstep,time(3),k,coord,xmin,xmax
      call flush( lp)
      call zaiowr(saln(1,1,k),ip,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'salin   ',nstep,time(3),k,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'salin   ',nstep,time(3),k,coord,xmin,xmax
      call flush( lp)
      if     (sigver.eq.0) then
        call zaiowr(th3d(1,1,k),ip,.true.,
     &              xmin,xmax, nop, .false.)
        write (nop,117) 'density ',nstep,time(3),k,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'density ',nstep,time(3),k,coord,xmin,xmax
        call flush( lp)
      endif !sigver==0
      if     (trcout) then
        do ktr= 1,ntracr
          call zaiowr(trcr(1,1,k,ktr),ip,.true.,
     &                xmin,xmax, nop, .false.)
          if     (ktr.le.99) then
c ---       older postprocessing code will not recognise this
            write(ctype,'(a6,i2.2)') 'tracer',ktr
          else
            ctype = 'tracer  '
          endif
          write (nop,117) ctype,nstep,time(3),k,coord,xmin,xmax
          call flush(nop)
          write ( lp,117) ctype,nstep,time(3),k,coord,xmin,xmax
          call flush( lp)
        enddo !ktr
      endif !trcout
 75   continue
c
 117  format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7)
c
      close (unit=nop)
      call zaiocl(nop)
      return
      end

      subroutine putdat_tide(flnm, artype,time,lsteric,
     &                       iexpt,iversn,yrflag, thbase)
      use mod_mom6  ! HYCOM mom6 array interface
      use mod_za    ! HYCOM array I/O interface
c
      character        flnm*(*)
      double precision time(3)
      real             thbase
      integer          artype,iexpt,iversn,yrflag
      logical          lsteric
c
c --- write model fields.
c --- HYCOM 2.0 array I/O archive file.
c
      character*8 ctype
c
      real      coord,xmin,xmax
      integer   i,j,k,ktr,jversn,l,nop
      data nop/24/
c
c
c ---   output is in "*.[ab]"
c
      l = len_trim(flnm)
      if     (flnm(l-1:l).eq.'.a' .or. flnm(l-1:l).eq.'.b') then
        open (unit=nop,file=flnm(1:l-2)//'.b',form='formatted',
     &          status='new',action='write')
        call zaiopf(flnm(1:l-2)//'.a','new', nop)
      else
        open (unit=nop,file=flnm(1:l)//'.b',form='formatted',
     &          status='new',action='write')
        call zaiopf(flnm(1:l)//'.a','new', nop)
      endif
c
c --- header.
c
      jversn = max(iversn,20)
      if     (artype.eq.1) then
        write(nop,116) ctitle,jversn,iexpt,yrflag,idm,jdm
        write( lp,  *)
        write( lp,116) ctitle,jversn,iexpt,yrflag,idm,jdm
 116    format (a80/a80/a80/a80/
     &   i5,4x,'''iversn'' = hycom version number x10'/
     &   i5,4x,'''iexpt '' = experiment number x10'/
     &   i5,4x,'''yrflag'' = days in year flag'/
     &   i5,4x,'''idm   '' = longitudinal array size'/
     &   i5,4x,'''jdm   '' = latitudinal  array size'/
     &   'field       time step  model day',
     &   '  k  dens        min              max')
      elseif (artype.eq.2) then
        write(nop,118) ctitle,jversn,iexpt,yrflag,idm,jdm
        write( lp,  *)
        write( lp,118) ctitle,jversn,iexpt,yrflag,idm,jdm
 118    format (a80/a80/a80/a80/
     &   i5,4x,'''iversn'' = hycom version number x10'/
     &   i5,4x,'''iexpt '' = experiment number x10'/
     &   i5,4x,'''yrflag'' = days in year flag'/
     &   i5,4x,'''idm   '' = longitudinal array size'/
     &   i5,4x,'''jdm   '' = latitudinal  array size'/
     &   'field         no. recs  mean day',
     &   '  k  dens        min              max')
      else
        write( lp,"(/ a /)")
     &    'error in putdat - only artpe==1 and artype==2 allowed'
        stop
      endif

c
c --- surface fields.
c
      coord=0.0
c
      call zaiowr(montg,ip,.true.,
     &            xmin,xmax, nop, .false.)
      if     (sigver.gt.0) then
c ---   identify the equation of state on the first record
        write (nop,117) 'montg1  ',nstep,time(1),sigver,thbase,xmin,xmax
        call flush(nop)
        write ( lp,117) 'montg1  ',nstep,time(1),sigver,thbase,xmin,xmax
        call flush( lp)
      else
        write (nop,117) 'montg1  ',nstep,time(1),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'montg1  ',nstep,time(1),0,coord,xmin,xmax
        call flush( lp)
      endif !sigver
      call zaiowr(srfht,ip,.true.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'srfhgt  ',nstep,time(2),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'srfhgt  ',nstep,time(2),0,coord,xmin,xmax
      call flush( lp)
c
      if(lsteric) then
        call zaiowr(steric,ip,.true., xmin,xmax, nop, .false.)
        write (nop,117) 'steric  ',nstep,time(3),0,coord,xmin,xmax
        call flush(nop)
        write ( lp,117) 'steric  ',nstep,time(3),0,coord,xmin,xmax
        call flush( lp)
      endif
c
c --- depth averaged fields (no mask, in case this is a subregion archive).
c
      call zaiowr(ubaro,iu,.false.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'u_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
      call zaiowr(vbaro,iv,.false.,
     &            xmin,xmax, nop, .false.)
      write (nop,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush(nop)
      write ( lp,117) 'v_btrop ',nstep,time(3),0,coord,xmin,xmax
      call flush( lp)
c
 117  format (a8,' =',i11,f11.3,i3,f7.3,1p2e16.7)
c
      close (unit=nop)
      call zaiocl(nop)
      return
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
Jump to Timeline #Navigation Navigation

  The Facilities And Campus Enhancement Initiative is passed by the UC Davis student body in February.
  Davis Botanical Society established.
  Neal Van Alfen becomes deans dean of the College of Agricultural and Environmental Sciences.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!
! darcy_weisbach_three_regimes.f90
! Contributors: Pi-Yueh Chuang <pychuang@gwu.edu>
!
! Distributed under terms of the BSD 3-Clause license.
!

!> @brief Calculate coefficients with the combination of three models.
module darcy_weisbach_three_regimes_module
    use darcy_weisbach_abstract_module
    implicit none
    private
    public:: DarcyWeisbachThreeRegimes

    !> @brief Darcy-Weisbach with cell-wide coefficients
    type, extends(DarcyWeisbachBase):: DarcyWeisbachThreeRegimes
        private ! variables
        !> @brief Keep the underlying roughness file name.
        character(len=255):: filename
        !> @brief Number of cells in x direction
        integer(kind=4):: mx
        !> @brief Number of cells in y direction
        integer(kind=4):: my
        !> @brief xlower of the coefficient grid provided.
        real(kind=8):: xlower
        !> @brief ylower of the coefficient grid provided.
        real(kind=8):: ylower
        !> @brief xupper of the coefficient grid provided.
        real(kind=8):: xupper
        !> @brief yupper of the coefficient grid provided.
        real(kind=8):: yupper
        !> @brief Cell size (assume uniform in both directions).
        real(kind=8):: cellsize
        !> @brief No data value.
        real(kind=8):: nodatavalue
        !> @brief Default roughness for regions uncovered by the file
        real(kind=8):: default_roughness = 0D0
        !> @brief Coefficient array.
        real(kind=8), allocatable, dimension(:, :):: roughness

        contains ! member functions
        !> @brief Initialize.
        procedure:: init_from_funit => init_from_funit_three_regimes
        !> @brief Underlying outputing.
        procedure:: write_data => write_data_three_regimes
        !> @brief Getting the coefficient of a single cell.
        procedure:: get_coefficient => get_coefficient_three_regimes
        !> @brief Read real coefficient file (in Esri ASCII format)
        procedure:: read_roughness_file
        !> @bried Destructor
        final:: destructor_three_regimes
    end type DarcyWeisbachThreeRegimes

contains

    ! implementation of init_from_funit_three_regimes
    subroutine init_from_funit_three_regimes(this, funit, kin_vis)
        class(DarcyWeisbachThreeRegimes), intent(inout):: this
        integer(kind=4), intent(in):: funit
        real(kind=8), intent(in):: kin_vis

        this%name = "Three-regime Darcy-Weisbach"
        this%nu = kin_vis
        read(funit, *) this%friction_tol
        read(funit, *) this%dry_tol
        read(funit, *) this%filename
        read(funit, *) this%default_roughness
        call this%read_roughness_file()
    end subroutine init_from_funit_three_regimes

    ! implementation of write_data_three_regimes
    subroutine write_data_three_regimes(this, iounit, iotype, v_list, stat, msg)
        class(DarcyWeisbachThreeRegimes), intent(in):: this
        integer(kind=4), intent(in):: iounit
        character(*), intent(in)::iotype
        integer(kind=4), intent(in):: v_list(:)
        integer(kind=4), intent(out):: stat
        character(*), intent(inout):: msg
        character:: n, t

        n = new_line(t) ! n is the "new line" character
        t = achar(9) ! t is the character for a tab

        write(iounit, *, iostat=stat, iomsg=msg) &
            n, this%friction_tol, t, t, "=: ", "friction_tol", &
            n, this%dry_tol, t, t, "=: ", "dry_tol", &
            n, this%filename, t, t, "=: filename", t, t, "# (coefficient file)", &
            n, this%default_roughness, t, t, "=: default_coefficien", t, t, &
            "# (coefficient for cells uncovered by the file)"

    end subroutine write_data_three_regimes

    ! implementation of get_coefficient_three_regimes
    function get_coefficient_three_regimes(this, x, y, q) result(coef)
        class(DarcyWeisbachThreeRegimes), intent(in):: this
        real(kind=8), intent(in):: x, y, q(3)
        real(kind=8):: coef, roughness, Re
        integer(kind=4):: i, j

        ! *********************************************************************
        ! Note: the following code is based on Re being defined by hydraulic
        !       radius, not hydraulic diameter
        ! *********************************************************************

        ! calculate local Reynolds number (defined by hydraulic radius, i.e., h)
        Re = dsqrt(q(2)**2+q(3)**2) / this%nu ! depth h is included in q(2) & q(3)

        ! static cell, no friction
        if (dabs(Re) <= 1D-7) return

        ! if it is laminar, we don't need roughness
        if (Re <= 5D2) then
            coef = 24D0 / Re
            return ! exit this function
        endif

        ! if it is not laminar but smooth turbulent regime, we don't need roughness
        if (Re <= 1250D0) then
            coef = 0.224D0 / (Re**0.25)
            return
        endif

        ! *********************************************************************
        ! For Re > 1250, we need roughness for Swamee & Jain model.
        ! *********************************************************************

        ! initialize roughness with default value
        roughness = this%default_roughness

        ! when the coordinate is covered by the provided roughness file
        if ((x >= this%xlower) .and. (x < this%xupper)) then
            if ((y >= this%ylower) .and. (y < this%yupper)) then
                ! TODO: should we use at least linear interpolation?

                i = int((x-this%xlower)/this%cellsize) + 1
                j = int((y-this%ylower)/this%cellsize) + 1
                roughness = this%roughness(i, j)
            endif
        endif

        ! Swamee & Jain (note we use hydraulic radius, not hydraulic diameter)
        coef = roughness/(14.8*q(1))+1.6483821394207454/(Re**0.9)
        coef = dlog10(coef)
        coef = coef * coef
        coef = 0.25D0 / coef

    end function get_coefficient_three_regimes

    ! implementation of read_roughness_file
    subroutine read_roughness_file(this, filename)
        use utility_module, only: parse_values
        class(DarcyWeisbachThreeRegimes), intent(inout):: this
        character(len=*), intent(in), optional:: filename
        integer(kind=4), parameter:: funit = 253
        character(len=255):: line
        integer(kind=4):: n_values, i, j
        real(kind=8):: values(10)

        if (present(filename)) then
            open(unit=funit, file=filename, action="read", &
                status="old", form="formatted")
        else
            open(unit=funit, file=this%filename, action="read", &
                status="old", form="formatted")
        endif

        ! mx
        read(funit, "(A)") line
        call parse_values(line, n_values, values)
        this%mx = idnint(values(1))

        ! my
        read(funit, "(A)") line
        call parse_values(line, n_values, values)
        this%my = idnint(values(1))

        ! xlower
        read(funit, "(A)") line
        call parse_values(line, n_values, values)
        this%xlower = values(1)

        ! ylower
        read(funit, "(A)") line
        call parse_values(line, n_values, values)
        this%ylower = values(1)

        ! cellsize
        read(funit, "(A)") line
        call parse_values(line, n_values, values)
        this%cellsize = values(1)

        ! nodatavalue
        read(funit, "(A)") line
        call parse_values(line, n_values, values)
        this%nodatavalue = values(1)

        ! calculate xupper and yupper
        this%xupper = this%xlower + this%cellsize * this%mx
        this%yupper = this%ylower + this%cellsize * this%my

        ! allocate coefficients
        allocate(this%roughness(this%mx, this%my))

        ! read coefficients
        do j = this%my, 1, -1
            read(funit, *) this%roughness(:, j)
        enddo

        ! handle missing data
        if (any(this%roughness == this%nodatavalue)) then
            where(this%roughness == this%nodatavalue) this%roughness = 0D0

            write(*, *) "WARNING: missing data found in the roughness file. &
                Set these data to zero automatically."
        endif

        close(funit)

        ! TODO: NetCDF file
    end subroutine read_roughness_file

    ! implementation of destructor_three_regimes
    subroutine destructor_three_regimes(this)
        type(DarcyWeisbachThreeRegimes), intent(inout):: this
        this%name = ''
        this%filename = ''
        this%mx = 0
        this%my = 0
        this%xlower = 0D0
        this%ylower = 0D0
        this%xupper = 0D0
        this%yupper = 0D0
        this%cellsize = 0D0
        this%nodatavalue = 0D0
        this%default_roughness = 0D0

        if (allocated(this%roughness)) deallocate(this%roughness)
    end subroutine destructor_three_regimes

end module darcy_weisbach_three_regimes_module

max_stars_repo_name: kbarnhart/geoclaw-landspill

max_issues_repo_name: kbarnhart/geoclaw-landspill

max_forks_repo_name: kbarnhart/geoclaw-landspill


text:
      SUBROUTINE ksone(data,n,func,d,prob)
      INTEGER n
      REAL d,data(n),func,prob
      EXTERNAL func
CU    USES probks,sort
      INTEGER j
      REAL dt,en,ff,fn,fo,probks
      call sort(n,data)
      en=n
      d=0.
      fo=0.
      do 11 j=1,n
        fn=j/en
        ff=func(data(j))
        dt=max(abs(fo-ff),abs(fn-ff))
        if(dt.gt.d)d=dt
        fo=fn
11    continue
      en=sqrt(en)
      prob=probks((en+0.12+0.11/en)*d)
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/ksone.for


text:
less_toxic(x1,v1).
less_toxic(cc1,d1).
less_toxic(q1,u1).
less_toxic(h1,c1).
less_toxic(v1,f1).
less_toxic(hh1,p1).
less_toxic(bb1,c1).
less_toxic(bb1,jj1).
less_toxic(y1,cc1).
less_toxic(m1,u1).
less_toxic(ff1,h1).
less_toxic(dd1,d1).
less_toxic(kk1,e1).
less_toxic(ee1,p1).
less_toxic(f1,d1).
less_toxic(ee1,h1).
less_toxic(bb1,d1).
less_toxic(aa1,d1).
less_toxic(ee1,a1).
less_toxic(ff1,l1).
less_toxic(n1,v1).
less_toxic(j1,jj1).
less_toxic(g1,c1).
less_toxic(r1,u1).
less_toxic(bb1,h1).
less_toxic(w1,i1).
less_toxic(q1,z1).
less_toxic(j1,aa1).
less_toxic(bb1,e1).
less_toxic(n1,aa1).
less_toxic(m1,d1).
less_toxic(ff1,f1).
less_toxic(g1,dd1).
less_toxic(ff1,e1).
less_toxic(b1,jj1).
less_toxic(j1,z1).
less_toxic(j1,f1).
less_toxic(x1,aa1).
less_toxic(kk1,b1).
less_toxic(x1,h1).
less_toxic(ee1,u1).
less_toxic(s1,c1).
less_toxic(o1,e1).
less_toxic(hh1,b1).
less_toxic(kk1,z1).
less_toxic(b1,u1).
less_toxic(h1,d1).
less_toxic(ff1,c1).
less_toxic(ee1,f1).
less_toxic(k1,v1).
less_toxic(a1,c1).
less_toxic(ll1,p1).
less_toxic(kk1,dd1).
less_toxic(x1,d1).
less_toxic(k1,a1).
less_toxic(bb1,v1).
less_toxic(n1,w1).
less_toxic(s1,t1).
less_toxic(n1,u1).
less_toxic(g1,cc1).
less_toxic(k1,h1).
less_toxic(w1,f1).
less_toxic(hh1,f1).
less_toxic(n1,e1).
less_toxic(ee1,d1).
less_toxic(m1,t1).
less_toxic(s1,f1).
less_toxic(kk1,l1).
less_toxic(t1,p1).
less_toxic(s1,i1).
less_toxic(g1,i1).
less_toxic(s1,p1).
less_toxic(ii1,u1).
less_toxic(b1,t1).
less_toxic(o1,t1).
less_toxic(r1,v1).
less_toxic(dd1,f1).
less_toxic(r1,aa1).
less_toxic(hh1,cc1).
less_toxic(x1,e1).
less_toxic(q1,i1).
less_toxic(n1,h1).
less_toxic(o1,v1).
less_toxic(ii1,f1).
less_toxic(n1,dd1).
less_toxic(y1,i1).
less_toxic(k1,d1).
less_toxic(j1,i1).
less_toxic(bb1,w1).
less_toxic(ee1,v1).
less_toxic(kk1,a1).
less_toxic(n1,p1).
less_toxic(x1,a1).
less_toxic(z1,f1).
less_toxic(r1,i1).
less_toxic(hh1,h1).
less_toxic(w1,p1).
less_toxic(w1,z1).
less_toxic(s1,z1).
less_toxic(jj1,d1).
less_toxic(b1,v1).
less_toxic(v1,p1).
less_toxic(r1,c1).
less_toxic(x1,i1).
less_toxic(y1,v1).
less_toxic(z1,d1).
less_toxic(r1,jj1).
less_toxic(ee1,w1).
less_toxic(r1,f1).
less_toxic(kk1,i1).
less_toxic(w1,u1).
less_toxic(m1,b1).
less_toxic(hh1,dd1).
less_toxic(k1,i1).
less_toxic(j1,b1).
less_toxic(m1,aa1).
less_toxic(ii1,v1).
less_toxic(m1,a1).
less_toxic(ll1,e1).
less_toxic(cc1,z1).
less_toxic(e1,p1).
less_toxic(x1,cc1).
less_toxic(ll1,u1).
less_toxic(s1,l1).
less_toxic(n1,i1).
less_toxic(w1,d1).
less_toxic(e1,c1).
less_toxic(q1,p1).
less_toxic(x1,t1).
less_toxic(kk1,jj1).
less_toxic(cc1,h1).
less_toxic(f1,c1).
less_toxic(h1,p1).
less_toxic(hh1,aa1).
less_toxic(bb1,aa1).
less_toxic(hh1,a1).
less_toxic(ee1,aa1).
less_toxic(ff1,aa1).
less_toxic(hh1,i1).
less_toxic(cc1,a1).
less_toxic(g1,v1).
less_toxic(ee1,jj1).
less_toxic(y1,b1).
less_toxic(o1,z1).
less_toxic(bb1,i1).
less_toxic(j1,t1).
less_toxic(s1,a1).
less_toxic(j1,e1).
less_toxic(p1,d1).
less_toxic(ff1,dd1).
less_toxic(o1,c1).
less_toxic(j1,l1).
less_toxic(p1,f1).
less_toxic(r1,d1).
less_toxic(bb1,f1).
less_toxic(ff1,jj1).
less_toxic(o1,aa1).
less_toxic(y1,w1).
less_toxic(bb1,cc1).
less_toxic(cc1,jj1).
less_toxic(j1,a1).
less_toxic(k1,z1).
less_toxic(k1,jj1).
less_toxic(cc1,e1).
less_toxic(j1,d1).
less_toxic(n1,f1).
less_toxic(m1,dd1).
less_toxic(s1,u1).
less_toxic(q1,aa1).
less_toxic(cc1,c1).
less_toxic(o1,b1).
less_toxic(cc1,t1).
less_toxic(ll1,z1).
less_toxic(hh1,c1).
less_toxic(aa1,f1).
less_toxic(o1,f1).
less_toxic(ii1,d1).
less_toxic(k1,dd1).
less_toxic(ii1,z1).
less_toxic(ff1,cc1).
less_toxic(y1,a1).
less_toxic(b1,e1).
less_toxic(w1,e1).
less_toxic(o1,u1).
less_toxic(bb1,t1).
less_toxic(q1,dd1).
less_toxic(y1,h1).
less_toxic(x1,p1).
less_toxic(g1,a1).
less_toxic(n1,c1).
less_toxic(ee1,t1).
less_toxic(j1,p1).
less_toxic(hh1,z1).
less_toxic(x1,w1).
less_toxic(ii1,a1).
less_toxic(y1,f1).
less_toxic(ii1,w1).
less_toxic(ll1,l1).
less_toxic(ff1,p1).
less_toxic(kk1,c1).
less_toxic(n1,b1).
less_toxic(p1,c1).
less_toxic(ff1,t1).
less_toxic(g1,aa1).
less_toxic(ee1,c1).
less_toxic(q1,w1).
less_toxic(s1,b1).
less_toxic(s1,cc1).
less_toxic(hh1,u1).
less_toxic(m1,h1).
less_toxic(bb1,dd1).
less_toxic(ll1,v1).
less_toxic(q1,e1).
less_toxic(j1,v1).
less_toxic(q1,d1).
less_toxic(aa1,p1).
less_toxic(ii1,dd1).
less_toxic(cc1,i1).
less_toxic(r1,w1).
less_toxic(l1,f1).
less_toxic(g1,u1).
less_toxic(ll1,d1).
less_toxic(w1,c1).
less_toxic(s1,h1).
less_toxic(ff1,u1).
less_toxic(ll1,aa1).
less_toxic(g1,e1).
less_toxic(jj1,f1).
less_toxic(ii1,l1).
less_toxic(g1,p1).
less_toxic(u1,c1).
less_toxic(k1,l1).
less_toxic(ll1,f1).
less_toxic(ff1,d1).
less_toxic(m1,l1).
less_toxic(q1,v1).
less_toxic(kk1,d1).
less_toxic(dd1,p1).
less_toxic(g1,f1).
less_toxic(m1,w1).
less_toxic(ff1,b1).
less_toxic(bb1,b1).
less_toxic(k1,p1).
less_toxic(r1,z1).
less_toxic(w1,a1).
less_toxic(w1,l1).
less_toxic(g1,d1).
less_toxic(o1,cc1).
less_toxic(x1,l1).
less_toxic(w1,v1).
less_toxic(kk1,p1).
less_toxic(t1,d1).
less_toxic(kk1,w1).
less_toxic(x1,u1).
less_toxic(kk1,h1).
less_toxic(ff1,i1).
less_toxic(q1,jj1).
less_toxic(n1,l1).
less_toxic(bb1,l1).
less_toxic(o1,a1).
less_toxic(ll1,dd1).
less_toxic(ii1,t1).
less_toxic(z1,c1).
less_toxic(y1,aa1).
less_toxic(n1,jj1).
less_toxic(m1,f1).
less_toxic(ii1,jj1).
less_toxic(q1,l1).
less_toxic(cc1,p1).
less_toxic(o1,p1).
less_toxic(k1,cc1).
less_toxic(ll1,h1).
less_toxic(k1,b1).
less_toxic(e1,f1).
less_toxic(ii1,cc1).
less_toxic(ll1,jj1).
less_toxic(bb1,z1).
less_toxic(u1,f1).
less_toxic(g1,l1).
less_toxic(s1,v1).
less_toxic(jj1,p1).
less_toxic(ff1,w1).
less_toxic(g1,z1).
less_toxic(cc1,v1).
less_toxic(ii1,aa1).
less_toxic(b1,p1).
less_toxic(m1,jj1).
less_toxic(e1,d1).
less_toxic(g1,h1).
less_toxic(ee1,z1).
less_toxic(g1,t1).
less_toxic(y1,p1).
less_toxic(v1,ii1).
less_toxic(d1,x1).
less_toxic(cc1,k1).
less_toxic(d1,p1).
less_toxic(t1,y1).
less_toxic(b1,j1).
less_toxic(l1,g1).
less_toxic(dd1,bb1).
less_toxic(p1,r1).
less_toxic(z1,cc1).
less_toxic(a1,x1).
less_toxic(c1,h1).
less_toxic(f1,y1).
less_toxic(jj1,b1).
less_toxic(u1,g1).
less_toxic(e1,q1).
less_toxic(p1,ff1).
less_toxic(u1,s1).
less_toxic(a1,j1).
less_toxic(w1,hh1).
less_toxic(v1,cc1).
less_toxic(b1,k1).
less_toxic(d1,y1).
less_toxic(b1,kk1).
less_toxic(d1,o1).
less_toxic(t1,ll1).
less_toxic(dd1,ff1).
less_toxic(z1,x1).
less_toxic(cc1,ii1).
less_toxic(i1,r1).
less_toxic(d1,f1).
less_toxic(aa1,cc1).
less_toxic(v1,m1).
less_toxic(c1,u1).
less_toxic(f1,ee1).
less_toxic(p1,aa1).
less_toxic(z1,y1).
less_toxic(h1,r1).
less_toxic(dd1,ii1).
less_toxic(jj1,j1).
less_toxic(u1,cc1).
less_toxic(c1,k1).
less_toxic(t1,w1).
less_toxic(p1,w1).
less_toxic(c1,x1).
less_toxic(l1,q1).
less_toxic(dd1,b1).
less_toxic(c1,ii1).
less_toxic(c1,y1).
less_toxic(c1,dd1).
less_toxic(v1,ll1).
less_toxic(i1,o1).
less_toxic(p1,dd1).
less_toxic(dd1,g1).
less_toxic(t1,q1).
less_toxic(c1,n1).
less_toxic(t1,j1).
less_toxic(cc1,s1).
less_toxic(t1,bb1).
less_toxic(c1,a1).
less_toxic(e1,x1).
less_toxic(z1,j1).
less_toxic(i1,k1).
less_toxic(jj1,kk1).
less_toxic(d1,m1).
less_toxic(jj1,ff1).
less_toxic(b1,ll1).
less_toxic(d1,ll1).
less_toxic(d1,h1).
less_toxic(aa1,q1).
less_toxic(v1,hh1).
less_toxic(u1,ii1).
less_toxic(p1,kk1).
less_toxic(d1,n1).
less_toxic(t1,m1).
less_toxic(l1,m1).
less_toxic(a1,g1).
less_toxic(w1,s1).
less_toxic(c1,v1).
less_toxic(c1,aa1).
less_toxic(i1,y1).
less_toxic(d1,r1).
less_toxic(t1,r1).
less_toxic(t1,g1).
less_toxic(p1,cc1).
less_toxic(p1,s1).
less_toxic(cc1,kk1).
less_toxic(a1,k1).
less_toxic(h1,w1).
less_toxic(dd1,kk1).
less_toxic(jj1,bb1).
less_toxic(u1,r1).
less_toxic(d1,k1).
less_toxic(e1,ll1).
less_toxic(z1,m1).
less_toxic(p1,e1).
less_toxic(u1,w1).
less_toxic(i1,s1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text: #redirect Varsity Theater


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

C=======================================================================
      SUBROUTINE CPUMVU(A,B,L)
      IMPLICIT INTEGER (A-Z)
      DIMENSION A(*),B(*)

      DO 2240 I = 1,L
         B(I) = A(I)
 2240 CONTINUE
      RETURN

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      SUBROUTINE WRITCP(LUNIT)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    WRITCP
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 2002-05-14
C
C ABSTRACT: THIS SUBROUTINE SHOULD ONLY BE CALLED WHEN LOGICAL UNIT
C   LUNIT HAS BEEN OPENED FOR OUTPUT OPERATIONS.  IT NOW SIMPLY CALLS
C   BUFR ARCHIVE LIBRARY SUBROUTINE CMPMSG TO TOGGLE ON MESSAGE
C   COMPRESSION, FOLLOWED BY A CALL TO WRITSB TO PACK UP THE CURRENT
C   SUBSET WITHIN MEMORY AND TRY TO ADD IT TO THE COMPRESSED BUFR
C   MESSAGE THAT IS CURRENTLY OPEN WITHIN MEMORY FOR THIS LUNIT,
C   FOLLOWED BY ANOTHER CALL TO CMPMSG TO TOGGLE OFF MESSAGE
C   COMPRESSION.  THIS SUBROUTINE USES THE SAME INPUT AND OUTPUT
C   PARAMETERS AS WRITSB.
C
C PROGRAM HISTORY LOG:
C 2002-05-14  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED
C                           DOCUMENTATION (INCLUDING HISTORY); OUTPUTS
C                           MORE COMPLETE DIAGNOSTIC INFO WHEN ROUTINE
C                           TERMINATES ABNORMALLY
C 2005-03-09  J. ATOR    -- MODIFIED TO USE CMPMSG AND WRITSB
C
C USAGE:    CALL WRITCP (LUNIT)
C   INPUT ARGUMENT LIST:
C     LUNIT    - INTEGER: FORTRAN LOGICAL UNIT NUMBER FOR BUFR FILE
C
C REMARKS:
C    THIS ROUTINE CALLS:        CMPMSG   WRITSB
C    THIS ROUTINE IS CALLED BY: None
C                               Normally called only by application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      CALL CMPMSG('Y')

      CALL WRITSB(LUNIT)

      CALL CMPMSG('N')

      RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak


text:
great(ii1,l1).
great(x1,s1).
great(h1,z1).
great(l1,ee1).
great(l1,kk1).
great(jj1,s1).
great(k1,w1).
great(ii1,dd1).
great(y1,s1).
great(ff1,w1).
great(j1,q1).
great(k1,j1).
great(h1,cc1).
great(y1,ff1).
great(d1,bb1).
great(h1,k1).
great(e1,o1).
great(ii1,u1).
great(i1,n1).
great(ii1,bb1).
great(q1,t1).
great(c1,n1).
great(hh1,u1).
great(l1,b1).
great(hh1,w1).
great(p1,o1).
great(w1,n1).
great(u1,kk1).
great(l1,p1).
great(i1,p1).
great(f1,jj1).
great(z1,b1).
great(jj1,u1).
great(l1,g1).
great(hh1,v1).
great(e1,x1).
great(e1,bb1).
great(cc1,v1).
great(m1,o1).
great(ii1,n1).
great(d1,m1).
great(ii1,e1).
great(ee1,b1).
great(w1,t1).
great(aa1,p1).
great(a1,j1).


text:
!-------------------------------------------------------------------------------!
!                                                                               !

subroutine phenology(nfile,row,col,s,y,yfit,nlocal,startcutoff,minimum)

!  This routine extracts phenological parameters from processed data            !
!                                                                               !
!   1.  time for which the left edge has increased to the 20 % level            !
!   2.  time for which the right edge has decreased to the 20 % level           !
!   3.  length of the season                                                    !
!   4.  average of left and right minimum values (base level)                   !
!   5.  time for the peak of the season                                         !
!   6.  value for the peak                                                      !
!   7.  seasonal amplitude                                                      !
!   8.  approximate left slope                                                  !
!   9.  approximate right slope                                                 !
!  10.  large seasonal integral                                                 !
!  11.  small seasonal integral                                                 !
!                                                                               !
!  Authors:                                                                     !
!  Per J\"onsson, Malm\"o University, Sweden                                    !
!  e-mail per.jonsson@ts.mah.se                                                 !
!                                                                               !
!  Lars Eklundh, Lund University, Sweden                                        !
!  e-mail lars.eklundh@nateko.lu.se                                             !
!                                                                               !
!  Updated 1/11-2005                                                            !
!                                                                               !
!-------------------------------------------------------------------------------!

implicit none

integer :: row,col,nyear,nptperyear,npt,nenvi,printflag,debugflag,win(3),nlocal
integer :: i,j,n,n1,n2,s(4*nyear+2),nfile,minposleft,minposright
integer :: maxposcenter,pos(1),maxminsum,minimum
real*4 :: dummy(11)
double precision :: wfact,y(npt),yfit(npt),minvalleft,minvalright,maxvalcenter,amplitude
double precision :: baseval,poscenter,posright20,posleft20,posright80,posleft80
double precision :: seasonlength,leftderivative,rightderivative,largeintegral,smallintegral
double precision :: valleft20,valleft80,valright20,valright80,startcutoff

common /parameters/ nyear,nptperyear,npt,nenvi,printflag,debugflag,wfact,win

n = nlocal + 2

if (minimum == 0) then       ! the position s(2) corresponds to a maximum
  n1 = 4
  n2 = n-2
else                         ! the position s(2) corresponds to a minimum
  n1 = 3
  n2 = n-3
end if

!---- Dump row, column and number of seasons on file ---------------------------
!    (n-4)/2 = number of seasons for which we have phenological metrics

write(nfile)  row,col,(n-4)/2

do i = n1,n2,2
  minvalleft = minval(yfit(s(i-2):s(i)))
  pos(1:1) = minloc(yfit(s(i-2):s(i)))
  minposleft = pos(1) + s(i-2) - 1

  minvalright = minval(yfit(s(i):s(i+2)))
  pos(1:1) = minloc(yfit(s(i):s(i+2)))
  minposright = pos(1) + s(i) - 1

  maxvalcenter = maxval(yfit(s(i-1):s(i+1)))
  pos(1:1) = maxloc(yfit(s(i-1):s(i+1)))
  maxposcenter = pos(1) + s(i-1) - 1

  amplitude = maxvalcenter - minvalleft

!---- Find left 80% level. Step from left min to max --------------------------

  valleft80 = minvalleft + 0.8d0*amplitude;
  posleft80 = 0.d0
  do j = minposleft,maxposcenter
    if (yfit(j) > valleft80) then
      posleft80 = j - 1.d0 + (valleft80 - yfit(j-1))/(yfit(j) - yfit(j-1))
!      posleft80 = j - 0.5d0
      exit
    end if
  end do

!---- Find left X % level. Step from max to left min --------------------------

  valleft20 = minvalleft + startcutoff*amplitude
  posleft20 = 0.d0
  do j = maxposcenter,minposleft,-1
    if (yfit(j) < valleft20) then
      posleft20 = j + (valleft20 - yfit(j))/(yfit(j+1) - yfit(j))
!      posleft20 = j + 0.5d0
	  exit
    end if
  end do

  amplitude = maxvalcenter - minvalright

!---- Find right 80% level. Step from right min to max ------------------------

  valright80 = minvalright + 0.8d0*amplitude
  posright80 = 0.d0
  do j = minposright,maxposcenter,-1
    if (yfit(j) > valright80) then
	  posright80 = j + 1.d0 - (valright80 - yfit(j+1))/(yfit(j) - yfit(j+1))
!      posright80 = j + 0.5d0
	  exit
    end if
  end do

!---- Find right X % level. Step from max to right min ------------------------

  valright20 = minvalright + startcutoff*amplitude
  posright20 = 0.d0
  do j = maxposcenter,minposright
    if (yfit(j) < valright20) then
      posright20 = j - (valright20 - yfit(j))/(yfit(j-1) - yfit(j));
!      posright20 = j - 0.5d0
	  exit
    end if
  end do

!--- Comment: We have added 0.01 in the denominator for the derivatives to ensure that
!    we never divide by zero

  if (posleft20*posleft80*posright20*posright80 > 0.d0) then
    seasonlength = posright20 - posleft20                                   ! 3) season length
    baseval = (minvalleft + minvalright)/2.d0                               ! 4) basevalue
    poscenter = (posleft80 + posright80)/2.d0                               ! 5) time for peak
    amplitude = maxvalcenter - baseval                                      ! 7) seasonal amplitude
    leftderivative = (valleft80 - valleft20)/(0.001 + posleft80 - posleft20) ! 8) left derivative
    rightderivative = -(valright80 - valright20)/(-0.001 + posright80 - posright20)  ! 9) right derivative
    largeintegral = sum(yfit(floor(posleft20):ceiling(posright20)))         ! 10) large integral
    smallintegral = sum(yfit(floor(posleft20):ceiling(posright20))-baseval) ! 11) small integral
  else
    posleft20 = 0.d0         ! 1)
    posright20 = 0.d0        ! 2)
    seasonlength = 0.d0      ! 3)
    baseval = 0.d0           ! 4)
    poscenter = 0.d0         ! 5)
    maxvalcenter = 0.d0      ! 6)
    amplitude = 0.d0         ! 7)
    leftderivative = 0.d0    ! 8)
    rightderivative = 0.d0   ! 9)
    largeintegral = 0.d0     ! 10)
    smallintegral = 0.d0     ! 11)
  end if
  dummy(1) = posleft20
  dummy(2) = posright20
  dummy(3) =   seasonlength
  dummy(4) =   baseval
  dummy(5) =   poscenter
  dummy(6) =   maxvalcenter
  dummy(7) =   amplitude
  dummy(8) =   leftderivative
  dummy(9) =   rightderivative
  dummy(10) =   largeintegral
  dummy(11) =   smallintegral
  write(nfile) dummy
end do

end subroutine phenology

max_stars_repo_name: bucricket/projectMASpreprocess

max_issues_repo_name: bucricket/projectMASpreprocess

max_forks_repo_name: bucricket/projectMASpreprocess


text:
!==============================================================================!
  subroutine Comm_Mod_Read_Int(fh, num, disp)
!------------------------------------------------------------------------------!
!   Read single integer for parallel runs.                                     !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  integer :: fh    ! file handle
  integer :: num   ! number to write out
  integer :: disp  ! diplacement in bytes
!-----------------------------------[Locals]-----------------------------------!
  integer :: error
!==============================================================================!

  ! Set it at position disp (same as in Write counterpart)
  call Mpi_File_Set_View(fh,             &
                         disp,           &
                         MPI_INTEGER8,   &
                         MPI_INTEGER8,   &
                         'native',       &
                         MPI_INFO_NULL,  &
                         error)

  ! Read integer value
  call Mpi_File_Read(fh,                 &
                     num,                &
                     1,                  &
                     MPI_INTEGER8,       &
                     MPI_STATUS_IGNORE,  &
                     error)

  disp = disp + SIZE_INT

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
module probdata_module

!     Combustor variables
      double precision, save :: pamb, Tinit, uinit
      double precision, save :: Tin
      double precision, save :: Rfrac
      double precision, save :: swrlang
      double precision, save :: pertmag

      double precision, save :: uin

      double precision :: rho_l, rhoe_l

!     These help specify which specific problem
      integer        , save ::  probtype,idir

      double precision, save :: split(3)

end module probdata_module

max_stars_repo_name: ndeak/PeleC

max_issues_repo_name: ndeak/PeleC

max_forks_repo_name: ndeak/PeleC


text:
      subroutine read_tcurve (tfilename, numpars, ierr)
!***********************************************************************
!  Copyright, 2004,  The  Regents  of the  University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
!***********************************************************************
!D1
!D1 PURPOSE
!D1
!D1 Read type curves for parallel fracture dispersion interpolation
!D1 calculations.
!D1
!***********************************************************************
!D2
!D2 REVISION HISTORY
!D2
!D2 FEHM Version 2.20 [10086-STN-2.20-00]
!D2 Initial implementation: 28-Oct-02, Programmer: Z. Dash
!D2
!D2 $Log:   /pvcs.config/fehm90/src/read_tcurve.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:43:48   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:13:30   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!***********************************************************************
!D3
!D3 REQUIREMENTS TRACEABILITY
!D3
!D3 2.3.5 Cell-based particle-tracking module
!D3 2.3.6 Streamline particle-tracking module
!D3
!***********************************************************************
!D4
!D4 SPECIAL COMMENTS AND REFERENCES
!D4
!D4 Sandia data curves
!D4
!D4 Requirements from SDN: 10086-RD-2.20-00
!D4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
!D4   FEHM Application Version 2.20
!D4
!***********************************************************************

      use compfrac
      use comai, only : neq
      use comdi, only : nspeci
      implicit none
      character*3 output_flag
      integer i, ierr, j, k, l, m, mm, n, dummy, numpars, nump_max
      integer ncurves
      real*8 dum(3)
      real*8 xconvert
      character*100 tfilename
      character*4 dummy_string
      character*80 input_msg
      integer msg(10)
      integer nwds
      real*8 xmsg(10)
      integer imsg(10),ditnumlines,i2,index
      character*32 cmsg(10)

      numparams = numpars
      open (121, file = tfilename, status = 'old', action = 'read')

      sigma_low = 1.d10
      sigma_high = 0.d0
      omega_low = 1.d10
      omega_high = 0.d0
      par3_low = 1.d10
      par3_high = 0.d0

      read(121,'(a4)') dummy_string
      if(dummy_string(1:4).eq.'free') then
         read(121,*) (log_flag(i),i=1,numparams)
         read(121,*) (normal_param(i),i=1,numparams)
         read(121,*) curve_structure

c     Input used to be 1 parameter, now it's 2. This section of code
c     checks for the number of parameters input. If only 1,
c     set the second parameter weight_factor to the old hardwired
c     value of 1.e-3 to make code behave as before. Otherwise,
c     set the value

         if(curve_structure.gt.1) then
            backspace 121
            read(121,'(a)') input_msg
            call parse_string(input_msg, imsg, msg, xmsg, cmsg, nwds)
            curve_structure = imsg(1)
            if(nwds.gt.1) then

c     Second number could be input as an integer, so make sure the
c     code accounts for that

               if(msg(2).eq.1) then
                  weight_factor = imsg(2)
               else
                  weight_factor = xmsg(2)
               end if
            else
               weight_factor = 1.e-3
            end if
         end if

         allocate(itf_curve(nspeci,neq,curve_structure))
         allocate(wt_curve(nspeci,neq,curve_structure))
         itf_curve = 0
         wt_curve = 0.
         read(121,*) ncurves
         allocate(param1(ncurves))
         allocate(param2(ncurves))
         param1 = 0.
         param2 = 0.
         nump1 = ncurves
         nump2 = 1
         if (numparams .gt. 2) then
            allocate (param3 (ncurves))
            param3 = 0.
            d4 = 4
         else
            d4 = 1
         end if
         nump3 = 1

      else
         backspace 121
         curve_structure = 0
         read (121, *) nump1
         allocate (param1 (nump1))
         read (121,*) (param1(i), i= 1, nump1)
         do i = 1, nump1
            sigma_low(1) = min (sigma_low(1), param1(i))
            sigma_high(1) = max (sigma_high(1), param1(i))
         end do
         read (121, *) nump2
         allocate (param2 (nump2))
         read (121, *) (param2(i), i= 1, nump2)
         do i = 1, nump2
            omega_low(1) = min(omega_low(1), param2(i))
            omega_high(1) = max(omega_high(1), param2(i))
         end do
         if (numparams .gt. 2) then
            read (121, *) nump3
            allocate (param3 (nump3))
            read (121, *) (param3(i), i= 1, nump3)
            do i = 1, nump3
               par3_low(1) = min(par3_low(1), param3(i))
               par3_high(1) = max(par3_high(1), param3(i))
            end do
            d4 = 4
         else
            nump3 = 1
            allocate (param3 (nump3))
            param3 = 0.
            par3_low = 0.
            par3_high = 0.
            d4 = 1
         end if
      end if
      read (121, *) nump_max
      allocate (nump(nump1,nump2,nump3,d4))
      allocate (dtime(nump1,nump2,nump3,d4,nump_max),
     .     conc(nump1,nump2,nump3,d4,nump_max))
      dtime = 0.
      conc = 0.
!      backspace (121)
      do l = 1, d4
         do i = 1, nump1
            do j = 1, nump2
               do k = 1, nump3
                  read (121, *) nump(i,j,k,l), (dum(n),
     .                 n = 1, numparams)
                  if(curve_structure.gt.0) then
                     param1(i) = dum(1)
                     sigma_low(1) = min (sigma_low(1), dum(1))
                     sigma_high(1) = max (sigma_high(1), dum(2))
                     param2(i) = dum(2)
                     omega_low(1) = min (omega_low(1), dum(1))
                     omega_high(1) = max (omega_high(1), dum(2))
                     if(numparams.gt.2) then
                        param3(i) = max(1.d-10,dum(3))
                        par3_low(1) = min (par3_low(1), dum(3))
                        par3_high(1) = max (par3_high(1), dum(3))
                     else
                        par3_low = 0.
                        par3_high = 0.
                     end if
                  end if
                  do m = 1, nump(i,j,k,l)
                     read (121, *) dtime(i,j,k,l,m), conc(i,j,k,l,m)
                  enddo
                  if(numparams.gt.2) then
c                     xconvert = (param1(i)+param2(i))/
                     xconvert = param2(i)*(1.-param3(i))/
     2                    (param1(i)+param2(i)*param3(i))
                     dtime(i,j,k,l,1:nump(i,j,k,l)) =
     2                    (dtime(i,j,k,l,1:nump(i,j,k,l))-1.)/xconvert
                  end if
                  mm=1
                  do m=2,nump(i,j,k,l)
                     if(conc(i,j,k,l,m)-conc(i,j,k,l,mm)
     2                    .lt.-1.e-20 )then
                        write(ierr,*)'Decreasing data found in type ',
     .                       'curve, stop'
                        write(ierr,*)'point:',m,' conc=',
     .                       conc(i,j,k,l,m)
                        stop
                     end if
                     if(conc(i,j,k,l,m).ne.conc(i,j,k,l,mm))then
                        mm=mm+1
                        if(mm.ne.m)then
                           conc(i,j,k,l,mm)=conc(i,j,k,l,m)
                           dtime(i,j,k,l,mm)=dtime(i,j,k,l,m)
                        end if
                     end if
                  end do
                  nump(i,j,k,l)=mm
               enddo
            enddo
         enddo
      enddo
      sigma_low(2) = sigma_high(1)
      sigma_high(2) = sigma_low(1)
      omega_low(2) = omega_high(1)
      omega_high(2) = omega_low(1)
      par3_low(2) = par3_high(1)
      par3_high(2) = par3_low(1)
      ipartout = 0
      read(121,'(a3)',end=2000) output_flag
      if(output_flag(1:3).eq.'out') then
         ipartout = 1
         allocate(param_density(nump1,nump2,nump3))
         param_density = 0
      end if
 2000 continue
      close (121)

      pfrac_read = .true.
      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM