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.

      module Grids

!          module contains everything for specification of multiple grids
!          created  17 October 2002 by Leo Postma
!          modified somewhere  2003 by Jan van Beek: special of layered bed
!                   May        2011 by Leo postma  : merge special with standard version

!     contains the following derived types:
!          GridPointer            ! a set of information with respect to one grid pointer
!          GridPointerColl        ! a collection of these grid pointers

!     contains the following functions:
!          GridPointerCollFind    ! to search a Grid in the GridPointerColl ; returns the index or zero if not found
!          GridPointerCollAdd     ! to add a GridPointer to the collection ; returns the current size

!     contains the following subroutine:

      integer, parameter :: NAME_SIZE       = 20                 ! size of descriptive names
      integer, parameter :: MAX_NUM         =  5                 ! allocated per bunch

      integer, parameter :: BaseGrid        =  1                 ! implementation of an enumeration
      integer, parameter :: ProcessGrid     =  2                 !               type in Fortran
      integer, parameter :: BottomGrid      =  3
      integer, parameter :: AggregationFile =  4
      integer, parameter :: NolayGrid       =  5
      integer, parameter ::              NrGridTypes    =  5
      character*20        GridTypes( NrGridTypes )
      DATA GridTypes / 'BASEGRID', 'PROCESSGRID', 'BOTTOMGRID' , 'AGGREGATIONFILE' , 'NOLAY' /

!          this is the grid pointer itself

      type GridPointer
         character(len=NAME_SIZE)         :: name               ! name of the grid
         integer                          :: noseg              ! number of segments
         integer                          :: noseg_lay          ! number of segments per layer / 2D
         integer                          :: iref               ! grid reference nr
         character(len=NAME_SIZE)         :: name_ref           ! name of the reference grid
         integer                          :: itype              ! type of grid
         integer, pointer                 :: iarray      (:)    ! the pointer to reference the reference grid
         integer, pointer                 :: finalpointer(:)    ! the pointer to the final grid
         logical                          :: space_var_nolay    ! switch for space varying nr of layers
         integer                          :: nolay              ! nr of expandable layers
         integer, pointer                 :: nolay_var   (:)    ! space varying nr of layers if any
      end type GridPointer

!          this is the collection of the grid pointers

      type GridPointerColl
         type(GridPointer), pointer       :: Pointers(:)  => NULL() ! array with gridpointers
         integer                          :: maxsize      = 0       ! maximum size of the current array
         integer                          :: cursize      = 0       ! filled up to this size
         integer                          :: base_grid    = 0       ! index base grid in collection
         integer                          :: bottom_grid  = 0       ! index bottom grid in collection
      end type GridPointerColl

      contains

!          function to find a grid name in a collection of GridPointers

      function GridPointerCollFind( aGridPointerColl, name ) result ( iret )

         type(GridPointerColl)            :: aGridPointerColl
         character(LEN=*)                 :: name
         integer                          :: iret

         iret = 0
         do i = 1 , aGridPointerColl%cursize         ! search by name
            if ( aGridPointerColl%Pointers(i)%name .eq. name ) then
               iret = i
               return
            endif
         end do

      end function GridPointerCollFind

!          function to add to a collection of fileproperties

      function GridPointerCollAdd( aGridPointerColl , aGridPointer ) result ( cursize )

         type(GridPointerColl)               :: aGridPointerColl      ! the collection of GridPointers
         type(GridPointer)                   :: aGridPointer          ! the GridPointer to add to the collection
         integer                             :: cursize               ! return value the new current collection size
                                                                      ! and the index of the added GridPointer

         type(GridPointer), pointer          :: aGridPointerPnts(:)   ! should be a pointer for the resize operation

         if ( aGridPointerColl%cursize .eq. aGridPointerColl%maxsize ) then
            allocate ( aGridPointerPnts ( aGridPointerColl%maxsize + MAX_NUM ) )
            do i = 1 , aGridPointerColl%maxsize
               aGridPointerPnts(i) = aGridPointerColl%Pointers(i)        ! copies the pointers
            enddo
            if ( aGridPointerColl%maxsize .ne. 0 ) deallocate ( aGridPointerColl%Pointers )
            aGridPointerColl%Pointers => aGridPointerPnts                   ! attaches this new array of pointers
            aGridPointerColl%maxsize = aGridPointerColl%maxsize + MAX_NUM
         endif

         aGridPointerColl%cursize = aGridPointerColl%cursize + 1
         aGridPointerColl%Pointers(aGridPointerColl%cursize ) = aGridPointer
         cursize = aGridPointerColl%cursize
         return

      end function GridPointerCollAdd

      function GridWrite( ilun, aGrid ) result ( ierror )

         integer, intent(in)                :: ilun         ! unit number binary file with data
         type(GridPointer), intent(in)      :: aGrid        ! datastructure to be written
         integer                            :: ierror       !

         ierror  = 0

         write(ilun, err = 100 ) aGrid%name
         write(ilun, err = 100 ) aGrid%noseg
         write(ilun, err = 100 ) aGrid%noseg_lay
         write(ilun, err = 100 ) aGrid%iref
         write(ilun, err = 100 ) aGrid%name_ref
         write(ilun, err = 100 ) aGrid%itype
         write(ilun, err = 100 ) aGrid%finalpointer
         write(ilun, err = 100 ) aGrid%space_var_nolay
         write(ilun, err = 100 ) aGrid%nolay
         if ( aGrid%space_var_nolay ) then
            write(ilun, err = 100 ) aGrid%nolay_var
         endif

         return

  100    continue

         ierror = 1
         return

      end function GridWrite

      function GridRead( ilun, aGrid, noseg ) result ( ierror )

         integer, intent(in)                :: ilun         ! unit number binary file with data
         type(GridPointer), intent(out)     :: aGrid        ! datastructure to be filled
         integer, intent(in)                :: noseg        ! number of segments in base grid (perhaps remove this dependency by adding the length of the pointer to the structure)
         integer                            :: ierror       !

         ierror  = 0

         read(ilun, err = 100 ) aGrid%name
         read(ilun, err = 100 ) aGrid%noseg
         read(ilun, err = 100 ) aGrid%noseg_lay
         read(ilun, err = 100 ) aGrid%iref
         read(ilun, err = 100 ) aGrid%name_ref
         read(ilun, err = 100 ) aGrid%itype
         allocate ( aGrid%finalpointer(noseg) , stat = ierr_alloc)
         if ( ierr_alloc .ne. 0 ) then
            write(*,*) 'ERROR : allocating array in GridRead'
            call srstop(1)
         endif
         read(ilun, err = 100 ) aGrid%finalpointer
         read(ilun, err = 100 ) aGrid%space_var_nolay
         read(ilun, err = 100 ) aGrid%nolay
         if ( aGrid%space_var_nolay ) then
            allocate ( aGrid%nolay_var(aGrid%noseg_lay) , stat = ierr_alloc)
            if ( ierr_alloc .ne. 0 ) then
               write(*,*) 'ERROR : allocating array in GridRead'
               call srstop(1)
            endif
            read(ilun, err = 100 ) aGrid%nolay_var
         endif

         return

  100    continue

         ierror = 1
         return

      end function GridRead

      end module Grids


text:
c     ****************************************************************
c     *                                                              *
c     *                      subroutine read_damage                  *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 10/16/94                   *
c     *                   last modified : 07/28/95 (rhd)             *
c     *                                                              *
c     *     this subroutine reads in the damage routines from        *
c     *     a file.                                                  *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine read_damage( status, fileno, prec_fact )
      use global_data ! old common.main
c
      use elem_extinct_data
      use node_release_data
      use damage_data
c
      implicit integer (a-z)
c
      go to (100,200,300,400,500,600,700,800,900,1000), status
c
c                               read dam_state, dam_ifv, dam_blk_killed
c
 100  continue
      call rd2d( fileno, dam_ifv, prec_fact*mxedof, prec_fact*mxedof,
     &           num_kill_elem )
      call rdbk( fileno, dam_state, num_kill_elem )
      call rdbk( fileno, dam_blk_killed, nelblk )
      go to 9999
c
c                               read dam_print list
c
 200  continue
      call rdbk( fileno, dam_print_list, num_print_list )
      call rdbk( fileno, old_mises, num_print_list*prec_fact )
      call rdbk( fileno, old_mean, num_print_list*prec_fact )
      go to 9999
c
c                               read kill_order_list
c
 300  continue
      call rdbk( fileno, kill_order_list, num_kill_order_list )
      go to 9999
c
c                               read dam_node_elecnt
c
 400  continue
      call rdbk( fileno, dam_node_elecnt, nonode )
      go to 9999
c
c                               read dam_face_nodes,
c                               dam_dbar_elems
c
 500  continue
      call rdbk( fileno, dam_face_nodes, 4*num_kill_elem )
      call rdbk( fileno, dam_dbar_elems, (2*prec_fact)*num_kill_elem )
      go to 9999
c
c                               read all of the node release crack
c                               growth variables
c
 600  continue
      call rdbk( fileno, crack_plane_nodes, num_crack_plane_nodes )
      call rdbk( fileno, inv_crkpln_nodes, nonode )
      call rdbk( fileno, num_neighbors, num_crack_plane_nodes )
      call rd2d( fileno, neighbor_nodes, mxconn, mxconn,
     &           num_crack_plane_nodes )
      call rdbk( fileno, crkpln_nodes_state, num_crack_plane_nodes )
      call rdbk( fileno, crkpln_nodes_react,
     &           num_crack_plane_nodes * prec_fact )
      call rd2d( fileno, crack_front_nodes, num_crack_plane_nodes,
     &           num_crack_plane_nodes,2 )
      go to 9999
c
c                               read the info for traction separation
c
 700  continue
      call rdbk( fileno, node_release_frac,
     &           num_crack_plane_nodes * prec_fact )
      go to 9999
c
c                               read the old angles for overshoot control
c
 800  continue
      call rd2d( fileno, old_angles_at_front, num_crack_plane_nodes
     &       * prec_fact, num_crack_plane_nodes * prec_fact, mxconn )
      go to 9999
c
c                               read the old damage values needed for
c                               load size control during crack growth
c
 900  continue
      if( crack_growth_type .eq. 1 ) then
         call rdbk( fileno, old_porosity, num_kill_elem * prec_fact )
         call rdbk( fileno, del_poros, mxstp_store * prec_fact )
      else if( crack_growth_type .eq. 3 ) then
         call rdbk( fileno, old_plast_strain, num_kill_elem *
     &              prec_fact )
      else if( crack_growth_type .eq. 4 ) then
         call rdbk( fileno, old_deff, num_kill_elem * prec_fact )
         call rdbk( fileno, del_deff, mxstp_store * prec_fact )
      end if
      go to 9999
c
c                               read the data for constant front growth
c
 1000 continue
      call rdbk( fileno, master_nodes, num_crack_fronts )
      call rd2d( fileno, crack_front_list,
     &           num_crack_fronts*num_nodes_grwinc,
     &           num_crack_fronts*num_nodes_grwinc, num_nodes_thick )
      call rd2d( fileno, master_lines, num_crack_fronts,
     &           num_crack_fronts, num_nodes_back + 1 )
      go to 9999
c
c
 9999 continue
c
      return
      end

ext: f


text:
-sv

../rtl/shared/params.sv
../rtl/shared/structs.sv
../rtl/vector/vmacros.sv
../rtl/vector/vstructs.sv

../rtl/shared/and_or_mux.sv
../rtl/shared/arbiter.sv
../rtl/shared/eb_buff_generic.sv
../rtl/shared/eb_one_slot.sv
../rtl/shared/eb_two_slot.sv
../rtl/shared/fifo_dual_ported.sv
../rtl/shared/fifo_duth.sv
../rtl/shared/fifo_flush.sv
../rtl/shared/fifo_overflow.sv
../rtl/shared/fifo_initialized.sv
../rtl/shared/lru.sv
../rtl/shared/lru_two.sv
../rtl/shared/lru_more.sv
../rtl/shared/onehot_detect.sv
../rtl/shared/rr_arbiter.sv
../rtl/shared/sram.sv

../rtl/shared/data_cache.sv
../rtl/shared/ld_st_buffer.sv
../rtl/shared/wait_buffer.sv
../rtl/shared/data_operation.sv
../rtl/vector/vdata_operation.sv
../rtl/vector/vld_st_buffer.sv

../rtl/vector/vector_top.sv

../rtl/vector/vrrm.sv
../rtl/vector/vrat.sv

../rtl/vector/vmu.sv
../rtl/vector/vmu_st_eng.sv
../rtl/vector/vmu_ld_eng.sv
../rtl/vector/vmu_tp_eng.sv

../rtl/vector/vis.sv
../rtl/vector/vrf.sv

../rtl/vector/vex.sv
../rtl/vector/vex_pipe.sv
../rtl/vector/v_int_alu.sv
../rtl/vector/v_fp_alu.sv

max_stars_repo_name: ic-lab-DUTH/RISC-V-Vector

max_issues_repo_name: ic-lab-DUTH/RISC-V-Vector

max_forks_repo_name: ic-lab-DUTH/RISC-V-Vector


text:
program mce
    ! Markov Chain Estimator
    ! Victhor S. Sartório

    ! Given an argument containing a filepath, it reads that file
    ! and estimates a transition matrix for it. Optionally, if given
    ! a second argument containing a new filepath, it saves the
    ! estimated transition matrix to it.

    use arguments, only: parse_arguments
    use file_handler, only: open_files
    use estimator, only: exploratory_analysis, estimate
    use simulator, only: simulate
    use reporter, only: report

    implicit none

    call parse_arguments

    call open_files

    call exploratory_analysis
    call estimate

    call simulate

    call report
end program mce

max_stars_repo_name: vct-sartor/MCE

max_issues_repo_name: vsart/MCE

max_forks_repo_name: vsart/MCE


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

!
!****************************************************************************************
!       For scribe I/O
!	Read in (x,y,z) from station.bp or station.sta (sta format; x,y, with z),
!       where z is either distance from F.S. or z-coord. If below bottom or above F.S.,
!       const. extrapolation is used, except if z=1.e10, in which case
!       depth averaged value will be calculated (for 3D vars).
!       Output time series for 3D variables (surface values for 2D variables), DEFINED AT NODES OR
!       ELEMENTS.

!       Inputs:
!              (1) screen;
!              (2) station.bp or station.sta
!              (3) vgrid.in: in this dir or ../
!              (4) out2d*.nc
!              (4) nc outputs for that variable(tri-quad)

!       Outputs: fort.18; ; fort.20 - local depth for each pt.
!       For ics=2 (e.g. for lon/lat), use nearest node for output
!
!   ifort -mcmodel=medium -assume byterecl -CB -O2 -o read_output10_xyz.exe ../UtilLib/extract_mod2.f90  ../UtilLib/compute_zcor.f90 ../UtilLib/pt_in_poly_test.f90 read_output10_xyz.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff
!****************************************************************************************
!
      program read_out
      use netcdf
      use extract_mod2
      use compute_zcor
      use pt_in_poly_test

      character(len=30) :: file63,file62,varname
      character(len=12) :: it_char
!      character(len=48) :: version,variable_nm,variable_dim
      logical :: lexist
      dimension swild(3)
      integer, allocatable :: i34(:),elnode(:,:),node3(:,:),kbp(:),iep(:),irank_read(:),kbp00(:)
      real*8,allocatable :: timeout(:),xnd(:),ynd(:)
      real,allocatable :: dp(:),ztot(:),sigma(:),sigma_lcl(:,:),outvar(:,:), &
     &out(:,:,:,:),out2(:,:),eta2(:),arco(:,:),ztmp(:),x00(:),y00(:), &
     &out3(:),z00(:),rl2min(:),dep(:),ztmp2(:,:)
      integer :: nodel(3),dimids(100),idims(100)
      integer :: char_len,start_2d(2),start_3d(3),start_4d(4), &
     &count_2d(2),count_3d(3),count_4d(4)
!      real*8 :: h0

!      print*, 'Do you work on uncombined (0) or combined (1) nc?'
!      read(*,*)icomb
!      if(icomb/=0.and.icomb/=1) stop 'Unknown icomb'

!...  Set max array size for system memory
!      print*, 'Recommendation: for uncombined nc, specify max array size (e.g., <=2.e9);'
!      print*, 'for combined nc, specify # of records to read each time.'
!      print*, 'Do you want to specify max array size (1) or # of records (2)'
!      read(*,*)ispec_max
!      if(ispec_max==1) then
!        print*, 'Input max array size (e.g., <=2.e9):'
!        read(*,*)max_array_size
!        print*, 'max_array_size read in=',max_array_size
!      else
!        print*, 'Input # of records:'
!        read(*,*)nrec3
!      endif
!      ispec_max=2
!      print*, 'Input # of records:'
!      read(*,*)nrec3
!      nrec3=1

      print*, 'Input extraction pts format (1: .bp; 2:.sta):'
      read(*,*)ibp
      if(ibp/=1.and.ibp/=2) stop 'Unknown format'

      print*, 'Input ics (1-linear interp; 2-nearest neighbor interp. 2 for node-based variables only! 2 is suggested for sub-meter resolution!):'
      read(*,*)ics

      print*, 'Input variable name to read from nc (e.g. elevation):'
      read(*,'(a30)')varname
      varname=adjustl(varname); len_var=len_trim(varname)

!      print*, 'Is the var 2D (1) or 3D (2)?'
!      read(*,*) i23d

      print*, 'Is the var node (1) or elem (2) based?'
      read(*,*) inode_elem

      print*, 'Input start and end file # to read:'
      read(*,*) iday1,iday2

      print*, 'Is the z-coord. in station.* relative to surface (1) or a fixed level (0)?'
      read(*,*) ifs
!'
      if(ibp==1) then !.bp format
        open(10,file='station.bp',status='old')
        read(10,*)
        read(10,*) nxy
      else !.sta format
        open(10,file='station.sta',status='old')
        read(10,*) nxy
      endif !ibp

      allocate(x00(nxy),y00(nxy),z00(nxy),rl2min(nxy),dep(nxy),stat=istat)
      if(istat/=0) stop 'Failed to allocate (1)'

      do i=1,nxy
        if(ibp==1) then !.bp format
          read(10,*)j,x00(i),y00(i),z00(i)
        else !.sta format
          read(10,*)
          read(10,*)x00(i),y00(i),z00(i)
        endif !ibp
      enddo !i
      close(10)

!...  Get basic info from out2d*.nc
      !Returned vars: ne,np,ns,nrec,[xnd ynd dp](np),
      !elnode,i34,nvrt,h0,dtout,kbp
      call get_dims(iday1,np,ne,ns,nvrt,h0)
      allocate(xnd(np),ynd(np),dp(np),kbp(np),i34(ne),elnode(4,ne),stat=istat)
      if(istat/=0) stop 'alloc (1)'
      call readheader(iday1,np,ne,ns,kbp,i34,elnode,nrec,xnd,ynd,dp,dtout)

      print*, 'After header:',ne,np,ns,nvrt,nrec,i34(ne), &
     &elnode(1:i34(ne),ne),h0,xnd(np),ynd(np),dp(np),dtout !,start_time

!...  Read in time records in segments for mem
      if(inode_elem==1) then !node based
        last_dim=np
      else !elem
        last_dim=ne
      endif

      allocate(timeout(nrec),ztot(nvrt),sigma(nvrt),sigma_lcl(nvrt,np),kbp00(np), &
     &outvar(nvrt,last_dim),eta2(np),out2(nxy,nvrt),out3(nxy), &
     &out(nxy,3,nvrt,2),node3(nxy,3),arco(nxy,3),iep(nxy),ztmp(nvrt),ztmp2(nvrt,3))
      outvar=-huge(1.0) !test mem

!     Read in vgrid
      call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b, &
     &theta_f,ztot,sigma,sigma_lcl,kbp)

!     Calculate kbp00
      if(ivcor==1) then
        kbp00=kbp
      else
        do i=1,np
          !Use large eta to get true bottom
          call zcor_SZ_single(dp(i),1.e8,h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma, &
     &ztmp(:),idry2,kbp00(i))
        enddo !i
      endif !ivcor

!...  Find parent element for (x00,y00) (global indices)
      iep=0
      if(ics==1) then !Cartesian
        do i=1,ne
          do l=1,nxy
            if(iep(l)/=0) cycle
            call pt_in_poly_single(i34(i),real(xnd(elnode(1:i34(i),i))), &
     &real(ynd(elnode(1:i34(i),i))),x00(l),y00(l),inside,arco(l,1:3),nodel)
            if(inside==1) then
              iep(l)=i
              !print*, 'Found:',l,arco(l,1:3),nodel
              node3(l,1:3)=elnode(nodel(1:3),i)
            endif !inside
          enddo !l; build pts

          ifl=0 !flag
          do l=1,nxy
            if(iep(l)==0) then
              ifl=1
              exit
            endif
          enddo !l
          if(ifl==0) exit
        enddo !i=1,ne
      else !lat/lon; needed for node-based only
        rl2min=1.e25 !min distance^2
        do ie=1,ne
          do j=1,i34(ie)
            i=elnode(j,ie)
            do l=1,nxy
              rl2=(xnd(i)-x00(l))**2+(ynd(i)-y00(l))**2
              if(rl2<rl2min(l)) then
                rl2min(l)=rl2
                iep(l)=ie
                node3(l,1:3)=i
                arco(l,1:3)=1./3
              endif
            enddo !l=1,nxy
          enddo !j
        enddo !i=1,np
      endif !ics

      do j=1,nxy
        if(iep(j)<=0) then
          print*, 'Cannot find a parent for pt:',j,x00(j),y00(j)
          stop
        endif
      enddo !j

!...  Time iteration
!...
      do iday=iday1,iday2
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      write(it_char,'(i12)')iday
      it_char=adjustl(it_char)
      leng=len_trim(it_char)
      file62='out2d_'//it_char(1:leng)//'.nc'
      iret=nf90_open(trim(adjustl(file62)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid4)
      !time is double
      iret=nf90_inq_varid(ncid4,'time',itime_id)
      iret=nf90_get_var(ncid4,itime_id,timeout,(/1/),(/nrec/))
      print*, 'time=',timeout !,trim(adjustl(file63))

      !Find nc file
      file63=varname(1:len_var)//'_'//it_char(1:leng)//'.nc'
      inquire(file=file63,exist=lexist)
      if(lexist) then
        i23d=2 !3D var
      else
        i23d=1 !2D
        file63=file62
      endif

      iret=nf90_open(trim(adjustl(file63)),OR(NF90_NETCDF4,NF90_NOWRITE),ncid)
      iret=nf90_inq_varid(ncid,varname(1:len_var),ivarid1)
      if(iret/=nf90_NoErr) stop 'Var not found'
      iret=nf90_Inquire_Variable(ncid,ivarid1,ndims=ndims,dimids=dimids)
      if(ndims>100) stop 'increase dimension of dimids & idims'
      do i=1,ndims
       iret=nf90_Inquire_Dimension(ncid,dimids(i),len=idims(i))
      enddo !i
      npes=idims(ndims-1) !np|ne|ns
      if(npes/=np.and.npes/=ne) stop 'can only handle node- or elem-based'
!'
      if(idims(ndims)/=nrec) stop 'last dim is not time'

!        iret=nf90_get_att(ncid,ivarid1,'i23d',i23d)
!        if(i23d<=0.or.i23d>6) stop 'wrong i23d'
!        if(i23d>3.and.ics==2) stop 'ics=2 with elem-based var'
!        iret=nf90_get_att(ncid,ivarid1,'ivs',ivs)
!        !print*, 'i23d:',i23d,ivs,idims(1:ndims)
!
      do irec=1,nrec
        !Get elev
        iret=nf90_inq_varid(ncid4,'elevation',itmp)
        start_2d(1)=1; start_2d(2)=irec
        count_2d(1)=np; count_2d(2)=1
        iret=nf90_get_var(ncid4,itmp,eta2,start_2d,count_2d)

!        if(ivs==1) then !scalar
        if(i23d==1) then !2D
          start_2d(1)=1; start_2d(2)=irec
          count_2d(1)=npes; count_2d(2)=1
          iret=nf90_get_var(ncid,ivarid1,outvar(1,1:npes),start_2d,count_2d)
        else !3D
          start_3d(1:2)=1; start_3d(3)=irec
          count_3d(1)=nvrt; count_3d(2)=npes; count_3d(3)=1
          iret=nf90_get_var(ncid,ivarid1,outvar(:,1:npes),start_3d,count_3d)
        endif

        !Available now: outvar(nvrt,np|ne), eta2(np)
        out2=0
        out3=0
        if(i23d==1) then !2D
          do i=1,nxy
            dep(i)=0
            do j=1,3 !nodes
              nd=node3(i,j)
              !Compute local depth
              dep(i)=dep(i)+arco(i,j)*dp(nd)
              if(inode_elem==1) then !node
                out2(i,1)=out2(i,1)+arco(i,j)*outvar(1,nd)
              else if (inode_elem==2) then !elem
                if(iep(i)<=0) stop 'iep(i)<=0'
                out2(i,1)=outvar(1,iep(i))
              endif
            enddo !j
          enddo !i
          write(18,'(e16.8,20000(1x,e14.6))')timeout(irec)/86400,(out2(i,1),i=1,nxy)
        else !3D
!         Do interpolation
          do i=1,nxy
            etal=0; dep(i)=0; idry=0
            do j=1,3
              nd=node3(i,j)
              if(eta2(nd)+dp(nd)<h0) idry=1
              etal=etal+arco(i,j)*eta2(nd)
              dep(i)=dep(i)+arco(i,j)*dp(nd)
!             Debug
!              write(11,*)i,j,nd,dp(nd),arco(i,j)
            enddo !j
            if(idry==1) then
              out3(i)=-99
!              write(65,*)'Dry'
            else !element wet
              !Compute z-coordinates
              if(ivcor==1) then !localized
                !Strictly speaking for elem-based vars, we need to use i34 nodes
                do j=1,3
                  nd=node3(i,j)
                  do k=kbp(nd)+1,nvrt-1
                    ztmp2(k,j)=(eta2(nd)+dp(nd))*sigma_lcl(k,nd)+eta2(nd)
                  enddo !k
                  ztmp2(kbp(nd),j)=-dp(nd) !to avoid underflow
                  ztmp2(nvrt,j)=eta2(nd) !to avoid underflow
                enddo !j

                ztmp=0
                kbpl=minval(kbp(node3(i,1:3)))
                do k=kbpl,nvrt
                  do j=1,3
                    nd=node3(i,j)
                    ztmp(k)=ztmp(k)+arco(i,j)*ztmp2(max(k,kbp(nd)),j)
                  enddo !j
                enddo !k
              else if(ivcor==2) then !SZ
                call zcor_SZ_single(dep(i),etal,h0,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma, &
     &ztmp(:),idry2,kbpl)
              endif

!             Horizontal interpolation
              if(inode_elem==1) then !node based
                do k=kbpl,nvrt
                  do j=1,3
                    nd=node3(i,j)
                    kin=max(k,kbp00(nd))
                    out2(i,k)=out2(i,k)+arco(i,j)*outvar(kin,nd) !out(i,j,kin,m)
                  enddo !j
                enddo !k
              endif !inode_elem

              if(abs(z00(i)-1.e10)<0.1) then !depth average
                total_dp=ztmp(nvrt)-ztmp(kbpl)
                if(total_dp<=0) then
                  write(*,*)'depth<=0:',total_dp,i
                  stop
                endif
                sum1=0
                do k=kbpl,nvrt-1
                  if(inode_elem==1) then !node
                    sum1=sum1+(ztmp(k+1)-ztmp(k))*(out2(i,k)+out2(i,k+1))/2
                  else if(inode_elem==2) then !elem
                    sum1=sum1+(ztmp(k+1)-ztmp(k))*outvar(k+1,iep(i))
                  endif
                enddo !k
                out3(i)=sum1/total_dp
              else !Interplate in vertical
                if(ifs==0) then !relative to MSL
                  z2=z00(i)
                else
                  z2=ztmp(nvrt)-z00(i)
                endif
                if(z2>=ztmp(nvrt)) then !above F.S.
                  k0=nvrt-1; rat=1
                else if(z2<=ztmp(kbpl)) then !below bottom; extrapolate
                  k0=kbpl; rat=0
                else !above bottom; cannot be above F.S.
                  k0=0
                  do k=kbpl,nvrt-1
                    if(z2>=ztmp(k).and.z2<=ztmp(k+1)) then
                      k0=k
                      rat=(z2-ztmp(k))/(ztmp(k+1)-ztmp(k))
                      exit
                    endif
                  enddo !k
                endif !ztmp

                if(k0==0) then
                  write(*,*)'read_output7b_xyz: failed to find a vertical level:',irec_real,i,ifs,z2,ztmp(:)
!'
                  stop
                else
                    if(inode_elem==1) then !node
                      out3(i)=out2(i,k0)*(1-rat)+out2(i,k0+1)*rat
                    else if(inode_elem==2) then !elem
                      if(iep(i)<=0) stop 'iep(i)<=0(2)'
                      out3(i)=outvar(k0+1,iep(i)) !FV
                    endif
                endif
              endif !depth average or not
            endif !dry/wet
          enddo !i=1,nxy
          write(18,'(e16.8,20000(1x,f14.6))')timeout(irec)/86400,(out3(i),i=1,nxy)

        endif !i23d
      enddo !irec
      iret=nf90_close(ncid)
      iret=nf90_close(ncid4)

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      enddo !iday

!     Output local depths info
      do i=1,nxy
        write(20,*)i,dep(i)
      enddo !i

      print*, 'Finished!'

      end program read_out

max_stars_repo_name: rustychris/schism

max_issues_repo_name: rustychris/schism

max_forks_repo_name: rustychris/schism


text:
      SUBROUTINE SUPLT (IZ,IY,X,U,GPLST,PEN,DEFORM)
C
C     TO CREATE A SET OF UNIQUE LINES TO BE PLOTTED.  IN ADDITION, TO
C     AVOID SKIPPING ALL OVER THE PLOT FOR EACH LINE.
C
C     INPUT (SIMULAR TO -GPCT-)
C       NGRID - NUMBER OF INTERNAL GRID POINTS.
C       IY    - 1 THRU NGRID - POINTERS TO FIRST CONNECTION OF THE
C                              INTERNAL GRID MATCHING THIS INDEX.
C                              IF THE GRID HAS NO ENTRIES IZ(I+1) WILL
C                              HAVE THE SAME POINTER VALUE.
C             - NGRID+1      - POINTER TO END-OF-RECORD.
C       IZ    - CONNECTING INTERNAL GRIDS.  POINTER FOR NEXT GRID
C                              DETERMINES LAST ENTRY.  ENTRIES PUSHED
C                              DOWN AND  -1  ADDED AT END AS EACH ENTRY
C                              IS USED.
C
C       NTAB  = TOTAL ENTRY COUNTER IN GPCT
C       ID1   = START OF CURRENT -LINE-
C       ID2   = END   OF CURRENT -LINE-
C       ID3   = START OF LAST -LINE-
C       ID4   = END   OF LAST -LINE-
C
      INTEGER         IZ(1),IY(1),PEN,DEFORM,GPLST(1),NM(5),M(71),
     1                M1(17),M2(17),M3(11),M4(17),M5(9),ERR(4),LM(5),
     2                LIMT(2)
      REAL            U(2,1),X(3,1)
      COMMON /SYSTEM/ SYSBUF,IOUT
      COMMON /BLANK / NGRID,SKP1(19),MERR
      EQUIVALENCE     (M1(1),M( 1)), (M2(1),M(18)), (M3(1),M(35)),
     1                (M4(1),M(46)), (M5(1),M(63))
      DATA    NM    / 17,17,11,17, 9 /,
     1        LM    /  1,18,35,46,63 /,
     2        M1    / 4H(35X, 4H,26H, 4HSUPL, 4HT RE, 4HJECT, 4HED P,
     3                4HLOT., 4H PIV, 4HOT,I, 4H8,26, 4HH IS, 4H ZER,
     4                4HO OR, 4H SAM, 4HE AS, 4H ENT, 4HRY.)/,
     5        M2    / 4H(35X, 4H,54H, 4HSUPL, 4HT RE, 4HJECT, 4HED P,
     6                4HLOT., 4H NEG, 4HATIV, 4HE NU, 4HMBER, 4H ENT,
     7                4HRIES, 4H - N, 4H3,N4, 4H =,2, 4HI10)/,
     8        M3    / 4H(35X, 4H,31H, 4HUNEX, 4HPECT, 4HED E, 4HOF I,
     9                4HN SU, 4HPLT , 4H- PI, 4HVOT,, 4HI10)/,
     O        M4    / 4H(35X, 4H,11H, 4HSUPL, 4HT-EN, 4HTRY,, 4HI10,,
     1                4H22H , 4HFOR , 4HPIVO, 4HT NO, 4HT FO, 4HUND ,
     2                4H(,2I, 4H6,8H, 4H) RA, 4HNGE., 4H)   /,
     3        M5    / 4H(35X, 4H,24H, 4HNO E, 4HLEME, 4HNTS , 4HIN T,
     4                4HHIS , 4HSET., 4H)   /
C
      LINESP = 0
C
C     LOCATE FIRST PIVOT (ID1) WITH ODD NUMBER OF ENTRIES
C
      ID2 = 0
      ID1 = 0
      DO 10 I = 1,NGRID
      LL  = IY(I+1) - IY(I)
      IF (LL .EQ. 0) GO TO 10
C
C     IN CASE AN ODD NO. ENTRIES ISN'T FOUND THE DEFAULT IS FIRST PIVOT
C
      IF (ID2 .EQ. 0) ID2 = I
      IF (MOD(LL,2) .EQ. 0) GO TO 10
      ID1 = I
      I1  = IY(I)
      ID2 = IZ(I1)
      GO TO 20
   10 CONTINUE
C
   20 NTAB = IY(NGRID+1) - IY(1)
C
C     NO ELEMENTS IN THE SET
C
      IF (NTAB .EQ. 0) GO TO 440
C
C     SEE IF ANY ODD ENTRIES FOUND
C
      LIMT(1) = 0
      LIMT(2) = NGRID + 1
      IF (ID1 .NE. 0) GO TO 140
      ID1 = ID2
      I1  = IY(ID1)
      ID2 = IZ(I1)
      GO TO 140
C
C     START OF LOOP AFTER FIRST -LINE-
C
   30 IF (N4) 410,260,40
C
C     LAST END HAS ENTRY TO CONTINUE FROM
C
   40 ID1 = ID4
C
C     INTERNAL SEARCH FOR FIRST GPCT ENTRY ABOVE AND BELOW THE PIVOT
C     VALUE FOR THE PIVOT           *** M I N / M A X ***
C
   50 I1 = IY(ID1)
      IL =-100000
      IH = 100000
      J1 = IY(ID1+1) - 1
C
      DO 80 I = I1,J1
      IF (IZ(I)    ) 100,400,60
   60 IF (IZ(I)-ID1)  70,400,90
   70 IL = IZ(I)
   80 CONTINUE
      GO TO 100
C
   90 IH = IZ(I)
C
C     DETERMINE  WHICH IS CLOSER TO PIVOT
C
  100 I = ID1 - IL
      J = IH  - ID1
      IF (J-I) 130,110,120
C
C     EQUAL DISTANT, GO TO SAME DIRECTION AS BEFORE
C
  110 IF (ID4-ID3) 120,400,130
C
C     ID2 IS LESSOR ID
C
  120 ID2 = IL
      GO TO 140
C
C     ID2 IS GREATER ID
C
  130 ID2 = IH
C
C     OUTPUT THE LINE -
C     NOTE THAT ID4 MAY BE RESET AT 320 SO DONT TAKE SHORTCUTS
C
  140 CONTINUE
      I  = IABS(GPLST(ID1))
      J  = IABS(GPLST(ID2))
      IF (DEFORM .NE. 0) GO TO 160
      X1 = X(2,I)
      Y1 = X(3,I)
      X2 = X(2,J)
      Y2 = X(3,J)
      GO TO 170
  160 X1 = U(1,I)
      Y1 = U(2,I)
      X2 = U(1,J)
      Y2 = U(2,J)
  170 CONTINUE
      CALL LINE (X1,Y1,X2,Y2,PEN,0)
      LINESP = LINESP + 1
C
C     REMOVE ENTRIES FROM CORE, LEFT SHIFT AS NEEDED AND PUT -1 AT THE
C     END OF THE TABLE.  PLACE THE NUMBER OF ENTRIES LEFT IN N3 AND N4.
C     DECREMENT THE TOTAL NUMBER OF ENTRIES BY 2. SET ID3 AND ID4.
C
      IF (NTAB .LE. 2) GO TO 460
      KK = 0
      J1 = I1
      J2 = IY(ID1+1) - 1
      LL = ID2
  180 IPAR = J1
      IL = 0
C
      DO 230 I = J1,J2
      IF (IZ(I)-LL) 210,200,190
  190 IZ(IPAR) = IZ(I)
      GO TO 220
C
C     COMPONENT TO BE ELIMINATED HAS BEEN FOUND
C
  200 IL = 1
      GO TO 230
  210 IF (IZ(I)) 240,400,220
  220 IPAR = IPAR + 1
  230 CONTINUE
C
  240 IZ(IPAR) = -LL
      IF (IL .EQ. 0) GO TO 430
      IF (KK .NE. 0) GO TO 250
      KK = 2
      ID3= ID1
      N3 = IPAR - J1
      LL = ID2
      J1 = IY(LL  )
      J2 = IY(LL+1) - 1
      LL = ID1
      GO TO 180
  250 NTAB= NTAB - 2
      ID4 = ID2
      N4  = IPAR - J1
      GO TO 30
C
C     CASE ID4 HAS NO MORE ENTRIES. CHECK IF ID3 CAN BE PIVOT
C
  260 IF (N3) 410,280,270
C
C     NONZERO - ID3 IS TO BE ID1
C
  270 ID1 = ID3
      GO TO 50
C
C     ID3 AND ID4 ARE NULL.  GO TO CLOSEST END OF TABLE FROM ID4
C
  280 I = NGRID - ID4
      J = 1
      IF (I .GT. ID4) J = -1
      L = (J+2)/2 + 1
      LIM = LIMT(L)
      LEN = ID4
C
      ASSIGN 310 TO IRET1
      KK = ID4
  290 KK = KK + J
      IF (KK .EQ. LIM) GO TO IRET1, (310,420)
      IPAR = 2
      ASSIGN 300 TO IRET
      GO TO 370
C
C     CHECK IF ANY ENTRIES FOUND
C
  300 IF (IPAR .EQ. 0) GO TO 290
C
C     ENTRY FOUND
C
      LEN = KK + J
      ID4 = KK
      N4  = IPAR
      GO TO 320
C
C     THAT END OF TABLE FAILED - TRY OTHER END
C
  310 J   = -J
      LIMT(L) = LEN
      L   = (J+2)/2 + 1
      LIM = LIMT(L)
      KK  = ID4
      ASSIGN 420 TO IRET1
      GO TO 290
C
C     AN ENTRY WAS FOUND - CHECK FOR ODD NUMBER OF ENTRIES FOR PIVOT
C
  320 IF (MOD(IPAR,2) .EQ. 1) GO TO 360
C
C     NOT AN ODD NUMBER OF ENTRIES FOR ID4.  CHECK GPCT ENTRIES
C     FOR ONLY ONE ENTRY.
C
      ASSIGN 340 TO IRET
C
      IH = J2
      IL = J1 - 1
  330 IL = IL + 1
      KK = IZ(IL)
      IF (KK .LE. 0) GO TO 40
      IPAR = 1
      GO TO 370
  340 CONTINUE
      IF (IPAR .EQ.  1) GO TO 360
      IF (IL   .LT. IH) GO TO 330
C
C     PIVOT NOW DETERMINED
C
      GO TO 40
  360 ID1 = KK
      GO TO 50
C
C
C     INTERNAL ROUTINE TO DETERMINE NUMBER OF ENTRIES FOR PIVOT
C
C     INPUT
C        IPAR = 1 -- 0,1 OR MORE THAN 1 ENTRY RETURN
C             = 2 -- ACTUAL NUMBER OF ENTRIES RETURN
C        KK   = ID OF PIVOT
C
C     OUTPUT
C        IPAR = DESIRED NUMBER OF ENTRIES
C        KK   = SAME AS INPUT
C        J1   = POINTER TO 1ST LOCATION
C        J2   = NOT NECESSARILY LAST LOCATION (I.E. IPAR INPUT AS 1)
C
  370 J1 = IY(KK  )
      J2 = IY(KK+1) - 1
      IF (IPAR .EQ. 1) J2 = MIN0(J1+2,J2)
      IPAR = 0
      IF (J2-J1 .LT. 0) GO TO 390
C
      DO 380 I = J1,J2
      IF (IZ(I)) 390,400,380
  380 IPAR = IPAR + 1
  390 GO TO IRET, (300,340)
C
C     ERROR MESSAGES
C
  400 ERR(1) = 1
      ERR(2) = ID1
      K = 1
      GO TO 450
  410 ERR(1) = 2
      ERR(2) = N3
      ERR(3) = N4
      K = 2
      GO TO 450
  420 ERR(1) = 1
      ERR(2) = ID4
      K = 3
      GO TO 450
  430 ERR(1) = 3
      ERR(2) = LL
      ERR(3) = J1
      ERR(4) = J2
      K = 4
      GO TO 450
  440 ERR(1) = 0
      K = 5
C
  450 I = LM(K)
      CALL WRTPRT (MERR,ERR,M(I),NM(K))
      IF (K .EQ. 5) GO TO 530
C
C     CONVERT TABLE TO ORIGINAL VALUES UNLESS THIS IS THE LAST CALL
C
  460 IL = NGRID  + 1
      I  = IY(IL) - 1
      IF (DEFORM .NE. 0) GO TO 530
      DO 480 J = 1,I
  480 IZ(J) = IABS(IZ(J))
C
      DO 520 J1 = 1,NGRID
      IF (IY(J1) .EQ. IY(J1+1)) GO TO 520
      I = IY(J1)
      L = I
      N = IY(J1+1) - 1
      IF (I+1 .GT. N) GO TO 520
C
C     SHUTTLE EXCHANGE
C     (NOTE FROM G.CHAN/UNISYS  10/1990
C     THERE ARE MORE THAN JUST A SHUTTLE SORTING HERE. REPLACING THE
C     SHUTTLE EXCHANGE METHOD BY SORT ROUTINE, WHICH USES A MUCH FASTER
C     TECHNEQUE, DOES NOT WORK HERE)
C
  490 IF (IZ(I) .LE. IZ(I+1)) GO TO 510
      K = IZ(I+1)
      IZ(I+1) = IZ(I)
      IZ(I  ) = K
      J = I
  500 IF (J .EQ. L) GO TO 510
      IF (IZ(J) .GE. IZ(J-1)) GO TO 510
      K = IZ(J)
      IZ(J  ) = IZ(J-1)
      IZ(J-1) = K
      J = J - 1
      GO TO 500
  510 IF (I .GE. N-1) GO TO 520
      I = I + 1
      GO TO 490
C
  520 CONTINUE
C
C     A NONSTANDARD RETURN COULD BE ADDED HERE.  BAD PLOT RESULTS IF
C     THIS ROUTINE FAILS.  THE FRAME WILL BE PRESENT HOWEVER
C
  530 CONTINUE
      RETURN
      END


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_types_regionCharges

IMPLICIT NONE

PRIVATE

   TYPE t_regionCharges

      REAL,    ALLOCATABLE :: sqal(:,:,:)
      REAL,    ALLOCATABLE :: ener(:,:,:)
      REAL,    ALLOCATABLE :: sqlo(:,:,:)
      REAL,    ALLOCATABLE :: enerlo(:,:,:)
      REAL,    ALLOCATABLE :: svac(:,:)
      REAL,    ALLOCATABLE :: pvac(:,:)

      CONTAINS
         PROCEDURE,PASS :: init => regionCharges_init
         PROCEDURE      :: sumBandsVac
   END TYPE t_regionCharges

PUBLIC t_regionCharges

CONTAINS

SUBROUTINE regionCharges_init(thisRegCharges,input,atoms)

   USE m_types_input
   USE m_types_atoms

   IMPLICIT NONE

   CLASS(t_regionCharges), INTENT(INOUT) :: thisRegCharges
   TYPE(t_input),          INTENT(IN)    :: input
   TYPE(t_atoms),          INTENT(IN)    :: atoms

   ALLOCATE(thisRegCharges%sqal(0:3,atoms%ntype,input%jspins))
   ALLOCATE(thisRegCharges%ener(0:3,atoms%ntype,input%jspins))

   ALLOCATE(thisRegCharges%sqlo(atoms%nlod,atoms%ntype,input%jspins))
   ALLOCATE(thisRegCharges%enerlo(atoms%nlod,atoms%ntype,input%jspins))
   ALLOCATE(thisRegCharges%svac(2,input%jspins))
   ALLOCATE(thisRegCharges%pvac(2,input%jspins))

   thisRegCharges%sqal = 0.0
   thisRegCharges%ener = 0.0
   thisRegCharges%sqlo = 0.0
   thisRegCharges%enerlo = 0.0
   thisRegCharges%svac = 0.0
   thisRegCharges%pvac = 0.0

END SUBROUTINE regionCharges_init

SUBROUTINE sumBandsVac(thisRegCharges,vacuum,dos,noccbd,ikpt,jsp_start,jsp_end,eig,we)

  USE m_types_vacuum
  USE m_types_dos

   USE m_types_dos

   IMPLICIT NONE

   CLASS(t_regionCharges), INTENT(INOUT) :: thisRegCharges
   TYPE(t_vacuum),         INTENT(IN)    :: vacuum
   TYPE(t_dos),            INTENT(IN)    :: dos
   INTEGER,                INTENT(IN)    :: noccbd
   INTEGER,                INTENT(IN)    :: ikpt
   INTEGER,                INTENT(IN)    :: jsp_start, jsp_end
   REAL,                   INTENT(IN)    :: eig(noccbd)
   REAL,                   INTENT(IN)    :: we(noccbd)

   INTEGER                               :: ispin, ivac

   ! perform Brillouin zone integration and summation over the bands in order to determine the vacuum energy parameters.
   DO ispin = jsp_start, jsp_end
      DO ivac = 1,vacuum%nvac
         thisRegCharges%pvac(ivac,ispin) = thisRegCharges%pvac(ivac,ispin)+&
            dot_product(eig(:noccbd)*dos%qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd))
         thisRegCharges%svac(ivac,ispin)=thisRegCharges%svac(ivac,ispin)+&
            dot_product(dos%qvac(:noccbd,ivac,ikpt,ispin),we(:noccbd))
      END DO
   END DO

END SUBROUTINE sumBandsVac

END MODULE m_types_regionCharges

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
      SUBROUTINE FRMLTD (IFILE,DZ,DY,ZM)
C
C     FEER MATRIX TRANSPOSE MULTIPLY  (DOUBLE PREC)
C               T
C     DY = IFILE * DZ        WHERE DZ IS A VECTOR ALREADY IN CORE
C                            IFILE IS A GINO MATIRX FILE
C
C     LAST REVISED  11/91, BY C.CHAN/UNISYS
C     ADDITION OF A NEW TRANSPOSE MULTIPLY METHOD WHICH IS MORE
C     EFFECIENT, AND IS ALREADY GOOD FOR VECTORIZATION
C
CDB   LOGICAL          DEBUG
      DOUBLE PRECISION DZ(1)    ,DY(1)    ,ZM(1)    ,DP     ,DSUM
      DIMENSION        IFILE(7) ,IDP(2)   ,NAM(2)
      COMMON  /UNPAKX/ ITYP     ,IP       ,NP       ,INCR
      COMMON  /SYSTEM/ IBUF     ,NOUT
      COMMON  /FEERXX/ DUM18(18),NZM
      COMMON  /ZZZZZZ/ IZ(1)
      EQUIVALENCE      (DP,IDP(1))
      DATA     NAM   / 4HFRML   ,4HTD    /
CDB   DATA     DEBUG , ITER     ,MAX     / .FALSE.  ,0      ,4       /
C
CDB   IF (.NOT.DEBUG) GO TO 20
C     ITER = ITER + 1
C     IF (ITER .GT. MAX) DEBUG = .FALSE.
C     IF (DEBUG) WRITE (NOUT,10) NZM,IFILE(5)
C  10 FORMAT ('  .... IN FRMLTD DEBUG.   NZM,IFILE(5) =',2I8)
C  20 CONTINUE
      N    = IFILE(2)
      IFL  = IFILE(1)
      IF (IFILE(7) .LT. 0) IFL = -IFILE(7)
      CALL REWIND (IFL)
      CALL SKPREC (IFL,1)
      IF (IFILE(7) .LT. 0) GO TO 50
      ITYP = IFILE(5)
C
C     NASTRAN ORIGIANL METHOD
C
      INCR  = 1
      DO 40 I = 1,N
      DY(I) = 0.0D+0
      IP    = 0
      CALL UNPACK (*40,IFL,ZM(1))
      DSUM  = 0.0D+0
      II    = 0
      DO 30 J = IP,NP
      II    = II + 1
   30 DSUM  = DSUM + ZM(II)*DZ(J)
      DY(I) = DSUM
   40 CONTINUE
      GO TO 200
C
C     NEW METHOD, READ ONLY AND NO UNPACK
C
C     UNLIKE FRMLTX, IFL WAS UNPACKED FORWARD BY UNPSCR
C
   50 NREC = 0
      NWDS = IFILE(5)
CDB   N20  = N - 20
C     IF (DEBUG) WRITE (NOUT,60) NWDS,NZM
C  60 FORMAT ('  /@60   NWDS,NZM =',2I8)
      LL2  = 0
      NEXT = 1
      DO 140 I = 1,N
      IF (NEXT .LT. LL2) GO TO 100
      NREC = NREC + 1
CDB   IF (DEBUG) WRITE (NOUT,70) NREC,I
C  70 FORMAT ('  ...READING RECORD',I5,'.  I =',I7)
      CALL READ (*150,*80,IFL,ZM,NZM,1,LL)
      CALL MESAGE (-8,0,NAM)
   80 LL2  = LL/NWDS
CDB   IF (DEBUG) WRITE (NOUT,90) LL,NREC,LL2
C  90 FORMAT (1X,I10,' WORDS READ FROM RECORD NO.',I5,'   LL2 =',I10)
      NEXT = 1
  100 DP   = ZM(NEXT)
      II   = IDP(1)
      JJ   = IDP(2)
CDB   IF (DEBUG .AND. (I.LT.20 .OR. I.GT.N20)) WRITE (NOUT,110) I,II,JJ,
C     1                                                         NEXT
C 110 FORMAT ('   @110  I,II,JJ,NEXT =',4I8)
      IF (II .EQ. JJ) GO TO 130
      DSUM = 0.0D+0
      LL   = NEXT
      DO 120 J = II,JJ
      LL   = LL + 1
  120 DSUM = DSUM + ZM(LL)*DZ(J)
      DY(I)= DSUM
      GO TO 140
  130 DY(I)= ZM(NEXT+1)*DZ(II)
  140 NEXT = NEXT + JJ - II + 2
      GO TO 200
C
  150 J = IFILE(4)/10
      WRITE  (NOUT,160) NREC,I,N,J
  160 FORMAT ('0*** TRY TO READ RECORD',I5,'.   I,N,IFILE(4) =',2I7,I5)
      CALL MESAGE (-2,IFL,NAM)
C
  200 RETURN
      END


text:
*DECK PCHNGS
      SUBROUTINE PCHNGS (II, XVAL, IPLACE, SX, IX, IRCX)
C***BEGIN PROLOGUE  PCHNGS
C***SUBSIDIARY
C***PURPOSE  Subsidiary to SPLP
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (PCHNGS-S, DPCHNG-D)
C***AUTHOR  Hanson, R. J., (SNLA)
C           Wisniewski, J. A., (SNLA)
C***DESCRIPTION
C
C     PCHNGS LIMITS THE TYPE OF STORAGE TO A SEQUENTIAL SCHEME.
C     SPARSE MATRIX ELEMENT ALTERATION SUBROUTINE.
C
C     SUBROUTINE PCHNGS() CHANGES ELEMENT II IN VECTOR +/- IRCX TO THE
C     VALUE XVAL.
C
C            II THE ABSOLUTE VALUE OF THIS INTEGER IS THE SUBSCRIPT FOR
C               THE ELEMENT TO BE CHANGED.
C          XVAL NEW VALUE OF THE MATRIX ELEMENT BEING CHANGED.
C     IPLACE POINTER INFORMATION WHICH IS MAINTAINED BY THE PACKAGE.
C   SX(*),IX(*) THE WORK ARRAYS WHICH ARE USED TO STORE THE SPARSE
C               MATRIX. THESE ARRAYS ARE AUTOMATICALLY MAINTAINED BY THE
C               PACKAGE FOR THE USER.
C          IRCX POINTS TO THE VECTOR OF THE MATRIX BEING UPDATED.
C               A NEGATIVE VALUE OF IRCX INDICATES THAT ROW -IRCX IS
C               BEING UPDATED.  A POSITIVE VALUE OF IRCX INDICATES THAT
C               COLUMN IRCX IS BEING UPDATED.  A ZERO VALUE OF IRCX IS
C               AN ERROR.
C
C     SINCE DATA ITEMS ARE KEPT SORTED IN THE SEQUENTIAL DATA STRUCTURE,
C     CHANGING A MATRIX ELEMENT CAN REQUIRE THE MOVEMENT OF ALL THE DATA
C     ITEMS IN THE MATRIX. FOR THIS REASON, IT IS SUGGESTED THAT DATA
C     ITEMS BE ADDED A COL. AT A TIME, IN ASCENDING COL. SEQUENCE.
C     FURTHERMORE, SINCE DELETING ITEMS FROM THE DATA STRUCTURE MAY ALSO
C     REQUIRE MOVING LARGE AMOUNTS OF DATA, ZERO ELEMENTS ARE EXPLICITLY
C     STORED IN THE MATRIX.
C
C     THIS SUBROUTINE IS A MODIFICATION OF THE SUBROUTINE LCHNGS,
C     SANDIA LABS. REPT. SAND78-0785.
C     MODIFICATIONS BY K.L. HIEBERT AND R.J. HANSON
C     REVISED 811130-1000
C     REVISED YYMMDD-HHMM
C
C***SEE ALSO  SPLP
C***ROUTINES CALLED  IPLOC, PRWPGE, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   811215  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900328  Added TYPE section.  (WRB)
C   910403  Updated AUTHOR and DESCRIPTION sections.  (WRB)
C***END PROLOGUE  PCHNGS
      DIMENSION IX(*)
      INTEGER IPLOC
      REAL SX(*),XVAL,ZERO,ONE,SXLAST,SXVAL
      SAVE ZERO, ONE
      DATA ZERO,ONE /0.E0,1.E0/
C***FIRST EXECUTABLE STATEMENT  PCHNGS
      IOPT=1
C
C     DETERMINE NULL-CASES..
      IF(II.EQ.0) RETURN
C
C     CHECK VALIDITY OF ROW/COL. INDEX.
C
      IF (.NOT.(IRCX.EQ.0)) GO TO 20002
      NERR=55
      CALL XERMSG ('SLATEC', 'PCHNGS', 'IRCX=0.', NERR, IOPT)
20002 LMX = IX(1)
C
C     LMX IS THE LENGTH OF THE IN-MEMORY STORAGE AREA.
C
      IF (.NOT.(IRCX.LT.0)) GO TO 20005
C
C     CHECK SUBSCRIPTS OF THE ROW. THE ROW NUMBER MUST BE .LE. M AND
C     THE INDEX MUST BE .LE. N.
C
      IF (.NOT.(IX(2).LT.-IRCX .OR. IX(3).LT.ABS(II))) GO TO 20008
      NERR=55
      CALL XERMSG ('SLATEC', 'PCHNGS',
     +   'SUBSCRIPTS FOR ARRAY ELEMENT TO BE ACCESSED WERE OUT OF ' //
     +   'BOUNDS.', NERR, IOPT)
20008 GO TO 20006
C
C     CHECK SUBSCRIPTS OF THE COLUMN. THE COL. NUMBER MUST BE .LE. N AND
C     THE INDEX MUST BE .LE. M.
C
20005 IF (.NOT.(IX(3).LT.IRCX .OR. IX(2).LT.ABS(II))) GO TO 20011
      NERR=55
      CALL XERMSG ('SLATEC', 'PCHNGS',
     +   'SUBSCRIPTS FOR ARRAY ELEMENT TO BE ACCESSED WERE OUT OF ' //
     +   'BOUNDS.', NERR, IOPT)
20011 CONTINUE
C
C     SET I TO BE THE ELEMENT OF ROW/COLUMN J TO BE CHANGED.
C
20006 IF (.NOT.(IRCX.GT.0)) GO TO 20014
      I = ABS(II)
      J = ABS(IRCX)
      GO TO 20015
20014 I = ABS(IRCX)
      J = ABS(II)
C
C     THE INTEGER LL POINTS TO THE START OF THE MATRIX ELEMENT DATA.
C
20015 LL=IX(3)+4
      II = ABS(II)
      LPG = LMX - LL
C
C     SET IPLACE TO START OUR SCAN FOR THE ELEMENT AT THE BEGINNING
C     OF THE VECTOR.
C
      IF (.NOT.(J.EQ.1)) GO TO 20017
      IPLACE=LL+1
      GO TO 20018
20017 IPLACE=IX(J+3)+1
C
C     IEND POINTS TO THE LAST ELEMENT OF THE VECTOR TO BE SCANNED.
C
20018 IEND = IX(J+4)
C
C     SCAN THROUGH SEVERAL PAGES, IF NECESSARY, TO FIND MATRIX ELEMENT.
C
      IPL = IPLOC(IPLACE,SX,IX)
      NP = ABS(IX(LMX-1))
      GO TO 20021
20020 IF (ILAST.EQ.IEND) GO TO 20022
C
C     THE VIRTUAL END OF DATA FOR THIS PAGE IS ILAST.
C
20021 ILAST = MIN(IEND,NP*LPG+LL-2)
C
C     THE RELATIVE END OF DATA FOR THIS PAGE IS IL.
C     SEARCH FOR A MATRIX VALUE WITH AN INDEX .GE. I ON THE PRESENT
C     PAGE.
C
      IL = IPLOC(ILAST,SX,IX)
      IL = MIN(IL,LMX-2)
20023 IF (.NOT.(.NOT.(IPL.GE.IL .OR. IX(IPL).GE.I))) GO TO 20024
      IPL=IPL+1
      GO TO 20023
C
C     SET IPLACE AND STORE DATA ITEM IF FOUND.
C
20024 IF (.NOT.(IX(IPL).EQ.I .AND. IPL.LE.IL)) GO TO 20025
      SX(IPL) = XVAL
      SX(LMX) = ONE
      RETURN
C
C     EXIT FROM LOOP IF ITEM WAS FOUND.
C
20025 IF(IX(IPL).GT.I .AND. IPL.LE.IL) ILAST = IEND
      IF (.NOT.(ILAST.NE.IEND)) GO TO 20028
      IPL = LL + 1
      NP = NP + 1
20028 GO TO 20020
C
C     INSERT NEW DATA ITEM INTO LOCATION AT IPLACE(IPL).
C
20022 IF (.NOT.(IPL.GT.IL.OR.(IPL.EQ.IL.AND.I.GT.IX(IPL)))) GO TO 20031
      IPL = IL + 1
      IF(IPL.EQ.LMX-1) IPL = IPL + 2
20031 IPLACE = (NP-1)*LPG + IPL
C
C     GO TO A NEW PAGE, IF NECESSARY, TO INSERT THE ITEM.
C
      IF (.NOT.(IPL.LE.LMX .OR. IX(LMX-1).GE.0)) GO TO 20034
      IPL=IPLOC(IPLACE,SX,IX)
20034 IEND = IX(LL)
      NP = ABS(IX(LMX-1))
      SXVAL = XVAL
C
C     LOOP THROUGH ALL SUBSEQUENT PAGES OF THE MATRIX MOVING DATA DOWN.
C     THIS IS NECESSARY TO MAKE ROOM FOR THE NEW MATRIX ELEMENT AND
C     KEEP THE ENTRIES SORTED.
C
      GO TO 20038
20037 IF (IX(LMX-1).LE.0) GO TO 20039
20038 ILAST = MIN(IEND,NP*LPG+LL-2)
      IL = IPLOC(ILAST,SX,IX)
      IL = MIN(IL,LMX-2)
      SXLAST = SX(IL)
      IXLAST = IX(IL)
      ISTART = IPL + 1
      IF (.NOT.(ISTART.LE.IL)) GO TO 20040
      K = ISTART + IL
      DO 50 JJ=ISTART,IL
      SX(K-JJ) = SX(K-JJ-1)
      IX(K-JJ) = IX(K-JJ-1)
50    CONTINUE
      SX(LMX) = ONE
20040 IF (.NOT.(IPL.LE.LMX)) GO TO 20043
      SX(IPL) = SXVAL
      IX(IPL) = I
      SXVAL = SXLAST
      I = IXLAST
      SX(LMX) = ONE
      IF (.NOT.(IX(LMX-1).GT.0)) GO TO 20046
      IPL = LL + 1
      NP = NP + 1
20046 CONTINUE
20043 GO TO 20037
20039 NP = ABS(IX(LMX-1))
C
C     DETERMINE IF A NEW PAGE IS TO BE CREATED FOR THE LAST ELEMENT
C     MOVED DOWN.
C
      IL = IL + 1
      IF (.NOT.(IL.EQ.LMX-1)) GO TO 20049
C
C     CREATE A NEW PAGE.
C
      IX(LMX-1) = NP
C
C     WRITE THE OLD PAGE.
C
      SX(LMX) = ZERO
      KEY = 2
      CALL PRWPGE(KEY,NP,LPG,SX,IX)
      SX(LMX) = ONE
C
C     STORE LAST ELEMENT MOVED DOWN IN A NEW PAGE.
C
      IPL = LL + 1
      NP = NP + 1
      IX(LMX-1) = -NP
      SX(IPL) = SXVAL
      IX(IPL) = I
      GO TO 20050
C
C     LAST ELEMENT MOVED REMAINED ON THE OLD PAGE.
C
20049 IF (.NOT.(IPL.NE.IL)) GO TO 20052
      SX(IL) = SXVAL
      IX(IL) = I
      SX(LMX) = ONE
20052 CONTINUE
C
C     INCREMENT POINTERS TO LAST ELEMENT IN VECTORS J,J+1,... .
C
20050 JSTART = J + 4
      JJ=JSTART
      N20055=LL
      GO TO 20056
20055 JJ=JJ+1
20056 IF ((N20055-JJ).LT.0) GO TO 20057
      IX(JJ) = IX(JJ) + 1
      IF(MOD(IX(JJ)-LL,LPG).EQ.LPG-1) IX(JJ) = IX(JJ) + 2
      GO TO 20055
C
C     IPLACE POINTS TO THE INSERTED DATA ITEM.
C
20057 IPL=IPLOC(IPLACE,SX,IX)
      RETURN
      END


text:
       integer function clcalc(iz, npts, energy, fp, fpp)
c cromer-libermann calculation of anomalous scattering factors
c arguments:
c   iz      atomic number of element                       [in]
c   npts    number of elements in energy array             [in]
c   energy  array of energies at which to calculate f'/f'' [in]
c   fp      real part of anomalous scattering   (f')       [out]
c   fpp     imag part of anomalous scattering   (f'')      [out]
c
c notes:
c   1  energy array is in eV
c   2  this code is based on, and modified from the cowan-brennan
c      routines.  data statements were simplified and rearranged,
c      code was cleaned up to be more in keeping with f77 standard
c
c  matthew newville oct 1996
       implicit none
       integer  iz, npts
       integer  i, j, k, nxpts, maxpts
       integer  nparms(92,24), norb(92), nparmz(24), norbz
       parameter  (maxpts = 65536)
       double precision  binden(92,24), xnrg(92,24,11), xsc(92,24,11)
       double precision  benaz(24),    xnrgz(24,11),   xscz(24,11)
       double precision  relcor(92), kpcor(92), xnrdat(5), kev2ry
       double precision  energy(maxpts), fp(maxpts), fpp(maxpts)
       double precision  ener, f1, f2
       parameter (kev2ry = 0.02721d0)
       include "cldata.f"

       clcalc = 0
c  initialize the trivial parts of xnrg:
c  (non-trivial parts are given in the data statements above)
       if (iz .gt. 92) return
       nxpts = min(maxpts, npts)
       norbz = norb(iz)
       do 50 k = 1, 5
          do 30 j = 1, norbz
             xnrg(iz,j,k) = xnrdat(k)
 30       continue
 50    continue
       if (iz.gt.3) then
c look-up coefficients for this element
          do 100 i = 1, 24
             nparmz(i) = nparms(iz,i)
             benaz(i)  = binden(iz,i)/ kev2ry
             do 80 j = 1, 11
                xnrgz(i,j) = xnrg(iz,i,j)
                xscz(i,j)  = xsc(iz,i,j)
 80          continue
 100      continue
c calculate fp and fpp for each energy point
          do 200 i = 1, nxpts
             ener = energy(i) / 1000
             call cromer(iz,ener,nparmz,norbz,benaz,xnrgz,xscz,f1,f2)
             fp(i)  = f1 - relcor(iz) + kpcor(iz)
             fpp(i) = f2
 200      continue
       end if
       return
       end


text:
      subroutine time_diff(time_fluid, nseed, b, retard_frac,
     2     fact_term, time_particle)
!***********************************************************************
!  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.
C***********************************************************************
CD1
CD1  PURPOSE
CD1
CD1  This subroutine is used in particle tracking to account for
CD1  diffusion.
CD1
C***********************************************************************
CD2
CD2  REVISION HISTORY
CD2
CD2 $Log:   /pvcs.config/fehm90/src/time_diff.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:44:22   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:21:32   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:28:46   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:28:20   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:13:00   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:47:46 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.2   Thu Jan 11 12:00:40 1996   hend
CD2 Added Prolog
CD2
CD2    Rev 1.1   03/16/95 09:48:50   llt
CD2 added PVCS log history
CD2
CD2    Rev 1.0   03/16/95 09:00:44   robinson
CD2 Initial revision.
CD2
C***********************************************************************
CD3
CD3  REQUIREMENTS TRACEABILITY
CD3
CD3  2.3.5 Cell-based particle-tracking module
CD3
C***********************************************************************
CD4
CD4  SPECIAL COMMENTS AND REFERENCES
CD4
CD4 Requirements from SDN: 10086-RD-2.20-00
CD4   SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD4   FEHM Application Version 2.20
CD4
C***********************************************************************

      implicit none

      real time_fluid,time_particle,xt,ran_sp,inverf,fact_term
      real retard_frac,b
      integer nseed

      xt = ran_sp(nseed)
      time_particle = retard_frac + fact_term *
     2     time_fluid / (b**2*inverf(xt)**2)

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
C*26-NOV-90 11:35:29 MDS$SOURCE:[KKLARE]ZBLAS.FOR;5
! Subprograms based on "LINPACK User's Guide" by J.J. Dongarra et al.
! The Basic Linear Algebra Subprograms portion.
!
! root    prefix-suffix of name     Function
! -AMAX   IS- ID-     IC-     IZ-   index of max |x|  or  max  |x|+|y|
! -ASUM   S-  D-      SC-     DZ-   w := sum |x|  or  sum |x|+|y|
! -AXPY   S-  D-      C-      Z-    w := a*x+y
! -COPY   S-  D-      C-      Z-    y := x
! -DOT-   S-  D-  C-U C-C Z-U Z-C   w:=sum x*y, for -C  w:=conjg(x)*y
! -NRM2   S-  D-      SC-     DZ-   w := sqrt(sum |x|**2)
! -ROT    S-  D-      CS-     ZD-   apply rotation
! -ROTG   S-  D-      C-      Z-    setup Givens rotation
! -SCAL   S-  D-  CS- C-  ZD- Z-    x := ax
! -SWAP   S-  D-      C-      Z-    y :=: x
! -ROTG used by -CHEX and -CHUD.
! S-ROTG and CS-ROT used by C-SVDC etc.
!
! S- are single-precision.
! D- are double-precision.
! C- are single-precision complex.
! Z- are double-precision complex.
!
! N is the active number of elements.
! X and Y are vectors stepped in increments INCX and INCY.
! Dimension 1 is used for these arrays.
! LINPACK uses INCX=INCY=1 and special code for this case.
! Negative INCX or INCY can mean reversed order.
! A, B, C, and S are scalars.
!
!+ index of absolute-value maximum
      Integer Function IZAMAX(n,x,incx)
      Integer i,incx,j,n
      Real*8 xmax
      Complex*16 x(*)
      If (n.le.0) Then
         Izamax=0
      ElseIf (n.eq.1) Then
         Izamax=1
      Else
         Izamax=1
         xmax=abs(dreal(x(1)))+abs(dimag(x(1)))
         If (incx.ne.1) Then
            j=1
            Do i=2,n
               j=j+incx
               If (abs(dreal(x(j)))+abs(dimag(x(j))).gt.xmax) Then
                xmax=abs(dreal(x(j)))+abs(dimag(x(j)))
                Izamax=i
               EndIf
            EndDo
         Else
            Do i=2,n
               If (abs(dreal(x(i)))+abs(dimag(x(i))).gt.xmax) Then
                  xmax=abs(dreal(x(i)))+abs(dimag(x(i)))
                  Izamax=i
               EndIf
            EndDo
         EndIf
      EndIf
      End
!+ Euclidean norm of vector
      Real*8 Function DZNRM2(n,x,incx)
      Real*8 DZASUM
      Integer i,incx,n
      Real*8 cutlo,cuthi,test
      Complex*16 x(*)
! machine precision      smallest-exp      largest-exp
! Prime            S 23      -128         127
! Prime            D 47      -32896(-32768-128) 32639(32767-128)
! Vax            S 24      -128         127
! Vax            D 56      -128         127
! DEC KL10      S 27      -128         127
! DEC KL10      D 62      -128         127 //worst case//
! DEC KI10      D 54      -128         127
! CDC-6600      S 48      -975(-1023+48)      1060(+1022+48) (+1023=infinity)
! CDC-6600      D 96      -975         1060
! The smallest number is 2**(smallest)*.5.
! The largest number is 2**(largest)*(1-2/2**precision).
! The expression eps+1 .gt. 1 is true for eps=2**(1-precision)
! unrounded or eps=2**(-precision) when rounded.
! This does not work on Prime where S-accumulator is 31 bits.
! For un-squared number halve the exponent.
! So we scale for maximum above largest/N (/2N for complex)
! And for maximum below CUTLO**2=smallest/eps.
! It is assumed that underflows are efficiently handled as zero.
! Round up CUTLO, down CUTHI.
!primecs          data cutlo,cuthi/1.570e-16,1.304e+19/
!primecd          data cutlo,cuthi/5.406e-4945,4.560e+4912/
!kl-10cs          data cutlo,cuthi/6.281e-16,1.304e+19/
!kl-10cd          data cutlo,cuthi/1.164e-10,1.304e+19/
!vaxcs          data cutlo,cuthi/8.881e-16,1.304e+19/
!vaxcd          data cutlo,cuthi/2.058e-11,1.304e+19/
      data cutlo,cuthi/8.232e-11,1.304e+19/
      DZNRM2=0
      If (n.le.0) Return
      test=DZASUM(n,x,incx)
      If (test.le.0) Return
!CIF
      If (test.le.cuthi.and.test*(n+n).ge.cutlo) Then
!CEND      If (test.le.cuthi.and.test*n.ge.cutlo) Then
         Do j=1,n*incx,incx
            DZNRM2=DZNRM2+dreal(x(j))**2
     1                  +dimag(x(j))**2
         EndDo
         DZNRM2=sqrt(DZNRM2)
      Else
         Do i=1,n*incx,incx
            DZNRM2=DZNRM2+(dreal(x(i))/test)**2
     1                  +(dimag(x(i))/test)**2
         EndDo
         DZNRM2=test*sqrt(DZNRM2)
      Endif
      End
!+ absolute values summed
      Real*8 Function DZASUM(n,x,incx)
      Integer i,incx,m,n
      Complex*16 x(*)
      DZASUM=0
      If (n.le.0) Then
      ElseIf (incx.ne.1) Then
         Do i=1,n*incx,incx
            DZASUM=DZASUM+abs(dreal(x(i)))+abs(dimag(x(i)))
         EndDo
      Else
         m=mod(n,6)
         Do i=1,m
            DZASUM=DZASUM+abs(dreal(x(i)))+abs(dimag(x(i)))
         EndDo
         Do i=m+1,n,6
            DZASUM=DZASUM+abs(dreal(x(i)))+abs(dimag(x(i)))
     1+abs(dreal(x(i+1)))+abs(dimag(x(i+1)))
     2+abs(dreal(x(i+2)))+abs(dimag(x(i+2)))
     3+abs(dreal(x(i+3)))+abs(dimag(x(i+3)))
     4+abs(dreal(x(i+4)))+abs(dimag(x(i+4)))
     5+abs(dreal(x(i+5)))+abs(dimag(x(i+5)))
         EndDo
      EndIf
      End
!+ y=a*x+y  scale and add, backwards if negative increment
      Subroutine ZAXPY(n,a,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      Complex*16 x(1),y(1),a
      If (n.le.0.or.a.eq.0) Then
      ElseIf (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            y(k)=y(k)+x(j)*a
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,4)
         Do i=1,m
            y(i)=y(i)+x(i)*a
         EndDo
         Do i=m+1,n,4
            y(i)=y(i)+x(i)*a
            y(i+1)=y(i+1)+x(i+1)*a
            y(i+2)=y(i+2)+x(i+2)*a
            y(i+3)=y(i+3)+x(i+3)*a
         EndDo
      EndIf
      End
!+ copy array, backwards if negative increment
      Subroutine ZCOPY(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      Complex*16 x(1),y(1)
      If (n.le.0) Then
      ElseIf (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            y(k)=x(j)
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,7)
         Do i=1,m
            y(i)=x(i)
         EndDo
         Do i=m+1,n,7
            y(i)=x(i)
            y(i+1)=x(i+1)
            y(i+2)=x(i+2)
            y(i+3)=x(i+3)
            y(i+4)=x(i+4)
            y(i+5)=x(i+5)
            y(i+6)=x(i+6)
         EndDo
      EndIf
      End
!+ dot product
      Complex*16 Function ZDOTU(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      Complex*16 x(*),y(*)
      ZDOTU=0
      If (n.le.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            ZDOTU=ZDOTU+x(j)*y(k)
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            ZDOTU=ZDOTU+x(i)*y(i)
         Enddo
         Do i=m+1,n,5
            ZDOTU=ZDOTU+x(i)*y(i)+y(i+1)*x(i+1)+
     1     x(i+2)*y(i+2)+y(i+3)*x(i+3)+x(i+4)*y(i+4)
         EndDo
      EndIf
      End
!CIF
!+ dot product
      Complex*16 Function ZDOTC(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      Complex*16 x(*),y(*)
      ZDOTC=0
      If (n.le.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            ZDOTC=ZDOTC+conjg(x(j))*y(k)
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            ZDOTC=ZDOTC+conjg(x(i))*y(i)
         Enddo
         Do i=m+1,n,5
            ZDOTC=ZDOTC+conjg(x(i))*y(i)
     1+conjg(x(i+1))*y(i+1)
     2+conjg(x(i+2))*y(i+2)
     3+conjg(x(i+3))*y(i+3)
     4+conjg(x(i+4))*y(i+4)
         EndDo
      Endif
      End
!CEND
!+ apply a plane rotation, backwards if negative increment
      Subroutine ZDROT(n,x,incx,y,incy,c,s)
      Integer i,incx,incy,j,k,n
      Complex*16 x(1),y(1),temp
      Real*8 c,s
      If (n.le.0.or.s.eq.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            temp=
     1dcmplx(dreal(x(j))*c+dreal(y(k))*s,dimag(x(j))*c+dimag(y(k))*s)
            y(k)=
     1dcmplx(-dreal(x(j))*s+dreal(y(k))*c,-dimag(x(j))*s+dimag(y(k))*c)
            x(j)=temp
            j=j+incx
            k=k+incy
         EndDo
      Else
         Do i=1,n
            temp=
     1 dcmplx(dreal(x(i))*c+dreal(y(i))*s,dimag(x(i))*c+dimag(y(i))*s)
            y(i)=
     1 dcmplx(-dreal(x(i))*s+dreal(y(i))*c,-dimag(x(i))*s+dimag(y(i))*c)
            x(i)=temp
         EndDo
      Endif
      End
!+ construct Givens plane rotation
      Subroutine ZROTG(a,b,c,s)
      Real*8 c,aa,bb
      Complex*16 a,b,s,r
      bb=abs(dreal(b))+abs(dimag(b))
      If (bb.gt.0) Then
         aa=abs(dreal(a))+abs(dimag(a))
         If (aa.gt.bb) Then
            r=a*sqrt(1+(bb/aa)**2)
            c=real(a/r)
            s=b/r
            b=s
            a=r
         Else
            r=b*sqrt(1+(aa/bb)**2)
            c=real(a/r)
            s=b/r
            a=r
            b=1
            If (c.ne.0) b=1/c
         Endif
      Else
            c=1
            s=0
      Endif
      End
!CIF
      Subroutine ZDSCAL(n,a,x,incx)
      Integer i,incx,m,n
      Real*8 a
      Complex*16 x(1)
      If (n.le.0) Return
      If (incx.ne.1) Then
         Do i=1,n*incx,incx
            x(i)=dcmplx(dreal(x(i))*a,dimag(x(i))*a)
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            x(i)=dcmplx(dreal(x(i))*a,dimag(x(i))*a)
         EndDo
         Do i=m+1,n,5
            x(i)=dcmplx(dreal(x(i))*a,dimag(x(i))*a)
            x(i+1)=dcmplx(dreal(x(i+1))*a,dimag(x(i+1))*a)
            x(i+2)=dcmplx(dreal(x(i+2))*a,dimag(x(i+2))*a)
            x(i+3)=dcmplx(dreal(x(i+3))*a,dimag(x(i+3))*a)
            x(i+4)=dcmplx(dreal(x(i+4))*a,dimag(x(i+4))*a)
         EndDo
      Endif
      End
!CEND
!+ x=a*x  scale
      Subroutine ZSCAL(n,a,x,incx)
      Integer i,incx,m,n
      Complex*16 x(1),a
      If (n.le.0) Return
      If (incx.ne.1) Then
         Do i=1,n*incx,incx
            x(i)=x(i)*a
         EndDo
      Else
         m=mod(n,5)
         Do i=1,m
            x(i)=x(i)*a
         EndDo
         Do i=m+1,n,5
            x(i)=x(i)*a
            x(i+1)=x(i+1)*a
            x(i+2)=x(i+2)*a
            x(i+3)=x(i+3)*a
            x(i+4)=x(i+4)*a
         EndDo
      EndIf
      End
!+ swap arrays, backwards if negative increment
      Subroutine ZSWAP(n,x,incx,y,incy)
      Integer i,incx,incy,j,k,m,n
      Complex*16 x(1),y(1),temp
      If (n.le.0) Return
      If (incx.ne.1.or.incy.ne.1) Then
         j=1
         k=1
         If (incx.lt.0) j=(1-n)*incx+1
         If (incy.lt.0) k=(1-n)*incy+1
         Do i=1,n
            temp=x(j)
            x(j)=y(k)
            y(k)=temp
            j=j+incx
            k=k+incy
         EndDo
      Else
         m=mod(n,3)
         Do i=1,m
            temp=x(i)
            x(i)=y(i)
            y(i)=temp
         EndDo
         Do i=m+1,n,3
            temp=x(i)
            x(i)=y(i)
            y(i)=temp
            temp=x(i+1)
            x(i+1)=y(i+1)
            y(i+1)=temp
            temp=x(i+2)
            x(i+2)=y(i+2)
            y(i+2)=temp
         EndDo
      EndIf
      End

max_stars_repo_name: fluffynukeit/mdsplus

max_issues_repo_name: fluffynukeit/mdsplus

max_forks_repo_name: fluffynukeit/mdsplus


text:
!=======================================================================
!  GrN_Ptl, Subroutine
!
!       Calculates N available for transfer to grain (g/plant) from each
!       plant part.  By definition, available grain N is set to 0.
!
!       Notes
!           N available for translocation to the grain is the sum of N
!              available in the stover and in the roots.
!           N available in stover is the difference of its N content and
!              the minimum it's allowed to fall to.
!           N available in roots is the product of root weight and the
!              fraction of root N available for translocation.
!-----------------------------------------------------------------------
!  Revision history
!  12-21-2011 FSR Created from NWheat subroutine "nwheats_gnptl"
!-----------------------------------------------------------------------
!  Called by: WH_GroSub
!=======================================================================
      SUBROUTINE GrN_Ptl (CONTROL, ISWITCH,
     &   mnc, nfact, nitmn, npot, optfr, part, pl_la, pl_nit,     !INPUT
     &   plantwt, sen_la,                                         !INPUT
     &   navl)                                                   !OUTPUT
!-----------------------------------------------------------------------
      USE ModuleDefs
      USE WH_module
      IMPLICIT NONE
      SAVE
!-----------------------------------------------------------------------
!                             Define Variables
!-----------------------------------------------------------------------
      INTEGER    DYNAMIC
      INTEGER    part   ! plant part number
      REAL       mnc(mxpart)  ! Minimum N concentration, by plant part
      REAL       navl(mxpart) ! N in each plant part available to grain?
      REAL       nfact(4)
      REAL       nitmn  ! nitrogen minimum level              (g/part)
      REAL       npot   ! potential N available for transfer  (g/part)
      REAL       optfmn ! optimum fraction minimum (0-1)
      REAL       optfmx ! optimum fraction maximum (0-1)
      REAL       optfr  ! optimum fraction due to stress      (0-1)
      REAL       pl_la  ! Plant leaf area mm2
      REAL       pl_nit(mxpart) !plant part nitrogen (g/plant)
      REAL       plantwt(mxpart) !plant part weights (g/plant)
      REAL       sen_la  ! Senesced leaf area

      PARAMETER (optfmn = 0.15) !*! should be a (species?) input
      PARAMETER (optfmx = 0.35) !*! should be a (species?) input

      REAL       g_gn_optfr ! this is a temp variable; remove ASAP

!     The variable "CONTROL" is of type "ControlType".
      TYPE (ControlType) CONTROL
      TYPE (ResidueType) SENESCE
      TYPE (SwitchType)  ISWITCH

      DYNAMIC = CONTROL % DYNAMIC
!-----------------------------------------------------------------------

!*! following moved to ModuleDefs.for
!*!   integer    mxpart !*! maximum number of plant parts
!*!   parameter (mxpart = 6) !*! counted root, stem, leaf sheath,
!*!             leaf, grain and seed: could not find any reference

*+  Local Variables	 *******OK*******
! -------------- get grain N potential (supply) ----------
!                get the fraction of optimum conditions

      optfr  =  optfmn + (optfmx - optfmn)*nfact(2)

cnh added for watch purposes
      g_gn_optfr = optfr

      ! Now find the available N of each part.
      DO 1000 part = 1, mxpart
         IF (pl_la - sen_la .gt. 0.0) THEN
      ! Calculate N available for transfer to grain
            nitmn = mnc(part) * plantwt(part)
            npot = MAX(pl_nit(part) - nitmn, 0.0)
            navl(part)  =  optfr * npot
         ELSE  ! No green leaf area. No further N transfer to grain
      navl(part) = 0.0
         ENDIF

1000  continue
!*!      ENDDO

      navl(grain_part) = 0.0
      navl(seed_part) = 0.0

      RETURN
      END SUBROUTINE GrN_Ptl
!-----------------------------------------------------------------------
C--------------------------------------------------------------------------------------------------
C                         Define Variables
C--------------------------------------------------------------------------------------------------

! ANDEM       !Crop N demand (kg N/ha)
! ANO3        !Total extractable nitrate N in soil profile (kg N/ha)
! ANH4        !Total extractable ammonium N in soil profile (kg N/ha)
! DLAYR(L)    Soil thickness in layer L (cm)
! DNG         !N demand of potential new growth of tops (g N/plant)
! DROOTN      !Daily change in plant root nitrogen content (g N/plant)
! DSTOVN      !Daily change in plant stover nitrogen content (g N/plant)
! FACTOR      !Ratio of root N uptake to plant N demand
! FNH4        !Unitless soil ammonium supply index
! FNO3        !Unitless soil nitrate supply index
! FON(20)     !Fresh organic nitrogen in soil layer L due to root senescence, kg N/ha
! KG2PPM(20)  !Factor that convergs mg elemental N/kg soil to kg N/ha for soil layer L
! L           !Index counter
! L1          !Lowest soil layer with roots
! LL(20)      !Lower limit of plant extractable water for soil layer L, cm3/cm3
! mnc(mxpart) ! Minimum N concentration, by plant part
! NDEM        !Plant nitrogen demand, g/plant
! NH4(20)     !Ammonium in soil layer L, mg elemental N/kg soil
! NLAYR       !Number of soil layer
! NO3(20)     !Nitrate in soil layer L (mg elemental N/kg soil)
! NUF         !Plant N supply to demand used to modify N uptake
! PDWI        !Potential increment in new shoot growth, g/plant
! PGRORT      !Potential increment in new root growth, g/plant
! pl_la       !Plant leaf area mm2
! PLTPOP      !Plant population, plants/m2
! PTF         !Ratio of above ground biomass to total biomass
! RANC        !Root actual N concentration, g N/g root
! RCNP        !Root critical nitrogen concentration, g N/g root dry weight
! RFAC        !Interim variable describing the effects of root length density on potential N uptake from a layer
! RLV(20)     !Root length density for soil layer L, cm root/cm2 soil
! RNDEM       !Plant root demand for nitrogen (g/plant)
! RNH4U(20)   !Potential ammonia uptake from layer L, kg N/ha
! RNLOSS      !Loss of N from the plant via root exudation in one layer (g N/m2)
! RNO3U(20)   !Potential nitrate uptake from layer L, kg N/ha
! ROOTN       !Root nitrogen content, g N/plant
! RTWT        !Root weight, g/plant
! sen_la      ! Senesced leaf area
! SHF(20)     !Relative root distribution in soil layer L (0-1)
! SMDFR       !Soil moisture deficit factor affecting N uptake
! SNH4(20)    !Ammonium nitrogen in soil layer L, kg N/ha
! SNO3(20)    !Nitrate content in soil layer L, kg N/ha
! STOVN       !Nitrogen content in stover, g N/plant
! STOVWT      !Stover weight (Stem + leaf+lfsheath), g/plant
! SAT(20)     !Saturated water holding capacity for soil layer L, cm3/cm3
! SW(20)      !Soil water content in layer L, cm3 water/cm3 soil
! TANC        !Nitrogen content in above ground biomass, decimal
! TCNP        !Critical nitrogen concentration in tops, g N/g dry weight
! TNDEM       !Plant tops demand for nitrogen (g N/plant)
! TRLV        !Total root length density, cm root/cm2 soil
! TRNLOS      !Total plant N lost by root exudation (g N/m2)
! TRNU        !Total potential root nitrogen uptake, kg N/ha
! UNH4(20)    !Plant uptake of ammonium from layer (kg N/ha/day)
! UNO3(20)    !Plant uptake of nitrate from a layer (kg N/ha/day)
! XMIN        !
! XNDEM       !Temporary variable
! xstag_nw    !Non-integer growth stage indicator (NWheats)
! XSTAGE      !Non-integer growth stage indicator

ext: for

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


      subroutine nut_nminrl

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine estimates daily nitrogen and phosphorus
!!    mineralization and immobilization considering fresh organic
!!    material (plant residue) and active and stable humus material

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ihru          |none          |HRU number
!!    rsdco_pl(:)   |none          |plant residue decomposition coefficient. The
!!                                 |fraction of residue which will decompose in
!!                                 |a day assuming optimal moisture,
!!                                 |temperature, C:N ratio, and C:P ratio
!!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Max, Exp, Sqrt, Min, Abs

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use septic_data_module
      use plant_data_module
      use basin_module
      use organic_mineral_mass_module
      use hru_module, only : rsdco_plcom, i_sep, ihru, ipl, isep
      use soil_module
      use plant_module
      use output_landscape_module, only : hnb_d

      implicit none

      integer :: j          !none          |HRU number
      integer :: k          !none          |counter (soil layer)
      integer :: kk         !none          |soil layer used to compute soil water and
                            !              |soil temperature factors
      integer :: idp
      real :: rmn1          !kg N/ha       |amount of nitrogen moving from fresh organic
                            !              |to nitrate(80%) and active organic(20%)
                            !              |pools in layer
      real :: rmp           !              |to labile(80%) and organic(20%) pools in layer
      real :: xx            !varies        |variable to hold intermediate calculation result
      real :: csf           !none          |combined temperature/soil water factor
      real :: rwn           !kg N/ha       |amount of nitrogen moving from active organic
                            !              |to stable organic pool in layer
      real :: hmn           !kg N/ha       |amount of nitrogen moving from active organic
                            !              |nitrogen pool to nitrate pool in layer
      real :: hmp           !kg P/ha       |amount of phosphorus moving from the organic
                            !              |pool to the labile pool in layer
      real :: r4            !              |
      real :: cnr           !              |
      real :: cnrf          !              |
      real :: cpr           !              |
      real :: cprf          !              |
      real :: ca            !              |
      real :: decr          !              |
      real :: rdc           !              |
      real :: wdn           !kg N/ha       |amount of nitrogen lost from nitrate pool in
                            !              |layer due to denitrification
      real :: cdg           !none          |soil temperature factor
      real :: sut           !none          |soil water factor
      real :: nactfr        !none          |nitrogen active pool fraction. The fraction
                            !              |of organic nitrogen in the active pool.

      j = ihru
      nactfr = .02
      !zero transformations for summing layers
      hnb_d(j)%rsd_nitorg_n = 0.
      hnb_d(j)%rsd_laborg_p = 0.
      hnb_d(j)%act_nit_n = 0.
      hnb_d(j)%org_lab_p = 0.
      hnb_d(j)%act_sta_n = 0.
      hnb_d(j)%denit = 0.

      !! mineralization can occur only if temp above 0 deg
      if (soil(j)%phys(1)%tmp > 0.) then
      !! compute residue decomp and mineralization of fresh organic n and p of flat residue
        do ipl = 1, pcom(j)%npl        !! we need to decompose each plant
          rmn1 = 0.
          rmp = 0.
          r4 = .58 * rsd1(j)%tot(ipl)%m
          if (rsd1(j)%tot(ipl)%n + rsd1(j)%mn%no3 > 1.e-4) then
            cnr = r4 / (rsd1(j)%tot(ipl)%n  + rsd1(j)%mn%no3)
            if (cnr > 500.) cnr = 500.
            cnrf = Exp(-.693 * (cnr - 25.) / 25.)
          else
            cnrf = 1.
          end if

          if (rsd1(j)%tot(ipl)%p + rsd1(j)%mp%lab > 1.e-4) then
            cpr = r4 / (rsd1(j)%tot(ipl)%p + rsd1(j)%mp%lab)
            if (cpr > 5000.) cpr = 5000.
            cprf = Exp(-.693 * (cpr - 200.) / 200.)
          else
            cprf = 1.
          end if

          !! compute soil water factor
	      if (soil(j)%phys(1)%st < 0.) soil(j)%phys(1)%st = .0000001
          sut = .1 + .9 * Sqrt(soil(j)%phys(1)%st / soil(j)%phys(1)%fc)
          sut = Max(.05, sut)

          !!compute soil temperature factor
          xx = soil(j)%phys(1)%tmp
          cdg = .9 * xx / (xx + Exp(9.93 - .312 * xx)) + .1
          cdg = Max(.1, cdg)

          !! compute combined factor
          xx = cdg * sut
          if (xx < 0.) xx = 0.
          if (xx > 1.e6) xx = 1.e6
          csf = Sqrt(xx)
          ca = Min(cnrf, cprf, 1.)
          !! compute residue decomp and mineralization for each plant
          if (pcom(j)%npl > 0) then
            idp = pcom(j)%plcur(ipl)%idplt
            decr = pldb(idp)%rsdco_pl * ca * csf
          else
            decr = 0.05
          end if
          decr = Max(bsn_prm%decr_min, decr)
          decr = Min(decr, 1.)
          rsd1(j)%tot(ipl)%m = Max(1.e-6, rsd1(j)%tot(ipl)%m)
          rdc = decr * rsd1(j)%tot(ipl)%m
          rsd1(j)%tot(ipl)%m = rsd1(j)%tot(ipl)%m - rdc
          if (rsd1(j)%tot(ipl)%m < 0.) rsd1(j)%tot(ipl)%m = 0.
          rmn1 = decr * rsd1(j)%tot(ipl)%n
          rsd1(j)%tot(ipl)%p = Max(1.e-6, rsd1(j)%tot(ipl)%p)
          rmp = decr * rsd1(j)%tot(ipl)%p

          rsd1(j)%tot(ipl)%p = rsd1(j)%tot(ipl)%p - rmp
          rsd1(j)%tot(ipl)%n  = Max(1.e-6, rsd1(j)%tot(ipl)%n)
          if (j==1 .and. rsd1(j)%tot(ipl)%n > 100) then
            decr = 0.05
          end if
          rsd1(j)%tot(ipl)%n  = rsd1(j)%tot(ipl)%n  - rmn1
          rsd1(j)%mn%no3 = rsd1(j)%mn%no3 + .8 * rmn1
          soil1(j)%hs(1)%n = soil1(j)%hs(1)%n + .2 * rmn1
          rsd1(j)%mp%lab = rsd1(j)%mp%lab + .8 * rmp
          soil1(j)%hp(1)%p = soil1(j)%hp(1)%p + .2 * rmp

          hnb_d(j)%rsd_nitorg_n = hnb_d(j)%rsd_nitorg_n + rmn1
          hnb_d(j)%rsd_laborg_p = hnb_d(j)%rsd_laborg_p + rmp

        end do      ! ipl = 1, pcom(j)%npl
      end if

      !! compute residue decomp and mineralization of fresh organic n and p
      !! root and incorporated residue
      do k = 1, soil(j)%nly

        if (k == 1) then
          kk = 2
        else
          kk = k
        end if

        !! mineralization can occur only if temp above 0 deg
        if (soil(j)%phys(kk)%tmp > 0.) then
          !! compute soil water factor
          sut = 0.
	      !! change for domain error 1/29/09 gsm check with Jeff !!!
	      if (soil(j)%phys(kk)%st < 0.) soil(j)%phys(kk)%st = .0000001
          sut = .1 + .9 * Sqrt(soil(j)%phys(kk)%st / soil(j)%phys(kk)%fc)
          sut = Max(.05, sut)

          !!compute soil temperature factor
          xx = soil(j)%phys(kk)%tmp
          cdg = .9 * xx / (xx + Exp(9.93 - .312 * xx)) + .1
          cdg = Max(.1, cdg)

          !! compute combined factor
          xx = cdg * sut
          if (xx < 0.) xx = 0.
          if (xx > 1.e6) xx = 1.e6
          csf = Sqrt(xx)

          !! compute flow from active to stable pools
          !rwn = .1e-4 * (soil1(j)%hs(k)%n * (1. / nactfr - 1.) - soil1(j)%hp(k)%n)
          rwn = .1e-4 * (soil1(j)%hs(k)%n - soil1(j)%hp(k)%n)
          if (rwn > 0.) then
            rwn = Min(rwn, soil1(j)%hs(k)%n)
          else
            rwn = -(Min(Abs(rwn), soil1(j)%hp(k)%n))
          endif
          soil1(j)%hp(k)%n = Max(1.e-6, soil1(j)%hp(k)%n + rwn)
          soil1(j)%hs(k)%n = Max(1.e-6, soil1(j)%hs(k)%n - rwn)
          hnb_d(j)%act_sta_n = hnb_d(j)%act_sta_n + rwn

          !! compute humus mineralization on active organic n
          hmn = bsn_prm%cmn * csf * soil1(j)%hs(k)%n
          hmn = Min(hmn, soil1(j)%hs(k)%n)
          !! compute humus mineralization on active organic p
          xx = soil1(j)%hp(k)%n + soil1(j)%hs(k)%n
          if (xx > 1.e-6) then
            hmp = 1.4 * hmn * soil1(j)%hp(k)%p / xx
          else
            hmp = 0.
          end if
          hmp = Min(hmp, soil1(j)%hp(k)%p)
          !! move mineralized nutrients between pools
          soil1(j)%hs(k)%n = Max(1.e-6, soil1(j)%hs(k)%n - hmn)
          soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 + hmn
          soil1(j)%hp(k)%p = soil1(j)%hp(k)%p - hmp
          soil1(j)%mp(k)%lab = soil1(j)%mp(k)%lab + hmp

          hnb_d(j)%act_nit_n = hnb_d(j)%act_nit_n + hmn
          hnb_d(j)%org_lab_p = hnb_d(j)%org_lab_p + hmp

          !! compute residue decomp and mineralization of
          !! fresh organic n and p (upper two layers only)
          rmn1 = 0.
          rmp = 0.
          if (k <= 2) then
            r4 = .58 * soil1(j)%tot(k)%m
            if (soil1(j)%tot(k)%n + soil1(j)%mn(k)%no3 > 1.e-4) then
              cnr = r4 / (soil1(j)%tot(k)%n  + soil1(j)%mn(k)%no3)
              if (cnr > 500.) cnr = 500.
              cnrf = Exp(-.693 * (cnr - 25.) / 25.)
            else
              cnrf = 1.
            end if

            if (soil1(j)%tot(k)%p + soil1(j)%mp(k)%lab > 1.e-4) then
              cpr = r4 / (soil1(j)%tot(k)%p + soil1(j)%mp(k)%lab)
              if (cpr > 5000.) cpr = 5000.
              cprf = Exp(-.693 * (cpr - 200.) / 200.)
            else
              cprf = 1.
            end if

            ca = Min(cnrf, cprf, 1.)

            !! compute residue decomp and mineralization for each plant
            decr = .05
            !do ipl = 1, pcom(j)%npl        !! we need to decomp each plant**

            if (pcom(j)%npl > 0) then
              decr = rsdco_plcom(j) * ca * csf
            else
              decr = 0.05
            end if
            decr = Max(bsn_prm%decr_min, decr)
            decr = Min(decr, 1.)
            soil1(j)%tot(k)%m = Max(1.e-6, soil1(j)%tot(k)%m)
            rdc = decr * soil1(j)%tot(k)%m
            soil1(j)%tot(k)%m = soil1(j)%tot(k)%m - rdc
            if (soil1(j)%tot(k)%m < 0.) soil1(j)%tot(k)%m = 0.
            rmn1 = decr * (soil1(j)%str(k)%n + soil1(j)%lig(k)%n + soil1(j)%meta(k)%n)
            soil1(j)%tot(k)%p = Max(1.e-6, soil1(j)%tot(k)%p)
            rmp = decr * (soil1(j)%str(k)%p + soil1(j)%lig(k)%p + soil1(j)%meta(k)%p)

            soil1(j)%tot(k)%p = soil1(j)%tot(k)%p - rmp
            soil1(j)%tot(k)%n  = Max(1.e-6, soil1(j)%tot(k)%n)
            soil1(j)%tot(k)%n  = soil1(j)%tot(k)%n  - rmn1
            soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 + .8 * rmn1
            soil1(j)%hs(k)%n = soil1(j)%hs(k)%n + .2 * rmn1
            soil1(j)%mp(k)%lab = soil1(j)%mp(k)%lab + .8 * rmp
            soil1(j)%hp(k)%p = soil1(j)%hp(k)%p + .2 * rmp

            hnb_d(j)%rsd_nitorg_n = hnb_d(j)%rsd_nitorg_n + rmn1
            hnb_d(j)%rsd_laborg_p = hnb_d(j)%rsd_laborg_p + rmp
          end if

          !!  compute denitrification
          wdn = 0.
	      if (i_sep(j) /= k .or. sep(isep)%opt  /= 1) then
	        if (sut >= bsn_prm%sdnco) then
	          wdn = soil1(j)%mn(k)%no3 * (1.-Exp(-bsn_prm%cdn * cdg * soil1(j)%cbn(k) / 100.))
	        else
	          wdn = 0.
	        endif
	        soil1(j)%mn(k)%no3 = soil1(j)%mn(k)%no3 - wdn
          end if
          hnb_d(j)%denit = hnb_d(j)%denit + wdn

          !call nut_denit(k,j,cdg,wdn,0.05)

        end if
      end do        ! k = 1, soil(j)%nly

      return
      end subroutine nut_nminrl


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 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 checkprojectgrad(nactiveold,nactive,ipoacti,ipoactiold,
     &   objectset,xlambd,nnlconst,iconstacti,iconstactiold,inameacti,
     &   inameactiold,g0,nobject,ndesi,nodedesi,dgdxglob,nk)
!
!     checks the lagrange multipliers and reduces the number of active
!     constraints if possible and update the values of the response
!     functions for the linear constraints
!
      implicit none
!
      character*81 objectset(5,*)
!
      integer nactiveold,nactive,ipoacti(*),ipoactiold(*),i,j,
     &   nnlconst,nnlconstold,iconstacti(*),iconstactiold(*),
     &   inameacti(*),inameactiold(*),iact,nobject,ndesi,
     &   nodedesi(*),node,nk
!
      real*8 xlambd(*),g0(*),dgdxglob(2,nk,*),val
!
!
      do i=1,nactive
         ipoactiold(i)=ipoacti(i)
         iconstactiold(i)=iconstacti(i)
         inameactiold(i)=inameacti(i)
      enddo
!
      nactiveold=nactive
      nnlconstold=nnlconst
      nactive=0
      nnlconst=0
!
      do i=1,nactiveold
!
!        check of all nonlinear constraints
!
         if(i.le.nnlconstold) then
            if(iconstactiold(i).eq.-1) then
               if(xlambd(i).gt.0.d0) then
                  nactive=nactive+1
                  nnlconst=nnlconst+1
                  ipoacti(nactive)=ipoactiold(i)
                  iconstacti(nactive)=iconstactiold(i)
                  inameacti(nactive)=inameactiold(i)
               endif
            elseif(iconstactiold(i).eq.1) then
               if(xlambd(i).lt.0.d0) then
                  nactive=nactive+1
                  nnlconst=nnlconst+1
                  ipoacti(nactive)=ipoactiold(i)
                  iconstacti(nactive)=iconstactiold(i)
                  inameacti(nactive)=inameactiold(i)
               endif
            endif
!
!        check of all geometric (linear) constraints
!
         else
!
!           MAXMEMBERSIZE and MINMEMBERSIZE
!
            if(objectset(1,inameactiold(i))(4:13).eq.'MEMBERSIZE') then
               node=nodedesi(ipoactiold(i))
               val=dgdxglob(2,node,inameactiold(i))
               if(iconstactiold(i).eq.-1) then
                  if(((xlambd(i).gt.0.d0).and.(val.lt.0.d0)).or.
     &               ((xlambd(i).lt.0.d0).and.(val.gt.0.d0))) then
                     nactive=nactive+1
                     ipoacti(nactive)=ipoactiold(i)
                     iconstacti(nactive)=iconstactiold(i)
                     inameacti(nactive)=inameactiold(i)
                  endif
               elseif(iconstactiold(i).eq.1) then
                  if(((xlambd(i).lt.0.d0).and.(val.lt.0.d0)).or.
     &               ((xlambd(i).gt.0.d0).and.(val.gt.0.d0))) then
                     nactive=nactive+1
                     ipoacti(nactive)=ipoactiold(i)
                     iconstacti(nactive)=iconstactiold(i)
                     inameacti(nactive)=inameactiold(i)
                  endif
               endif
!
!           FIXGROWTH and FIXSHRINKAGE
!
            else
               if(iconstactiold(i).eq.-1) then
                  if(xlambd(i).gt.0.d0) then
                     nactive=nactive+1
                     ipoacti(nactive)=ipoactiold(i)
                     iconstacti(nactive)=iconstactiold(i)
                     inameacti(nactive)=inameactiold(i)
                  endif
               elseif(iconstactiold(i).eq.1) then
                  if(xlambd(i).lt.0.d0) then
                     nactive=nactive+1
                     ipoacti(nactive)=ipoactiold(i)
                     iconstacti(nactive)=iconstactiold(i)
                     inameacti(nactive)=inameactiold(i)
                  endif
               endif
            endif
         endif
      enddo
!
!     update the values of the response functions and
!     the sensitivities for the linear constraint
!
      do i=1,nobject-1
         if(objectset(5,i)(81:81).eq.'G') then
            iact=0
            do j=1,ndesi
               node=nodedesi(j)
               dgdxglob(2,node,i)=0.d0
            enddo
            do j=1,nactive
               if(inameacti(j).eq.i) then
                  iact=iact+1
                  node=nodedesi(ipoacti(j))
                  dgdxglob(2,node,i)=1.d0
               endif
            enddo
            g0(i)=iact
         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:
	LOGICAL FUNCTION RS_NMER ( lets, numb )
C************************************************************************
C* RS_NMER								*
C*									*
C* This function determines whether the input letters are ASCII digits.	*
C*									*
C* LOGICAL RS_NMER  ( LETS, NUMB )					*
C*									*
C* Input parameters:							*
C*	LETS		CHAR*		String				*
C*	NUMB		INTEGER		Number of letters		*
C*									*
C* Output parameters:							*
C*	RS_NMER		LOGICAL		Return value			*
C*					= .TRUE. if letters are digits, *
C*						between '0' and '9'	*
C*					= .FALSE. if any letter is not 	*
C*						a digit			*
C**									*
C* Log:									*
C* J. Nielsen/MIT	10/86						*
C* J. Nielsen/TAMU	 2/92	Gempacized				*
C************************************************************************
	CHARACTER*(*)	lets
	INTEGER		numb
C------------------------------------------------------------------------
	RS_NMER = .TRUE.
	DO  i = 1, numb
	    IF  ( ( lets(i:i) .lt. '0' ) .or.
     +		  ( lets(i:i) .gt. '9' ) ) THEN
	        RS_NMER = .FALSE.
	        RETURN
	    ENDIF
	END DO
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C **********************************************************************
C *                                                                    *
C * ISNUM - RETURN .TRUE. IF NUMERIC                                   *
C *                                                                    *
C **********************************************************************
C
      LOGICAL FUNCTION ISNUM(C)
      INTEGER C
C
      ISNUM = .FALSE.
      IF (C .GE. 48 .AND. C .LE. 57) ISNUM = .TRUE.
      RETURN
      END


      ! Compute the external source terms. This is called twice in the evolution:
      !
      ! 1) For the predictor, it is called with (old, old) states
      !    This is also used in the first pass of the conservative update
      !    (adding dt * S there)
      !
      ! 2) Next we correct the source terms in the conservative update to
      ! time-center them.  Here we call ext_src(old, new), and then
      ! in time_center_source_terms we subtract off 1/2 of the first S
      ! and add 1/2 of the new S.
      !
      ! Therefore, to get a properly time-centered source, generally
      ! speaking, you always want to use the "new" state here.  That
      ! will be the time n state in the first call and the n+1 in the
      ! second call.
      !

      subroutine ext_src(lo,hi,&
                         old_state,old_state_l1,old_state_l2,old_state_l3,old_state_h1,old_state_h2,old_state_h3,&
                         new_state,new_state_l1,new_state_l2,new_state_l3,new_state_h1,new_state_h2,new_state_h3,&
                         old_diag ,old_diag_l1,old_diag_l2,old_diag_l3,old_diag_h1,old_diag_h2,old_diag_h3,&
                         new_diag ,new_diag_l1,new_diag_l2,new_diag_l3,new_diag_h1,new_diag_h2,new_diag_h3,&
                         src,src_l1,src_l2,src_l3,src_h1,src_h2,src_h3,problo,dx,time,z,dt) &
                         bind(C, name="fort_ext_src")

      use amrex_fort_module, only : rt => amrex_real
      use meth_params_module, only : NVAR

      implicit none

      integer ,intent(in   ) :: lo(3),hi(3)
      integer ,intent(in   ) :: old_state_l1,old_state_l2,old_state_l3,old_state_h1,old_state_h2,old_state_h3
      integer ,intent(in   ) :: new_state_l1,new_state_l2,new_state_l3,new_state_h1,new_state_h2,new_state_h3
      integer ,intent(in   ) :: old_diag_l1,old_diag_l2,old_diag_l3,old_diag_h1,old_diag_h2,old_diag_h3
      integer ,intent(in   ) :: new_diag_l1,new_diag_l2,new_diag_l3,new_diag_h1,new_diag_h2,new_diag_h3
      integer ,intent(in   ) :: src_l1,src_l2,src_l3,src_h1,src_h2,src_h3
      real(rt),intent(in   ) :: old_state(old_state_l1:old_state_h1,old_state_l2:old_state_h2, &
                                                 old_state_l3:old_state_h3,NVAR)
      real(rt),intent(in   ) :: new_state(new_state_l1:new_state_h1,new_state_l2:new_state_h2, &
                                                 new_state_l3:new_state_h3,NVAR)
      real(rt),intent(in   ) :: old_diag(old_diag_l1:old_diag_h1,old_diag_l2:old_diag_h2, &
                                                 old_diag_l3:old_diag_h3,NVAR)
      real(rt),intent(in   ) :: new_diag(new_diag_l1:new_diag_h1,new_diag_l2:new_diag_h2, &
                                                 new_diag_l3:new_diag_h3,NVAR)
      real(rt),intent(  out) :: src(src_l1:src_h1,src_l2:src_h2,src_l3:src_h3,NVAR)
      real(rt),intent(in   ) :: problo(3),dx(3),time,z,dt

      call ext_src_force(lo,hi,&
                         old_state,old_state_l1,old_state_l2,old_state_l3,old_state_h1,old_state_h2,old_state_h3,&
                         new_state,new_state_l1,new_state_l2,new_state_l3,new_state_h1,new_state_h2,new_state_h3,&
                         old_diag ,old_diag_l1,old_diag_l2,old_diag_l3,old_diag_h1,old_diag_h2,old_diag_h3,&
                         new_diag ,new_diag_l1,new_diag_l2,new_diag_l3,new_diag_h1,new_diag_h2,new_diag_h3,&
                         src,src_l1,src_l2,src_l3,src_h1,src_h2,src_h3,problo,dx,time,z,dt)

      end subroutine ext_src