text:
	SUBROUTINE OAGOPN  ( filnam, source, dattim, iflno,
     +			     time, ivert, iret )
C************************************************************************
C* OAGOPN								*
C*									*
C* This subroutine opens the file to be used for computing the 		*
C* station spacing.							*
C*									*
C* OAGOPN  ( FILNAM, SOURCE, DATTIM, IFLNO, TIME, IVERT,		*
C*           IRET )							*
C*									*
C* Input parameters:							*
C*	FILNAM		CHAR*		File name			*
C*	SOURCE		CHAR*		Data type ( SN or SF )		*
C*	DATTIM		CHAR*		Input date/time			*
C*									*
C* Output parameters:							*
C*	IFLNO		INTEGER		File number			*
C*	TIME		CHAR*		Date/time used			*
C*	IVERT		INTEGER		File vertical coordinate	*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					-11 = file not opened		*
C*					-12 = invalid time		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 8/85						*
C* M. desJardins/GSFC	11/88	GEMPAK 4.1				*
C* S. Schotz/GSC	 6/90	Removed respnd flag			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	filnam, source, dattim, time
C*
	LOGICAL		mrgdat
	CHARACTER	timlst (LLMXTM)*20, chparm (MMPARM)*4
C------------------------------------------------------------------------
	iret = 0
C
C*	Open new file.
C
	IF  ( source .eq. 'SN' )  THEN
	    CALL SN_OPNF  ( filnam, .false., iflno, isrc, nparm,
     +			    chparm, ivert, mrgdat, iret )
	    IF  ( iret .ne. 0 ) THEN
		iret = -11
		CALL ER_WMSG  ( 'OAGRID', iret, filnam, ier )
		RETURN
	      ELSE
	        CALL PC_INIT  ( ivert, nparm, chparm, ier )
	        CALL SN_BEGS  ( iflno, ier )
	    END IF
C*
	  ELSE
	    CALL SF_OPNF  ( filnam, .false., iflno, isrc, nparm, chparm,
     +			    iret )
	    ivert = 0
	    IF  ( iret .ne. 0 ) THEN
		iret = -11
		CALL ER_WMSG  ( 'OAGRID', iret, filnam, ier )
		RETURN
	      ELSE
	        CALL PC_INIT  ( 0, nparm, chparm, ier )
	        CALL SF_BEGS  ( iflno, ier )
	    END IF
	END IF
C
C*	Get the times in the file and set the time correctly.
C
	IF  ( source .eq. 'SN' )  THEN
	    CALL SN_GTIM  ( iflno, LLMXTM, numtim, timlst, ier )
	    CALL TI_FIND  ( dattim, numtim, timlst, dattim,
     +			    nt, time, ier )
	    IF  ( ( ier .ne. 0 ) .or. ( nt .lt. 1 ) )  THEN
		iret = -12
		CALL ER_WMSG  ( 'OAGRID', iret, dattim, ier )
		RETURN
	    END IF
C
C*	    Set the time.
C
	    CALL SN_STIM  ( iflno, time, iret )
	  ELSE
	    CALL SF_GTIM  ( iflno, LLMXTM, numtim, timlst, iret )
	    CALL TI_FIND  ( dattim, numtim, timlst, dattim,
     +			    nt, time, ier )
	    IF  ( ( ier .ne. 0 ) .or. ( nt .lt. 1))  THEN
		iret = -12
		CALL ER_WMSG  ( 'OAGRID', iret, dattim, ier )
	    END IF
C
C*	    Set the time.
C
	    CALL SF_STIM (iflno, time, iret )
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine jac(neq,t,y,ml,mu,pd,nrowd)
      use stel_kinds
c
c     dummy jacobian subroutine (not used) for LSODE
c
      real(rprec), dimension(neq) :: y(neq)
      real(rprec), dimension(nrowd,neq) :: pd
      real(rprec) t
      integer neq, ml, mu, nrowd

      end subroutine jac

ext: f

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


! Module providing mesh interaction routines for AGN feedback.
!
! Paul Ricker (2/2013)

module agn_mesh

!==============================================================================

use agn_geometry
use agn_params

    use amrex_fort_module, only : rt => amrex_real
    implicit none

public

! Derived type to hold Cartesian mesh coordinate information

type agn_mesh_data
  integer                       :: Nx, Ny, Nz          ! # of zones in each coordinate
  real(rt), allocatable :: x(:), y(:), z(:)    ! zone center arrays
  real(rt), allocatable :: dx(:), dy(:), dz(:) ! zone spacing arrays
  real(rt)              :: xmin, xmax, ymin, ymax, zmin, zmax
                                           ! min, max values for each coord.
end type agn_mesh_data

contains

!==============================================================================

! average_over_sphere: Average a mesh quantity over a spherical volume of a
!                      given radius about a given center. The radius is
!                      specified in units of the mean zone spacing.
!
! arguments: A                     Mesh quantity to be averaged
!            m                     Mesh data structure
!            x, y, z               Coordinates of sphere center
!            R                     Radius of sphere

real(rt) function average_over_sphere(A, A_l1, A_l2, A_l3, A_h1, A_h2, A_h3, &
                                              x, y, z, R, dx)

integer         , intent(in) :: A_l1, A_l2, A_l3, A_h1, A_h2, A_h3
real(rt), intent(in) :: A(A_l1:,A_l2:,A_l3:)
real(rt), intent(in) :: x, y, z, R
real(rt), intent(in) :: dx(:)

integer             :: i, j, k, Nr, Nmu, Nphi, is, js, ks
real(rt)    :: d_zone, rs, mus, phis, xs, ys, zs, Asum, dr, dmu, dphi, &
                       sint, sinp, cosp

Asum = 0.

Nr   = max(2, int(2*R))
Nmu  = max(5, int(1+2*R))
Nphi = max(10, int(5*R))

d_zone = (dx(1) + dx(2) + dx(3)) / 3.d0
dr   = R*d_zone / (Nr-1)
dmu  = 2. / (Nmu-1)
dphi = 2.*pi / (Nphi-1)

do k = 1, Nphi
  phis = (k-1)*dphi
  sinp = dsin(phis)
  cosp = dcos(phis)
  do j = 1, Nmu
    mus = (j-1)*dmu - 1.d0
    sint = sqrt(1.d0- mus**2)
    do i = 1, Nr
      rs = (i-1)*dr
      xs = x + rs*sint*cosp
      ys = y + rs*sint*sinp
      zs = z + rs*mus

      is = floor(xs/dx(1))
      js = floor(ys/dx(2))
      ks = floor(zs/dx(3))

      if (is .lt. A_l1 .or. is .gt. A_h1 .or. &
          js .lt. A_l2 .or. js .gt. A_h2 .or. &
          ks .lt. A_l3 .or. ks .gt. A_h3) then
          print *,'ACK -- sphere going outside of array bounds '
          print *,'DX         ',dx(:)
          print *,'(x,y,z) is ',xs,ys,zs
          print *,'(i,j,k) is ',is,js,ks
          print *,'   A_lo is ',A_l1,A_l2,A_l3
          print *,'   A_hi is ',A_h1,A_h2,A_h3
          stop
!     else
!        print *,'ADDING TO ',is,js,ks
      end if

      Asum = Asum + A(is,js,ks)
    enddo
  enddo
enddo

average_over_sphere = Asum / (Nr*Nmu*Nphi)

print*,'SUM ',Asum, Nr, Nmu, Nphi, average_over_sphere

return
end function average_over_sphere

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

! map_sphere_onto_mesh: Map a function onto a (possibly) non-Cartesian region
!                       of a Cartesian mesh.  The routine abstracts out the size
!                       and shape of the geometrical region.
!
! arguments: x, y, z            Characteristic coordinates of region
!            R                  Radius of sphere
!            norm               Weighting function normalization (eg. mass)
!            A                  Mesh quantity to map (eg. density)

subroutine map_sphere_onto_mesh(x, y, z, R, norm, &
                                A, A_l1, A_l2, A_l3, A_h1, A_h2, A_h3, &
                                dx)

real(rt), intent(in   ) :: x, y, z, R, norm
real(rt), intent(in   ) :: dx(:)

integer         , intent(in   ) :: A_l1, A_l2, A_l3, A_h1, A_h2, A_h3
real(rt), intent(inout) :: A(A_l1:,A_l2:,A_l3:)

real(rt)                :: vol_sphere

integer, parameter  :: Nsub = 3
integer             :: i, j, k, ii, jj, kk, region
real(rt)    :: Aavg, f, xx, yy, zz, dxx, dyy, dzz
real(rt)    :: xgl, ygl, zgl, xgr, ygr, zgr

vol_sphere = 4.d0/3.d0*pi*R**3
Aavg = norm / vol_sphere

dxx = dx(1) / Nsub
dyy = dx(2) / Nsub
dzz = dx(3) / Nsub

do k = A_l1, A_h1
  zgl = (dble(k)-0.5d0)*dx(3) - z
  zgr = zgl + dx(3)

  do j = A_l2, A_h2
    ygl = (dble(j)-0.5d0)*dx(2) - y
    ygr = ygl + dx(2)

    do i = A_l3, A_h3
      xgl = (dble(i)-0.5d0)*dx(1) - x
      xgr = xgl + dx(1)

      region = loc_sphere(xgl, xgr, ygl, ygr, zgl, zgr, R)

      if (region /= OUTSIDE) then
        f = 0.
        do kk = 1, Nsub
          zz = zgl + (kk-0.5)*dzz
          do jj = 1, Nsub
            yy = ygl + (jj-0.5)*dyy
            do ii = 1, Nsub
              xx = xgl + (i-0.5)*dxx
              if (loc_sphere(xx, xx, yy, yy, zz, zz, R) /= OUTSIDE) &
                f = f + Aavg
            enddo
          enddo
        enddo
        A(i,j,k) = A(i,j,k) + f/Nsub**3
      endif

    enddo
  enddo
enddo

return
end subroutine map_sphere_onto_mesh

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

! map_cylinder_onto_mesh: Map a cylinder onto a (possibly) non-Cartesian region
!                         of a Cartesian mesh.  The routine abstracts out the size
!                         and shape of the geometrical region.
!
! arguments: x, y, z            Characteristic coordinates of region
!            weight             Weighting function
!            wparams            List of weight function parameters (passed)
!            norm               Weighting function normalization (eg. mass)
!            location           Function that determines where a point is with
!                               respect to mapping region
!            volume             Function that returns volume of mapping region
!            A                  Mesh quantity to map (eg. density)
!            m                  Mesh data object

subroutine map_cylinder_onto_mesh(x, y, z, R, height, nxc, nyc, nzc, norm, &
                                  A, A_l1, A_l2, A_l3, A_h1, A_h2, A_h3, &
                                  weight, wparams, dx)

real(rt)   , intent(in   ) :: x, y, z, norm
real(rt)   , intent(in   ) :: R, height, nxc, nyc, nzc
real(rt)   , intent(in   ) :: wparams(:)
integer            , intent(in   ) :: A_l1, A_l2, A_l3, A_h1, A_h2, A_h3
real(rt)   , intent(inout) :: A(A_l1:,A_l2:,A_l3:)

real(rt)                   :: weight, vol_cylinder
real(rt)                   :: dx(:)

integer, parameter  :: Nsub = 3
integer             :: i, j, k, ii, jj, kk, region
real(rt)    :: Aavg, f, xx, yy, zz, dxx, dyy, dzz
real(rt)    :: xgl, ygl, zgl, xgr, ygr, zgr

vol_cylinder = pi * R**2 * height
Aavg = norm / vol_cylinder

dxx = dx(1) / Nsub
dyy = dx(2) / Nsub
dzz = dx(3) / Nsub

do k = A_l3, A_h3
  zgl = (dble(k)-0.5d0)*dx(3) - z
  zgr = zgl + dx(3)
  do j = A_l2, A_h2
    ygl = (dble(j)-0.5d0)*dx(2) - y
    ygr = ygl + dx(2)
    do i = A_l1, A_h1
      xgl = (dble(i)-0.5d0)*dx(1) - x
      xgr = xgl + dx(1)

      region = loc_cylinder(xgl, xgr, ygl, ygr, zgl, zgr, R, height, nxc, nyc, nzc)

      if (region /= OUTSIDE) then
        f = 0.
        do kk = 1, Nsub
          zz = zgl + (kk-0.5)*dzz
          do jj = 1, Nsub
            yy = ygl + (jj-0.5)*dyy
            do ii = 1, Nsub
              xx = xgl + (i-0.5)*dxx
              if (loc_cylinder(xx, xx, yy, yy, zz, zz, R, height, nxc, nyc, nzc) /= OUTSIDE) &
                f = f + Aavg*weight(xx, yy, zz, wparams)
            enddo
          enddo
        enddo
        A(i,j,k) = A(i,j,k) + f/Nsub**3
      endif

    enddo
  enddo
enddo

return
end subroutine map_cylinder_onto_mesh

end module agn_mesh


text:
      SUBROUTINE MOLVAL(C,P,NOCC,RHFUHF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INCLUDE 'SIZES'
      DIMENSION C(NORBS,NORBS), P(*)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),NMIDLE(NUMATM),
     1                NLAST(NUMATM), NORBS, NELECS,NALPHA,NBETA,
     2                NCLOSE,NOPEN,NDUMY,FRACT
      DIMENSION VAL(MAXORB)
      DO 40 I=1,NOCC
         SUM=0.D0
         DO 30 JJ=1,NUMAT
            JL=NFIRST(JJ)
            JU=NLAST(JJ)
            DO 30 J=JL,JU
               DO 30 KK=1,NUMAT
                  IF(KK.EQ.JJ) GOTO 20
                  KL=NFIRST(KK)
                  KU=NLAST(KK)
                  DO 10 K=KL,KU
                     L1=MAX(J,K)
                     L2=J+K-L1
                     L=(L1*(L1-1))/2+L2
                     SUM=SUM+C(J,I)*C(K,I)*P(L)
   10             CONTINUE
   20             CONTINUE
   30    CONTINUE
         VAL(I)=SUM*RHFUHF
   40 CONTINUE
      WRITE(6,'(10F8.4)')(VAL(I),I=1,NOCC)
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
! THIS VERSION: GALAHAD 2.1 - 22/03/2007 AT 09:00 GMT.

!-*-*-*-*-*-*-*-*-  G A L A H A D   R U N L A N B _ S I F  *-*-*-*-*-*-*-*-*-*-

!  Nick Gould, Dominique Orban and Philippe Toint, for GALAHAD productions
!  Copyright reserved
!  March 14th 2003

   PROGRAM RUNLANB_SIF
   USE GALAHAD_USELANB_double

!  Main program for the SIF interface to LANCELOT B, an augmented Lagrangian
!  algorithm for nonlinear programming

!  Problem input characteristics

   INTEGER, PARAMETER :: input = 55
   CHARACTER ( LEN = 16 ) :: prbdat = 'OUTSDIF.d'

!  Open the data input file

   OPEN( input, FILE = prbdat, FORM = 'FORMATTED', STATUS = 'OLD'  )
   REWIND input

!  Call the SIF interface

   CALL USE_LANB( input )

!  Close the data input file

   CLOSE( input  )
   STOP

!  End of RUNLANB_SIF

   END PROGRAM RUNLANB_SIF

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
  subroutine cross_section_Rayleigh(   &
    nsmax,       wavelength_ns,  relN2,    &
    relO2,       relAr,          eps_N2,   &
    eps_O2,      eps_Ar,         csray_N2, &
    csray_O2,    csray_Ar)

!-----------------------------------------------------------
! Subroutine:
!    csray
!
! Purpose:
!    Calculate Rayleigh cross section for
!    N2, O2, Ar, and air
!------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision

! nr of wavelengths
  integer,intent(in):: &
     nsmax
! wavelengths around incident light (ns_elas) [nm]
  real(kind=dp):: &
    wavelength_ns(0:nsmax)

  real(kind=dp),intent(in):: &
     relN2, & !N2 volume mixing ratio
     relO2, & !O2 volume mixing ratio
     relAr    !Ar volume mixing ratio

  real(kind=dp),dimension(0:nsmax):: &
    rn_N2,      & !refraction index of N2
    rn_O2,      & !refraction index of O2
    rn_Ar         !refraction index of Ar

  real(kind=dp):: &
    Wcm,        & !wavelengthlength [cm]
    WNmu,       & !1/wavelengthlength [um-1]
    frac_rn_N2, &
    frac_rn_O2, &
    frac_rn_Ar, &
    rn_N2m1,    &
    rn_O2m1,    &
    rn_Arm1,    &
    fk_N2,      & !King factor
    fk_O2,      &
    fk_Ar,      &
    dens

! parametrization of refraction index of N2 and O2 ( Bates 1984 )
  real(kind=dp):: &
     a_N2(3),   &
     b_N2(3),   &
     a_O2(4),   &
     b_O2(4)
  data &
     a_N2 / 6855.200d0,  5989.242d0,  6998.749d0             /, &
     b_N2 / 3243157.0d0, 3363266.3d0, 3233582.0d0            /, &
     a_O2 / 21351.1d0,   20564.8d0,   22120.4d0,  23796.7d0  /, &
     b_O2 / 218567.0d0,  248089.9d0,  203187.6d0, 168988.4d0 /

  real(kind=dp),parameter::    &
     pi       = 3.141592653589793115997963468544d0, &
     temp_std = 273.15d0,    & !standard temperature [K]
     Dair     = 2.54743D+19    !molar nrdensity air [part. cm-3]

  integer::ni

  real(kind=dp),dimension(0:nsmax):: &
     csray_air    !Rayleigh cross sections for air

  real(kind=dp),dimension(0:nsmax),intent(out):: &
     eps_N2,   &  !epsilon for N2
     eps_O2,   &  !epsilon for O2
     eps_AR,   &  !epsilon for Ar = 0
     csray_O2, &  !Rayleigh cross sections for O2 molecules
     csray_N2, &  !Rayleigh cross sections for N2 molecules
     csray_Ar     !Rayleigh cross sections for Ar atoms

  integer::ns

  do ns=0,nsmax

     WNmu = (1.d+07/wavelength_ns(ns)) * 1.d-4

!    King correction factor
     fk_N2 = 1.034d0 + 3.17d-4 * WNmu**2
     fk_O2 = 1.096d0 + 1.385d-03*WNmu**2 + 1.448d-4*WNmu**4
     fk_Ar = 1.d0

!    anisotropic polarizability factor eps
     eps_N2(ns) = 9.d0/2.d0 * (fk_N2-1.d0)
     eps_O2(ns) = 9.d0/2.d0 * (fk_O2-1.d0)
     eps_Ar(ns) = 0.d0

!    refraction index of N2  (n(N2)-1) [Bates1984]
     ni = 1
     if(wavelength_ns(ns).lT.468.d0) ni = 2
     if(wavelength_ns(ns).lT.254.d0) ni = 3
     rn_N2m1  = (a_N2(ni) + b_N2(ni)/(144.d0-WNmu**2))*1.d-08
     rn_N2(ns) = 1.d0 + rn_N2m1

!    refraction index of O2  (n(O2)-1) [Bates1984]
     ni = 1
     if(wavelength_ns(ns).lT.546.d0)  ni = 2
     if(wavelength_ns(ns).lT.288.d0)  ni = 3
     if(wavelength_ns(ns).lT.221.d0)  ni = 4
     rn_O2m1  = (a_O2(ni) + b_O2(ni)/(40.90d0-WNmu**2))*1.d-08
     rn_O2(ns) = 1.d0 + rn_O2m1

!    refraction index of Ar (n(Ar)-1)  [Peck and Fisher1964]
     rn_Arm1  = (678.6711d0 + 301829.43d0/(144.d0-WNmu**2))*1.d-07
     rn_Ar(ns) = 1.d0 + rn_Arm1

!    (n**2-1)/(n**2+2)
     frac_rn_N2 = (rn_N2(ns)**2-1.d0)**2/(rn_N2(ns)**2+2.d0)**2
     frac_rn_O2 = (rn_O2(ns)**2-1.d0)**2/(rn_O2(ns)**2+2.d0)**2
     frac_rn_Ar = (rn_Ar(ns)**2-1.d0)**2/(rn_Ar(ns)**2+2.d0)**2

!    Dair is given for 15 degree Celius
!    scale to standard temp T = 273.15 K because of refractive index

     dens = (288.16d0/temp_std) * Dair

     Wcm = wavelength_ns(ns) * 1.d-7 !change units nm -> cm

     csray_N2(ns) = 24.d0*pi**3/(Wcm**4*dens**2)*frac_rn_N2*fk_N2
     csray_O2(ns) = 24.d0*pi**3/(Wcm**4*dens**2)*frac_rn_O2*fk_O2
     csray_Ar(ns) = 24.d0*pi**3/(Wcm**4*dens**2)*frac_rn_Ar*fk_Ar

     csray_air(ns) = relN2 * csray_N2(ns) + &
                     relO2 * csray_O2(ns) + &
                     relAr * csray_Ar(ns)

  enddo

  return

end subroutine cross_section_Rayleigh

max_stars_repo_name: rutgervandeelen/das

max_issues_repo_name: rutgervandeelen/das

max_forks_repo_name: rutgervandeelen/das


! --- Reference -------------------------------------
! Physical Review B 51, 14057 (1995)
! Ari P. Seitsonen, M. J. Puska, and R. M. Nieminen
!----------------------------------------------------
MODULE cgpc_seitsonen_module

  use localpot_module, only: Vloc
  use kinetic_module, only: op_kinetic

  implicit none

  PRIVATE
  PUBLIC :: init_cgpc_seitsonen, cgpc_seitsonen

  integer :: ML_0, ML_1
  integer :: comm_grid
  real(8) :: dV
  real(8) :: A = 1.d0

CONTAINS

  SUBROUTINE init_cgpc_seitsonen(n1,n2,comm_in,dV_in)
    implicit none
    integer,intent(IN) :: n1,n2,comm_in
    real(8),intent(IN) :: dV_in
    dV = dV_in
    comm_grid = comm_in
    ML_0 = n1
    ML_1 = n2
  END SUBROUTINE init_cgpc_seitsonen

  SUBROUTINE cgpc_seitsonen(k,s,mm,nn,E,wf,Pgk)
    implicit none
    integer,intent(IN) :: k,s,mm,nn
    real(8),intent(IN) :: E(nn)
#ifdef _DRSDFT_
    real(8),intent(IN)    :: wf(mm,nn)
    real(8),intent(INOUT) :: Pgk(mm,nn)
    real(8),allocatable :: Hwf(:,:)
    real(8),parameter :: zero=0.0d0
#else
    complex(8),intent(IN)    :: wf(mm,nn)
    complex(8),intent(INOUT) :: Pgk(mm,nn)
    complex(8),allocatable :: Hwf(:,:)
    complex(8),parameter :: zero=(0.0d0,0.0d0)
#endif
    integer :: i,n,ierr
    real(8) :: Kr,x
    real(8),allocatable :: T0(:),T1(:)
    include 'mpif.h'
    real(8),parameter :: c1=2.7d1, c2=1.8d1, c3=1.2d1, c4=8.0d0, c5=1.6d1

    allocate( T0(nn) ) ; T0=0.0d0
    allocate( T1(nn) ) ; T1=0.0d0
    allocate( Hwf(mm,nn) ) ; Hwf=zero

    call op_kinetic( wf, Hwf, k )

    do n=1,nn
#ifdef _DRSDFT_
       T0(n) = sum( wf(:,n)*Hwf(:,n) )*dV
#else
       T0(n) = sum( conjg(wf(:,n))*Hwf(:,n) )*dV
#endif
    end do

    deallocate( Hwf )

    call MPI_ALLREDUCE(T0,T1,nn,MPI_REAL8,MPI_SUM,comm_grid,ierr)

    T1(:) = 1.0d0/T1(:)

    do n=1,nn
       do i=1,mm
          x  = A*abs( E(n) - Vloc(i-1+ML_0,s) )*T1(n)
          Kr = ( c1 + c2*x + c3*x**2 + c4*x**3 ) &
              /( c1 + c2*x + c3*x**2 + c4*x**3 + c5*x**4 )
          Pgk(i,n) = Kr*Pgk(i,n)
       end do
    end do

    deallocate( T1 )
    deallocate( T0 )

  END SUBROUTINE cgpc_seitsonen

END MODULE cgpc_seitsonen_module


MODULE arguments
USE mod_mpg
CONTAINS
!--------------------------------------------------------------------
SUBROUTINE get_option(mpg, fn_in)
IMPLICIT NONE
TYPE (mpeg_parameters), INTENT(   OUT) :: mpg
CHARACTER (LEN = *)   , INTENT(IN OUT) :: fn_in
mpg%ibit_rate = 9       ! bit_rate 9 = 128kbps
fn_in         = 'fort'  ! input "fn_in.wav"/ output "fn_in.mp3"
icut          = 25      ! cut subband(icut + 1:)
!qms_stereo    = .TRUE.  ! Mid-Side stereo modo on / off
!q_alias       = .TRUE.  ! take anti-alias in Subband 0-1 in mixed mode
!q_mask        = .TRUE.  ! ATH masking on / off
!q_emphasis    = .TRUE.  ! pre-emphasis
!q_sm          = .TRUE.  ! use both of short & mixed
!q_old_ath     = .FALSE. ! use old ATH function
!q_vbr         = .FALSE. ! VBR mode
!q_rio500      = .FALSE. ! avoid RIO500 VBR bug
!q_info        = .TRUE.  ! print info at the end
!mblock_type_param = 20 ! default short block! long = 0, short = 20, mixed = 21
RETURN
END SUBROUTINE get_option
!--------------------------------------------------------------------
END MODULE arguments
!=====================================================================
MODULE raw_io
IMPLICIT NONE
PRIVATE
PUBLIC :: open_wav_file, close_wav_file, open_mpg_file, close_mpg_file, &
          rd_move, rd_word32, rd_word16, rd_int32, rd_int16, wr_cha8
INTEGER, SAVE :: ir, iw, ir_record = 0, iw_record = 0
LOGICAL, SAVE :: qlittle_endian
CONTAINS
!---------------------------------------------------------------------
SUBROUTINE abort(text)
IMPLICIT NONE
CHARACTER (LEN = *), INTENT(IN) :: text
WRITE(*, *) 'Abort:: ', text
STOP
END SUBROUTINE abort
!--------------------------------------------------------------------
SUBROUTINE check_endian()
IMPLICIT NONE
CHARACTER(LEN = 4) :: ctmp
INTEGER            :: itmp
ctmp = 'ABCD'
itmp = TRANSFER(ctmp, itmp)
IF      ( CHAR( MOD(itmp, 256) ) == 'A' ) THEN
 qlittle_endian = .TRUE.
ELSE IF ( CHAR( MOD(itmp, 256) ) == 'D' ) THEN
 qlittle_endian = .FALSE.
ELSE
 CALL abort('endian indeterminable')
END IF
RETURN
END SUBROUTINE check_endian
!------------------------------------------------------------------
SUBROUTINE open_wav_file(iread, fname)
IMPLICIT NONE
INTEGER            , INTENT(IN    ) :: iread
CHARACTER (LEN = *), INTENT(IN    ) :: fname
INTEGER :: io
ir = iread
CALL check_endian()
OPEN(ir, FILE = fname, STATUS = 'OLD', IOSTAT = io, ACCESS = 'DIRECT', RECL = 1)  ! RECL = 1byte assumed
IF (io /= 0) THEN
 WRITE(*, '(a, i3, a, i3, 2a)' ) ' I/O error ', io, ' occuerred. file =', iread, ' file name ', fname
 CALL abort('Check input file! Suggestion: Is file name correct?')
END IF
RETURN
END SUBROUTINE open_wav_file
!------------------------------------------------------------------
SUBROUTINE close_wav_file
IMPLICIT NONE
CLOSE(ir)
RETURN
END SUBROUTINE close_wav_file
!------------------------------------------------------------------
SUBROUTINE open_mpg_file(iwrite, fname)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iwrite
CHARACTER (LEN = *) :: fname
INTEGER :: io
iw = iwrite
OPEN(iw, FILE = fname, STATUS = 'REPLACE', IOSTAT = io, ACCESS = 'DIRECT', RECL = 1)
IF (io /= 0) THEN
 WRITE(*, '(a, i3, a, i3, 2a)') ' I/O error ', io, ' occuerred. file =', iw, ' file name ', fname
 CALL abort('Check output file! Suggestion: File may be in use.')
END IF
RETURN
END SUBROUTINE open_mpg_file
!-------------------------------------------------------------------
SUBROUTINE close_mpg_file
IMPLICIT NONE
CLOSE(iw)
RETURN
END SUBROUTINE close_mpg_file
!------------------------------------------------------------------
SUBROUTINE rd_move(n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
ir_record = ir_record + n
RETURN
END SUBROUTINE rd_move
!------------------------------------------------------------------
FUNCTION rd_word32() RESULT(res)
IMPLICIT NONE
CHARACTER (LEN = 4) :: res
INTEGER             :: i, io
DO i = 1, 4
 ir_record = ir_record + 1
 READ(ir, REC = ir_record, IOSTAT = io) res(i:i)
END DO
RETURN
END FUNCTION rd_word32
!---------------------------------------------------------------------
FUNCTION rd_word16() RESULT(res)
IMPLICIT NONE
CHARACTER (LEN = 4) :: res
INTEGER             :: i, io
DO i = 1, 2
 ir_record = ir_record + 1
 READ(ir, REC = ir_record, IOSTAT = io) res(i:i)
END DO
RETURN
END FUNCTION rd_word16
!---------------------------------------------------------------------
FUNCTION rd_int32() RESULT(ires)
IMPLICIT NONE
CHARACTER (LEN = 4) :: tmp
INTEGER :: i, io, ires
IF (qlittle_endian) THEN
 DO i = 1, 4
  ir_record = ir_record + 1
  READ(ir, REC = ir_record, IOSTAT = io) tmp(i:i)
 END DO
ELSE
 DO i = 4, 1, -1
  ir_record = ir_record + 1
  READ(ir, REC = ir_record, IOSTAT = io) tmp(i:i)
 END DO
END IF
ires = TRANSFER(tmp, ires)
RETURN
END FUNCTION rd_int32
!---------------------------------------------------------------------
FUNCTION rd_int16() RESULT(ires)
IMPLICIT NONE
CHARACTER (LEN = 4) :: tmp
INTEGER :: i, io, ires
IF (qlittle_endian) THEN
 DO i = 3, 4
  ir_record = ir_record + 1
  READ(ir, REC = ir_record, IOSTAT = io) tmp(i:i)
 END DO
ELSE
 DO i = 2, 1, -1
  ir_record = ir_record + 1
  READ(ir, REC = ir_record, IOSTAT = io) tmp(i:i)
 END DO
END IF
ires = TRANSFER(tmp, ires) / 2**16
RETURN
END FUNCTION rd_int16
!---------------------------------------------------------------------
SUBROUTINE wr_cha8(n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
CHARACTER (LEN = 4) :: cha
INTEGER :: io
cha = TRANSFER(n, cha)
iw_record = iw_record + 1
IF (qlittle_endian) THEN
 WRITE(iw, REC = iw_record, IOSTAT = io) cha(1:1)
ELSE
 WRITE(iw, REC = iw_record, IOSTAT = io) cha(4:4)
END IF
RETURN
END SUBROUTINE wr_cha8
!---------------------------------------------------------------------
END MODULE raw_io
!==================================================================================
MODULE wav_io
USE mod_mpg
USE raw_io
IMPLICIT NONE
PRIVATE
PUBLIC :: riff_chunk, fmt_chunk, fact_chunk, data_chunk              ! type
PUBLIC :: check_riff_chunk, open_wav_file, close_wav_file, read_pcm_1frame, read_pcm0 ! subroutine
!
TYPE :: fmt_chunk
 CHARACTER (LEN = 4):: chunk_id  != 'fmt '
 INTEGER :: ichunk_size
 INTEGER :: iformat_type, ichannels
 INTEGER :: isamples_per_sec
 INTEGER :: ibytes_per_sec
 INTEGER :: iblock_size, ibits_per_sample
END TYPE fmt_chunk
!
TYPE :: fact_chunk
 CHARACTER (LEN = 4) :: chunk_id  != 'fact'
 INTEGER :: ichunk_size !( bytes )
!! same unknown data
END TYPE fact_chunk
!
TYPE :: data_chunk
 CHARACTER (LEN = 4) :: chunk_id  != 'data'
 INTEGER :: ichunk_size !( bytes )
!! PCM data follows INTEGER (KIND = 2) :: idata(isize / 2)
END TYPE data_chunk
!
TYPE :: riff_chunk
 CHARACTER (LEN = 4) :: chunk_id  != 'RIFF'
 INTEGER :: ichunk_size
 CHARACTER (LEN = 4) :: format_type ! = 'WAVE'
 TYPE ( fmt_chunk) :: fmt
 TYPE (fact_chunk) :: fct
 TYPE (data_chunk) :: dat
END TYPE riff_chunk
!
CONTAINS
!--------------------------------------------------------------------
SUBROUTINE init_chunk_names(riff)
IMPLICIT NONE
TYPE (riff_chunk), INTENT(OUT) :: riff
riff%chunk_id     = 'RIFF'
riff%format_type  = 'WAVE'
riff%fmt%chunk_id = 'fmt '
riff%fct%chunk_id = 'fact'
riff%dat%chunk_id = 'data'
RETURN
END SUBROUTINE init_chunk_names
!------------------------------------------------------------------
SUBROUTINE check_riff_chunk(riff)
IMPLICIT NONE
TYPE (riff_chunk), INTENT(IN OUT) :: riff
CALL init_chunk_names(riff)
IF ( rd_word32() == riff%chunk_id ) THEN    ! 'RIFF'?
 WRITE(*, '(a)', ADVANCE = 'NO') ' MS RIFF '
ELSE
 CALL abort('This is not MS-RIFF file!')
END IF
riff%ichunk_size = rd_int32()
IF ( rd_word32() == riff%format_type ) THEN ! 'WAVE'?
 WRITE(*, '(a)', ADVANCE = 'NO') 'WAV audio '
ELSE
 WRITE(*, *)
 CALL abort('This is not WAV file!')
END IF
CALL check_fmt_chunk(riff%fmt)
CALL check_fct_chunk(riff%fct)
CALL check_dat_chunk(riff%dat)
RETURN
END SUBROUTINE check_riff_chunk
!------------------------------------------------------------------
SUBROUTINE check_fmt_chunk(fmt)
IMPLICIT NONE
TYPE (fmt_chunk), INTENT(IN OUT) :: fmt
IF ( rd_word32() /= fmt%chunk_id ) CALL abort('Cannot find format chunk!')
fmt%ichunk_size      = rd_int32()
fmt%iformat_type     = rd_int16()
fmt%ichannels        = rd_int16()
fmt%isamples_per_sec = rd_int32()
fmt%ibytes_per_sec   = rd_int32()
fmt%iblock_size      = rd_int16()
fmt%ibits_per_sample = rd_int16()
IF ( fmt%iformat_type     /=  1) CALL abort('Unknown WAVE format!') !linear PCM
IF ( fmt%ibits_per_sample /= 16) CALL abort('Not 16bit data!')
SELECT CASE ( fmt%ichannels )
 CASE (1)
  WRITE(*, '(a, i3, a, i6, a)', ADVANCE = 'NO') &
   'Monoral', fmt%ibits_per_sample, 'bit Sampling rate', fmt%isamples_per_sec, 'Hz '
 CASE (2)
  WRITE(*, '(a, i3, a, i6, a)', ADVANCE = 'NO') &
   'Stereo' , fmt%ibits_per_sample, 'bit Sampling rate', fmt%isamples_per_sec, 'Hz '
 CASE DEFAULT
  WRITE(*, '(a, i1)') ' Number of wave channels is ', fmt%ichannels
  CALL abort('Wave channel must be 1 or 2!')
END SELECT
RETURN
END SUBROUTINE check_fmt_chunk
!------------------------------------------------------------------
SUBROUTINE check_fct_chunk(fct)
IMPLICIT NONE
TYPE (fact_chunk), INTENT(IN OUT) :: fct
IF ( rd_word32() == fct%chunk_id ) THEN
 fct%ichunk_size = rd_int32()
 CALL rd_move(fct%ichunk_size) ! skip to next chunk
ELSE
 CALL rd_move(-4) ! back to beginning of this chunk
END IF
RETURN
END SUBROUTINE check_fct_chunk
!------------------------------------------------------------------
SUBROUTINE check_dat_chunk(dat)
IMPLICIT NONE
TYPE (data_chunk), INTENT(IN OUT) :: dat
IF ( rd_word32() /= dat%chunk_id ) CALL abort('Cannot find data chunk!')
dat%ichunk_size = rd_int32()
RETURN
END SUBROUTINE check_dat_chunk
!------------------------------------------------------------------
SUBROUTINE wav_read(pcm)
IMPLICIT NONE
REAL (KIND = 8), INTENT(OUT) :: pcm(:, :)
REAL (KIND = 8), PARAMETER   :: denom = 32768.0d0 !32768 = 2^15
INTEGER        , PARAMETER   :: maxbuff = 1152 * 2
INTEGER :: i, nchannel, ndat
ndat     = SIZE(pcm, 1)
nchannel = SIZE(pcm, 2)
IF (ndat * nchannel > maxbuff) CALL abort('check maxbuff: subroutine wav_get')
SELECT CASE (nchannel)
 CASE (1) !mono
  DO i = 1, ndat
   pcm(i, 1) = REAL( rd_int16(), KIND = 8) / denom
   pcm(i, 2) = 0.0d0
  END DO
 CASE (2) !stereo
  DO i = 1, ndat
   pcm(i, 1) = REAL( rd_int16(), KIND = 8) / denom
   pcm(i, 2) = REAL( rd_int16(), KIND = 8) / denom
  END DO
 CASE DEFAULT
  CALL abort('ichannel must be 1 or 2: subroutine wav_get')
END SELECT
RETURN
END SUBROUTINE wav_read
!------------------------------------------------------------------
SUBROUTINE read_pcm_1frame(pcm)
IMPLICIT NONE
REAL (KIND = 8), INTENT(OUT) :: pcm(:, :)
pcm = EOSHIFT(pcm, 1152, 0.0d0, 1)
CALL wav_read(pcm(481:1632, :))
RETURN
END SUBROUTINE read_pcm_1frame
!------------------------------------------------------------------
SUBROUTINE read_pcm0(pcm)
IMPLICIT NONE
REAL (KIND = 8), INTENT(OUT) :: pcm(:, :)
CALL wav_read(pcm(1153:1632, :))
RETURN
END SUBROUTINE read_pcm0
!------------------------------------------------------------------
END MODULE wav_io
!======================================================================
MODULE bit_io
USE raw_io
USE mod_mpg
USE mod_crc
USE mod_layer3
USE mod_huffman
IMPLICIT NONE
INTEGER,SAVE :: ip
CHARACTER (LEN = 10000) :: bit_string
CONTAINS
!---------------------------------------------------------------------
SUBROUTINE clear_bit_buff()
IMPLICIT NONE
ip = 1
bit_string = REPEAT(' ', LEN(bit_string))
RETURN
END SUBROUTINE clear_bit_buff
!-------------------------------------------------------------------
SUBROUTINE put_bits(n, inp)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n, inp
INTEGER :: i, m
IF (n > 32) CALL abort('out of range: n must be < 32: subroutine put_bits')
DO i = 1, n
 m = 2**(n - i)
 IF (MOD(inp / m, 2) == 1) THEN
  bit_string(ip:ip) = '1'
 ELSE
  bit_string(ip:ip) = '0'
 END IF
 ip = ip + 1
END DO
RETURN
END SUBROUTINE put_bits
!-------------------------------------------------------------------
SUBROUTINE put_bits_dim(n, inp)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n, inp(:)
INTEGER :: i
DO i = 1, SIZE(inp)
 CALL put_bits(n, inp(i))
END DO
RETURN
END SUBROUTINE put_bits_dim
!-------------------------------------------------------------------
SUBROUTINE put_bits_dim2(n, inp)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n, inp(:, :)
INTEGER :: i, j
DO i = 1, SIZE(inp, 1)
 DO j = 1, SIZE(inp, 2)
  CALL put_bits(n, inp(i, j))
 END DO
END DO
RETURN
END SUBROUTINE put_bits_dim2
!-------------------------------------------------------------------
SUBROUTINE put_bits_c(str)
IMPLICIT NONE
CHARACTER (LEN = *) :: str
INTEGER :: i
DO i = 1, LEN_TRIM(str)
 IF (str(i:i) /= '0' .AND. str(i:i) /= '1') CALL abort('invalid string: subroutine put_bit_c')
 bit_string(ip:ip) = str(i:i)
 ip = ip + 1
END DO
RETURN
END SUBROUTINE put_bits_c
!-------------------------------------------------------------------
SUBROUTINE write_bits_1frame(n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
INTEGER :: i, j, ipos, m
IF (MOD(n, 8) /= 0) CALL abort('input error: n must be multiples of 8: SUBROUTINE write_bits_1frame')
ipos = 0
DO i = 1, n, 8
 m = 0
 DO j = 1, 8
  ipos = ipos + 1
  IF (ipos > LEN(bit_string)) EXIT
  IF (bit_string(ipos:ipos) == '1') m = m + 2**(8 - j)
 END DO
 CALL wr_cha8(m)
END DO
RETURN
END SUBROUTINE write_bits_1frame
!-------------------------------------------------------------------
END MODULE bit_io


text:
c***********************************************************************
c	                    COORD_VH2B_TP.F
c***********************************************************************
*     PURPOSE: Converts test part from Heliocentric to Barycentric coords.
*              Velocity only
*     ARGUMENTS:  Input is
*                              ntp ==> number of test part (<= NTPMAX)
*                                              (integer)
*		     vxht,vyht,vzht ==> heliocentric particle velocities
*                                             (real array)
*		     vxhs,vyhs,vzhs ==> bary vel of the Sun
*                                          (real scalar)
*                 Returned are
*                    vxbt,vybt,vzbt ==> bary. particle velocities
*                                            (real array)
*
*     Authors:  Hal Levison
*     ALGORITHM: Obvious
*     WRITTEN:  11/14/96
*     REVISIONS:

	subroutine coord_vh2b_tp(ntp,vxht,vyht,vzht,vxs,vys,vzs,
     &      vxbt,vybt,vzbt)

      include '../swift.inc'

c...  Inputs:
	integer ntp
	real*8 vxht(NTPMAX),vyht(NTPMAX),vzht(NTPMAX)
        real*8 vxs,vys,vzs

c...  Outputs:
	real*8 vxbt(NTPMAX),vybt(NTPMAX),vzbt(NTPMAX)

c...  Internals:
	integer i

c----
c...  Executable code
	do i=1,ntp
	  vxbt(i) = vxht(i) + vxs
	  vybt(i) = vyht(i) + vys
	  vzbt(i) = vzht(i) + vzs
	enddo

	return
	end     ! coord_vh2b_tp
c--------------------------------------------------------------------------

max_stars_repo_name: Simske/exostriker


text:
      subroutine rpotph (io, nhead0, head0, lhead0,
     1             nat, nph, nfr, ihole, gamach, iafolp, intclc,
     1             ixc, vr0, vi0, rs0, iphat, rat, iatph, ifrph,
     1             xnatph, novr,
     2             iphovr, nnovr, rovr, folp, ion, iz, iprint,
     2             ixanes, nemax, xkmin, xkmax, potlbl)
      implicit double precision (a-h, o-z)

c     Notes:
c        nat   number of atoms in problem
c        nph   number of unique potentials
c        nfr   number of unique free atoms
c        ihole hole code of absorbing atom
c        iph=0 for central atom
c        ifr=0 for central atom
c        xkmin, xkmax  min and max energy mesh points to consider

      include 'dim.h'

      character*(*) head0(nhead0)
      dimension lhead0(nhead0)

c     End of line comments removed -- see include file arrays.h for
c     comments.
c     Specific atom input data
      dimension iphat(natx)
      dimension rat(3,natx)

c     Unique potential input data
      dimension iatph(0:nphx)
      dimension ifrph(0:nphx)
      dimension xnatph(0:nphx)
      character*6  potlbl(0:nphx)

      dimension folp(0:nphx)
      dimension novr(0:nphx)
      dimension iphovr(novrx,0:nphx)
      dimension nnovr(novrx,0:nphx)
      dimension rovr(novrx,0:nphx)

c     Free atom data
      dimension ion(0:nfrx)
      dimension iz(0:nfrx)

c     read and save header from old file, has carriage control char
      head0(1) = ' '
      call rdhead (io, nhead0, head0, lhead0)
      read(io,*) ihole, gamach, iprint, iafolp, intclc
      read(io,*) ixc, vr0, vi0, rs0
      read(io,*) ixanes, nemax, xkmin, xkmax
      read(io,*) nfr
      do 710  ifr = 0, nfr
         read(io,*)  index, iz(ifr), ion(ifr)
  710 continue
      read(io,*) nat
      do 720  iat = 1, nat
         read(io,*) index, iphat(iat), (rat(j,iat),j=1,3)
  720 continue
      read(io,*) nph
      do 740  iph = 0, nph
         read(io,*) index, iatph(iph), ifrph(iph), xnatph(iph),
     1                folp(iph), novr(iph)
         read(io,*) potlbl(iph)
         do 730  iovr = 1, novr(iph)
            read(io,*) iphovr(iovr,iph), nnovr(iovr,iph),
     1                   rovr(iovr,iph)
  730    continue
  740 continue

      return
      end

max_stars_repo_name: keechul/ifeffit

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: bruceravel/ifeffit

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: bruceravel/ifeffit

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
!*********************************************************
!>
! A number of constants used in math and physics
!

module nfConstants
  implicit none

  integer, parameter :: sp = selected_real_kind(6, 37)
  integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: qp = selected_real_kind(33, 4931)

  ! Math constants

  ! Pi
  real(dp), parameter :: pi_ = 2.d0*acos(0.d0)

  ! Eulers number
  real(dp), parameter :: e_ = exp(1.d0)

  ! Euler-Mascheroni constant
  real(dp), parameter :: em_ = 0.57721566490153286d0

  ! Complex i
  complex(dp), parameter :: i_ = (0,1)

  ! Physics constants

  ! Speed of light in m/s
  real(dp), parameter :: c_ = 2.99792458d8
  ! Speed of light in cm/s (cgs units)
  real(dp), parameter :: c_cgs_ = 2.99792458d10

  ! Planck constant in J*s
  real(dp), parameter :: h_ = 6.626070040d-34
  ! Reduced Planck constant, h/2*pi, in J*s
  real(dp), parameter :: hbar_ = 1.054571800d-34
  ! Planck constant in eV*s
  real(dp), parameter :: h_eV_ = 4.135667662d-15
  ! Reduced Planck constant, \(h/2\pi\), in \(eV\) \(s\)
  real(dp), parameter :: hbar_eV_ = 6.582119514d-16
  ! Planck constant in \(erg\) \(s\) (cgs units)
  real(dp), parameter :: h_cgs_ = 6.626070040d-27
  ! Reduced Planck constant, \(h/2\pi\), in \(erg\) \(s\) (cgs units)
  real(dp), parameter :: hbar_cgs_ = 1.054571800d-27
  ! Reduced Planck constant, \(hc/2\pi\), in \(MeV\) \(fm\)
  real(dp), parameter :: hbarc_ = 197.3269788d0

  ! Gravitational constant in m^3 kg^-1 s^-2
  real(dp), parameter :: G_ = 6.6743d-11
  ! Gravitational constant in cm^3 g^-1 s^-2 (dyne cm^2 g^-2) (cgs units)
  real(dp), parameter :: G_cgs_ = 6.6743d-8

  ! Mass of electron in g (cgs units)
  real(dp), parameter :: me_ = 9.10938356d-28
  ! Mass of muon in g (cgs units)
  real(dp), parameter :: mmu_ = 1.883531594d-25
  ! Mass of tau in g (cgs units)
  real(dp), parameter :: mtau_ = 3.16747d-24
  ! Mass of proton in g (cgs units)
  real(dp), parameter :: mp_ = 1.672621898d-24
  ! Mass of neutron in g (cgs units)
  real(dp), parameter :: mn_ = 1.674927471d-24
  ! Mass of alpha particle in g (cgs units)
  real(dp), parameter :: malpha_ = 6.644657230d-24
  ! Atomic mass unit in g (cgs units)
  real(dp), parameter :: amu_ = 1.660539d-24

  ! Bohr radius in cm (cgs units)
  real(dp), parameter :: bohrr_ = 5.2917721067d-9
  ! Fine-structure constant
  real(dp), parameter :: alpha_ = 7.2973525664d-3
  ! Hartree energy in J
  real(dp), parameter :: eHartree_ = 4.359744650d-18
  ! Hartree energy in eV
  real(dp), parameter :: eHartree_eV_ = 27.21138602d0

  ! Avagadro's number
  real(dp), parameter :: na_ = 6.0221367d23
  ! Boltzmann constant in J K^-1
  real(dp), parameter :: k_ = 1.38064852d-23
  ! Faraday constant in C mol^-1
  real(dp), parameter :: F_ = 96485.33289d0
  ! Rydberg constant in m^-1
  real(dp), parameter :: ryd_ = 10973731.568508d0
  ! Stefan-Boltzmann constant in W m^-2 K^-4
  real(dp), parameter :: stefboltz_ = 5.670367d-8
  ! Wien displacement law constant in m K
  real(dp), parameter :: wiend_ = 2.8977729d-3

  ! Electric charge in C
  real(dp), parameter :: ec_ = 1.6021766208d-19
  ! Electric charge in cm^{3/2} g^{1/2} s^-1 (cgs units)
  real(dp), parameter :: ec_cgs_ = 4.8032d-10
  ! Electron volt in J
  real(dp), parameter :: ev_ = 1.6021766208d-19
  ! Electron volt in erg (cgs units)
  real(dp), parameter :: ev_cgs_ = 1.60217733d-12

  ! Electric constant (vacuum permittivity) in F m^-1
  real(dp), parameter :: epsilon0_ = 8.854187817d-12
  ! Magnetic constant (vacuum permeability) in N A^-2
  real(dp), parameter :: mu0_ = 1.2566370614d-6

  ! Solar mass in g
  real(dp), parameter :: msolar_ = 1.988435d33
  ! Solar radius in cm
  real(dp), parameter :: rsolar_ = 6.955d10
  ! Solar luminosity in erg s^-1
  real(dp), parameter :: lumsolar_ = 3.839d33
  ! Astronomical Unit in cm
  real(dp), parameter :: au_ = 1.49597870700d13

end module nfConstants

max_stars_repo_name: joshhooker/numFortran

max_issues_repo_name: joshhooker/numFortran

max_forks_repo_name: joshhooker/numFortran


text:
cvk  1/2008 ------------------------------------------------------------------
cvk  1/2008 Added option to normalize soil moisture
cvk  1/2008 To select this option, change the statement in the subr. header to
cvk  1/2008         PARAMETER (NORMALIZED = 1)
cvk  1/2008 Default not normalized option is
cvk  1/2008         PARAMETER (NORMALIZED = 0)
cvk  1/2008 ------------------------------------------------------------------

C MEMBER FLAND1
C  (from old member FCFLAND1)
C
      SUBROUTINE FLAND1(PXV,EDMND,TA,WE,AESC,SH,DT,SACST,FRZST,SACPAR,
     +                FRZPAR,NSOIL,NUPL,NSAC,IVERS,SURF,GRND,TET,SMC,
     +                SH2O,SACST_PRV,DTFRZ,IDTFRZ,FRZDUP,FRZDBT,FROST,
     +             TSINT, SWINT, SWHINT, DSINT, NDSINT, DSINTW, NDINTW,
     +             NORMALIZE, HRAPX, HRAPY,
     +             ERROR )

cvk  1/2008  Introduced option to generate normalized SM
cvk  1/2008  To select this option, change the next parameter from 0 to 1
cvk  1/2008
cvk Normalized soil moisture content
cc      PARAMETER (NORMALIZE = 1)
cvk Not normalized soil moisture content
cc      PARAMETER (NORMALIZE = 0)

c  DT is in days here
C  DTFRZ IN SEC., IDTFRZ IS # FRZ_STEPS
C.......................................
C     THIS SUBROUTINE EXECUTES THE 'SAC-SMA ' OPERATION FOR ONE TIME
C         PERIOD.
C.......................................
C     SUBROUTINE INITIALLY WRITTEN BY. . .
C            ERIC ANDERSON - HRL     APRIL 1979     VERSION 1
C.......................................

CVK  FROZEN GROUND CHANGES
CVK  UZTWC,UZFWC,LZTWC,LZFSC,LZFPC ARE TOTAL WATER STORAGES
CVK  UZTWH,UZFWH,LZTWH,LZFSH,LZFPH ARE UNFROZEN WATER STORAGES

      PARAMETER (T0 = 273.16)
      REAL SACPAR(*),FRZPAR(*),SACST(*),FRZST(*),SACST_PRV(*)
      real smc(*),sh2o(*)
c SACPAR() is array of original SAC parameters, and FRZPAR() is array
c of frozen ground parameters and calculated constants
c SACST() and FRZST() same for states

c  delited real FGCO(6),ZSOIL(6),TSOIL(8),FGPM(11)
      REAL LZTWM,LZFSM,LZFPM,LZSK,LZPK,LZTWC,LZFSC,LZFPC
      REAL LZTWH,LZFSH,LZFPH

CVK----------------------------------------------------------------
CVK_02  NEW COMMON STATEMENT FOR DESIRED SOIL LAYERS
CVK     THIS VERSION HAS HARD CODED OUTPUT SOIL LAYERS
CVK     LATER ON IT SHOULD BE CHANGED TO MAKE THEM VARIABLE
c      INTEGER NINT/5/,NINTW/5/
      INTEGER NDSINT,NDINTW, NORMALIZE, ERROR
      REAL TSINT(*),SWINT(*),SWHINT(*)
ck      REAL DSINT(10)/0.075,0.15,0.35,0.75,1.5,0.0,0.0,0.,0.,0./
c      REAL DSINT(10)/0.10,0.40,0.6,0.75,1.5,0.0,0.0,0.,0.,0./
ck      REAL DSINTW(10)/0.075,0.15,0.35,0.75,1.5,0.0,0.0,0.,0.,0./
      REAL DSINT(*), DSINTW(*)
c      REAL DSINTW(10)/0.10,0.40,0.6,0.75,1.5,0.0,0.0,0.,0.,0./
      REAL TSTMP(10),DSMOD(10),SWTMP(10),SWHTMP(10)
c      SAVE DSINT,DSINTW,NINT,NINTW
CVK----------------------------------------------------------------

cc      DIMENSION EPDIST(24)
C     COMMON BLOCKS
cc      COMMON/FSMPM1/UZTWM,UZFWM,UZK,PCTIM,ADIMP,RIVA,ZPERC,REXP,LZTWM,
cc     1              LZFSM,LZFPM,LZSK,LZPK,PFREE,SIDE,SAVED,PAREA
CVK
CVK      COMMON/FPMFG1/FGPM(10)
c--      COMMON/FPMFG1/itta,FGPM(15),ivers,ifrze

CVK_02  NEW COMMON BLOCK FOR INTERPOLATED SOIL TEMP AND SOIL PARAMETERS
c--      COMMON/TSLINT/TSINT(10),NINT,SWINT(10),SWHINT(10),NINTW
c--      COMMON/FRDSTFG/SMAX,PSISAT,BRT,SWLT,QUARTZ,STYPE,NUPL,NSAC,
c--     +               RTUZ,RTLZ,DZUP,DZLOW
cc      COMMON/FRZCNST/ FRST_FACT,CKSOIL,ZBOT
c--      COMMON/FRZCNST/ FRST_FACT,ZBOT

CVK
CVK  NEW FG VERSION PARAMETERS & SOIL LAYER DEFINITION:
CVK          FGPM(1) - SOIL TEXTURE CLASS
CVK          FGPM(2) - OPTIONAL, SOIL TEMPERATURE AT THE 3M DEPTH
CVK          FGPM(3) - OPTIONAL, POROSITY OF RESIDUE LAYER
CVK          PAR(18) [if no calb=FGPM(4)] - RUNOFF REDUCTION PARAMETER 1
CVK          PAR(19) [if no calb=FGPM(5)] - RUNOFF REDUCTION PARAMETER 2
CVK          PAR(20) [if no calb=FGPM(6)] - RUNOFF REDUCTION PARAMETER 3
CVK          FGPM(6) - RUNOFF REDUCTION PARAMETER 3 (FOR ERIC'S VERSION ONLY)
CVK          FGPM(7) - NUMBER OF SOIL LAYERS
CVK          FGPM(8)-FGPM(15) - DEPTHS OF SOIL LAYERS (M), NEGATIVE.
CVK                             FIRST LAYER (RESIDUE) DEPTH=-0.03M

c--      COMMON/FSMCO1/UZTWC,UZFWC,LZTWC,LZFSC,LZFPC,ADIMC,FGCO(6),RSUM(7),
c--     1   PPE,PSC,PTA,PWE,PSH,TSOIL(8)

c--      COMMON/FSUMS1/SROT,SIMPVT,SRODT,SROST,SINTFT,SGWFP,SGWFS,SRECHT,
c--     1              SETT,SE1,SE3,SE4,SE5

C
C    ================================= RCS keyword statements ==========
c--      CHARACTER*68     RCSKW1,RCSKW2
c--      DATA             RCSKW1,RCSKW2 /                                 '
c--     .$Source: /fs/hsmb5/hydro/CVS_root/gorms/sac/fland1.f,v $
c--     . $',                                                             '
c--     .$Id: fland1.f,v 3.1 2012-02-02 21:32:13 zcui Exp $
c--     . $' /
C    ===================================================================

CVK  ADDITION FOR FROZEN DEPTH ESTIMATION
c--      SAVE IIPREV,FRSTPREV,FRZPREV
c--      IF(FGCO(1) .EQ. 0.) THEN
c--       IIPREV=0
c--       FRSTPREV=0.
c--       FRZPREV=0.
c--      ENDIF
CVK--------------------------------------

c       write(*,*) ' FRZPAR:',(frzpar(ii),ii=1,6)
c define major parameters from the array
c      WRITE(*,*) 'FLAND1: SACPAR: ', (SACPAR(i), i=1, 17 )
      do i=1,6
       if(sacst(i) .lt. -1.0) then
         write(*,*) 'Passed SAC state#',i,'<-1.',sacst(i)
       endif
      enddo

      ERROR = 0
      UZTWM=SACPAR(1)
      UZFWM=SACPAR(2)
      ADIMP=SACPAR(5)
      LZTWM=SACPAR(9)
      LZFSM=SACPAR(10)
      LZFPM=SACPAR(11)
      PAREA=1.0-SACPAR(4)-ADIMP
      IF(IVERS .NE. 0) THEN
          CKSL=FRZPAR(4)
          SMCDRY = FRZPAR( 9 )
      ELSE
         SMCDRY = 0.0
      ENDIF

c define states from the array
      UZTWC=SACST(1)
      UZFWC=SACST(2)
      LZTWC=SACST(3)
      LZFSC=SACST(4)
      LZFPC=SACST(5)
      ADIMC=SACST(6)

c      WRITE(*,*) 'FLAND1: UZFWC= ', UZFWC
      IF(IVERS .EQ. 0) THEN
CVK  OLD FROZEN GROUND VERSION: KEEP UNFROZEN WATER = TOTAL
C--       RUZICE=UZK
C--       RLZICE=LZSK
C--       RUZPERC=1.0
       UZTWH=UZTWC
       UZFWH=UZFWC
       LZTWH=LZTWC
       LZFSH=LZFSC
       LZFPH=LZFPC
      ELSE
CVK  NEW FROZEN GROUND VERSION: USE ESTIMATED UNFROZEN WATER
CVK  REDEFINE UNFROZEN WATER VARIABLES IF THEY ARE NEGATIVE
C
       UZTWH=FRZST(6)
       IF(UZTWH .LT. 0.) UZTWH=0.
       UZFWH=FRZST(7)
       IF(UZFWH .LT. 0.) UZFWH=0.
       LZTWH=FRZST(8)
       IF(LZTWH .LT. 0.) LZTWH=0.
       LZFSH=FRZST(9)
       IF(LZFSH .LT. 0.) LZFSH=0.
       LZFPH=FRZST(10)
       IF(LZFPH .LT. 0.) LZFPH=0.

CVK  RUZICE & RLZICE ARE REDUCTION OF FREE WATER MOVEMENT
CVK  BASED ON KULIK'S THEORY: Kfrz = Kunfrz/(1+FGPM(4)*ICE)**2
       ZUP=FRZPAR(9+NUPL)
       ZLW=FRZPAR(9+NSAC)
       RUZICE=0.001*FRZPAR(6)*(UZTWC-UZTWH+UZFWC-UZFWH)/
     +        (FRZPAR(10)-ZUP)
       RLZICE=0.001*FRZPAR(7)*(LZTWC-LZTWH+LZFSC-LZFSH)/
     +        (ZUP-ZLW)
       RUZPERC=1.0
       IF(RUZICE .EQ. 0.) THEN
        RUZICE = SACPAR(3)
       ELSE
        RUZPERC=1.0/((1.+CKSL*RUZICE)**2)
        RUZICE=1.-EXP(LOG(1.-SACPAR(3))*RUZPERC)
       ENDIF
       IF(RLZICE .EQ. 0.) THEN
        RLZICE = SACPAR(12)
       ELSE
        RLZICE=1.0/((1.+CKSL*RLZICE)**2)
        RLZICE=1.-EXP(LOG(1.-SACPAR(12))*RLZICE)
       ENDIF
      ENDIF
c      if(uztwc .ne. uztwh) WRITE(*,*) 'ST1=',uztwc,uztwh
c      if(uzfwc .ne. uzfwh) WRITE(*,*) 'ST2=',uzfwc,uzfwh
c      if(lztwc .ne. lztwh) WRITE(*,*) 'ST3=',lztwc,lztwh
c      if(lzfsc .ne. lzfsh) WRITE(*,*) 'ST4=',lzfsc,lzfsh
c      if(lzfpc .ne. lzfph) WRITE(*,*) 'ST5=',lzfpc,lzfph

C.......................................
C     COMPUTE EVAPOTRANSPIRATION LOSS FOR THE TIME INTERVAL.
C        EDMND IS THE ET-DEMAND FOR THE TIME INTERVAL
cc      EDMND=EP*EPDIST(KINT)
cVK ADJUST EDMND FOR EFFECT OF SNOW & FOREST COVER.
c from EX1 OFS subroutine
      EDMND=(1.-(1.0-SACPAR(17))*AESC)*EDMND
C
C     COMPUTE ET FROM UPPER ZONE.
CVK      E1=EDMND*(UZTWC/UZTWM)
CVK  ONLY UNFROZEN WATER CAN BE EVAPORATED
      E1=EDMND*(UZTWH/UZTWM)

      RED=EDMND-E1
C     RED IS RESIDUAL EVAP DEMAND
CVK      UZTWC=UZTWC-E1
      UZTWH=UZTWH-E1

      E2=0.0
CV.K      IF(UZTWC.GE.0.) THEN
      IF(UZTWH.GE.0.) THEN
CV.K    SUBTRACT ET FROM TOTAL WATER STORAGE
       UZTWC=UZTWC-E1
       GO TO 220
      ENDIF

C     E1 CAN NOT EXCEED UZTWC
CV.K      E1=E1+UZTWC
CV.K      UZTWC=0.0
      E1=E1+UZTWH
      UZTWH=0.0
CV.K   REDUCE TOTAL TENSION WATER BY ACTUAL E1
      UZTWC=UZTWC-E1
      IF(UZTWC .LT. 0.0) UZTWC=0.0

      RED=EDMND-E1
CV.K      IF(UZFWC.GE.RED) GO TO 221
      IF(UZFWH.GE.RED) GO TO 221

C     E2 IS EVAP FROM UZFWC.
CV.K      E2=UZFWC
CV.K      UZFWC=0.0
      E2=UZFWH
      UZFWH=0.0
CV.K   REDUCE TOTAL FREE WATER BY ACTUAL E2
      UZFWC=UZFWC-E2
      IF(UZFWC .LT. 0.0) UZFWC=0.0

      RED=RED-E2
      GO TO 225
  221 E2=RED
CVK   SUBTRACT E2 FROM TOTAL & UNFROZEN FREE WATER STORAGES
      UZFWC=UZFWC-E2
      UZFWH=UZFWH-E2
      RED=0.0
  220 IF((UZTWC/UZTWM).GE.(UZFWC/UZFWM)) GO TO 225
C     UPPER ZONE FREE WATER RATIO EXCEEDS UPPER ZONE
C     TENSION WATER RATIO, THUS TRANSFER FREE WATER TO TENSION
      UZRAT=(UZTWC+UZFWC)/(UZTWM+UZFWM)

CV.K  ACCOUNT FOR RATIO OF UNFROZEN WATER ONLY
CV.K  AND ADJUST FOUR SOIL STATES
CV.K      UZTWC=UZTWM*UZRAT
CV.K      UZFWC=UZFWM*UZRAT
      DUZTWC=UZTWM*UZRAT-UZTWC
      IF(DUZTWC .GT. UZFWH) DUZTWC=UZFWH
CV.K  TRANSFERED WATER CAN NOT EXCEED UNFROZEN FREE WATER
      UZTWC=UZTWC+DUZTWC
      UZTWH=UZTWH+DUZTWC
      UZFWC=UZFWC-DUZTWC
      UZFWH=UZFWH-DUZTWC

CV.K  CHECK UNFROZEN WATER STORAGES TOO
  225 IF (UZTWC.LT.0.00001) THEN
       UZTWC=0.0
       UZTWH=0.0
      ENDIF
      IF (UZFWC.LT.0.00001) THEN
       UZFWC=0.0
       UZFWH=0.0
      ENDIF
C
C     COMPUTE ET FROM THE LOWER ZONE.
C     COMPUTE ET FROM LZTWC (E3)
CV.K      E3=RED*(LZTWC/(UZTWM+LZTWM))
CV.K      LZTWC=LZTWC-E3
CV.K      IF(LZTWC.GE.0.0) THEN
CV.K  ONLY UNFROZEN WATER CAN BE EVAPORATED
      E3=RED*(LZTWH/(UZTWM+LZTWM))
      LZTWH=LZTWH-E3
      IF(LZTWH.GE.0.0) THEN
       LZTWC=LZTWC-E3
       GO TO 226
      ENDIF

C     E3 CAN NOT EXCEED LZTWC
CV.K      E3=E3+LZTWC
CV.K      LZTWC=0.0
      E3=E3+LZTWH
      LZTWH=0.0
CV.K   REDUCE TOTAL TENSION WATER BY E3
       LZTWC=LZTWC-E3

  226 RATLZT=LZTWC/LZTWM
      RATLZ=(LZTWC+LZFPC+LZFSC-SACPAR(16))/(LZTWM+LZFPM+LZFSM
     +       -SACPAR(16))
      IF(RATLZT.GE.RATLZ) GO TO 230
C     RESUPPLY LOWER ZONE TENSION WATER FROM LOWER
C     ZONE FREE WATER IF MORE WATER AVAILABLE THERE.
      DEL=(RATLZ-RATLZT)*LZTWM
CV.K  ONLY UNFROZEN WATER CAN BE TRANSFERED
c       if(lzfsc .ne. lzfsh) write(*,*) 'BST4=',lzfsc,lzfsh
      SFH=LZFSH+LZFPH
      IF(DEL .GT. SFH) DEL=SFH
      LZFSH=LZFSH-DEL
      IF(LZFSH .GE. 0.0) THEN
C     TRANSFER FROM LZFSC TO LZTWC.
       LZFSC=LZFSC-DEL
c         if(lzfsc .lt. lzfsh) then
c          write(*,*) ' lzfsc1: ',lzfsc,lzfsh,del
c          stop
c         endif
      ELSE
C     IF TRANSFER EXCEEDS LZFSC THEN REMAINDER COMES FROM LZFPC
       LZFPC=LZFPC+LZFSH
       LZFPH=LZFPH+LZFSH
       xx=LZFSH+DEL
       LZFSC=LZFSC-xx
c         if(lzfsc .lt. lzfsh) then
c          write(*,*) ' lzfsc2: ',lzfsc,lzfsh,del,xx
c          stop
c         endif
       LZFSH=0.0
      ENDIF
      LZTWC=LZTWC+DEL
      LZTWH=LZTWH+DEL

CV.K      LZTWC=LZTWC+DEL
CV.K      LZFSC=LZFSC-DEL
CV.K      IF(LZFSC.GE.0.0) GO TO 230
CV.K      LZFPC=LZFPC+LZFSC
CV.K      LZFSC=0.0

CV.K  CHECK UNFROZEN WATER STORAGE
  230 IF (LZTWC.LT.0.00001) THEN
       LZTWC=0.0
       LZTWH=0.0
      ENDIF
C
C     COMPUTE ET FROM ADIMP AREA.-E5
      E5=E1+(RED+E2)*((ADIMC-E1-UZTWC)/(UZTWM+LZTWM))
C      ADJUST ADIMC,ADDITIONAL IMPERVIOUS AREA STORAGE, FOR EVAPORATION.
      ADIMC=ADIMC-E5
      IF(ADIMC.GE.0.0) GO TO 231
C     E5 CAN NOT EXCEED ADIMC.
      E5=E5+ADIMC
      ADIMC=0.0
  231 E5=E5*ADIMP
C     E5 IS ET FROM THE AREA ADIMP.
C.......................................
C     COMPUTE PERCOLATION AND RUNOFF AMOUNTS.
      TWX=PXV+UZTWC-UZTWM
C     TWX IS THE TIME INTERVAL AVAILABLE MOISTURE IN EXCESS
C     OF UZTW REQUIREMENTS.
      IF(TWX.GE.0.0) GO TO 232
C     ALL MOISTURE HELD IN UZTW--NO EXCESS.
      UZTWC=UZTWC+PXV
CV.K  ADJUST UNFROZEN TENSION WATER
      UZTWH=UZTWH+PXV

      TWX=0.0
      GO TO 233
C      MOISTURE AVAILABLE IN EXCESS OF UZTWC STORAGE.
CV.K  232 UZTWC=UZTWM
  232 UZTWH=UZTWH+(UZTWM-UZTWC)
      UZTWC=UZTWM

  233 ADIMC=ADIMC+PXV-TWX
C
C     COMPUTE IMPERVIOUS AREA RUNOFF.
      ROIMP=PXV*SACPAR(4)
C      ROIMP IS RUNOFF FROM THE MINIMUM IMPERVIOUS AREA.
      SIMPVT=SIMPVT+ROIMP
C
C     INITIALIZE TIME INTERVAL SUMS.
      SBF=0.0
      SSUR=0.0
      SIF=0.0
      SPERC=0.0
      SDRO=0.0
      SPBF=0.0
C
C     DETERMINE COMPUTATIONAL TIME INCREMENTS FOR THE BASIC TIME
C     INTERVAL
CV.K      NINC=1.0+0.2*(UZFWC+TWX)
CV.K  PERCOLATE UNFROZEN WATER ONLY
      NINC=1.0+0.2*(UZFWH+TWX)
C     NINC=NUMBER OF TIME INCREMENTS THAT THE TIME INTERVAL
C     IS DIVIDED INTO FOR FURTHER
C     SOIL-MOISTURE ACCOUNTING.  NO ONE INCREMENT
C     WILL EXCEED 5.0 MILLIMETERS OF UZFWC+PAV
      DINC=(1.0/NINC)*DT
C     DINC=LENGTH OF EACH INCREMENT IN DAYS.
      PINC=TWX/NINC
C     PINC=AMOUNT OF AVAILABLE MOISTURE FOR EACH INCREMENT.
C      COMPUTE FREE WATER DEPLETION FRACTIONS FOR
C     THE TIME INCREMENT BEING USED-BASIC DEPLETIONS
C      ARE FOR ONE DAY
CVK INTRODUCED REDUCTION (RUZICE & RLZICE) DUE FROZEN GROUND
CVK HOWEVER, PRIMARY RUNOFF IS UNCHANGED
CVK      DUZ=1.0-((1.0-UZK)**DINC)
CVK      DLZS=1.0-((1.0-LZSK)**DINC)
CVK  Linear transformation for frozen ground
cc      DUZ=1.0-((1.0-UZK*RUZICE)**DINC)
cc      DLZS=1.0-((1.0-LZSK*RLZICE)**DINC)
CVK  Non-linear (correct) transformation for frozen ground
      IF(IVERS .EQ. 0) THEN
       DUZ =1.0-((1.0-SACPAR(3))**DINC)
       DLZS=1.0-((1.0-SACPAR(12))**DINC)
      ELSE
       DUZ=1.0-((1.0-RUZICE)**DINC)
       DLZS=1.0-((1.0-RLZICE)**DINC)
      ENDIF
      DLZP=1.0-((1.0-SACPAR(13))**DINC)
C
C CVK  ADJUSTMENT TO DEPLETIONS DUE TO FROZEN WATER

C.......................................
C     START INCREMENTAL DO LOOP FOR THE TIME INTERVAL.
      DO 240 I=1,NINC
      ADSUR=0.0
C     COMPUTE DIRECT RUNOFF (FROM ADIMP AREA).
      RATIO=(ADIMC-UZTWC)/LZTWM
      IF (RATIO.LT.0.0) RATIO=0.0
      ADDRO=PINC*(RATIO**2)
C     ADDRO IS THE AMOUNT OF DIRECT RUNOFF FROM THE AREA ADIMP.
C
C     COMPUTE BASEFLOW AND KEEP TRACK OF TIME INTERVAL SUM.
CV.K      BF=LZFPC*DLZP
CV.K      LZFPC=LZFPC-BF
CV.K      IF (LZFPC.GT.0.0001) GO TO 234
CV.K      BF=BF+LZFPC
CV.K      LZFPC=0.0
CV.K  BASEFLOW FROM UNFROZEN WATER ONLY
      BF=LZFPH*DLZP
      LZFPH=LZFPH-BF
      IF (LZFPH.GT.0.0001) THEN
       LZFPC=LZFPC-BF
       GO TO 234
      ENDIF
      BF=BF+LZFPH
      LZFPH=0.0
      LZFPC=LZFPC-BF
      IF(LZFPC .LE. 0.0001) LZFPC=0.0
CV.K-------------------------------------
C
  234 SBF=SBF+BF
      SPBF=SPBF+BF
CV.K  SUPPLAMENTAL FLOW FROM UNFROZEN WATER ONLY (NOTE, DLZS
CV.K  NOTE, DLZS IS REDUCED DUE FROZEN GROUND
CV.K      BF=LZFSC*DLZS
CV.K      LZFSC=LZFSC-BF
CV.K      IF(LZFSC.GT.0.0001) GO TO 235
CV.K      BF=BF+LZFSC
CV.K      LZFSC=0.0
      BF=LZFSH*DLZS
      LZFSH=LZFSH-BF
      IF(LZFSH.GT.0.0001) THEN
cc?      IF(LZFSH.GT.0.0) THEN
       LZFSC=LZFSC-BF
c         if(abs(lzfsc-lzfsh) .gt. 0.000001) then
c         if(abs(lzfsc-lzfsh) .gt. 0.000001) then
c          write(*,*) ' lzfsc3: ',lzfsc,lzfsh,bf
c         endif
       GO TO 235
      ENDIF
      BF=BF+LZFSH
      LZFSH=0.0
      LZFSC=LZFSC-BF
      IF(LZFSC .LE. 0.0001) LZFSC=0.0
CV.K--------------------------------------------
C
  235 SBF=SBF+BF
C
C      COMPUTE PERCOLATION-IF NO WATER AVAILABLE THEN SKIP
ccvk      IF((PINC+UZFWC).GT.0.01) GO TO 251
      xx1=PINC+UZFWH
      IF(xx1.GT.0.01) GO TO 251
      UZFWC=UZFWC+PINC
CV.K  ADD TO UNFROZEN WATER ALSO
      UZFWH=UZFWH+PINC
      GO TO 249
  251 PERCM=LZFPM*DLZP+LZFSM*DLZS
CVK      PERC=PERCM*(UZFWC/UZFWM)
CV.K  USE ONLY UNFROZEN WATER RATIOS
ccvk  new change: PERCOLATION REDUCED BY RUZPERC
CC       PERC=PERCM*(UZFWH/UZFWM)*RUZICE
      PERC=PERCM*(UZFWH/UZFWM)
      IF(IVERS .NE. 0) PERC=PERC*RUZPERC
C--      PERC=PERCM*(UZFWH/UZFWM)*RUZPERC

CV.K      DEFR=1.0-((LZTWC+LZFPC+LZFSC)/(LZTWM+LZFPM+LZFSM))
cvk 6/22/00      DEFR=1.0-((LZTWH+LZFPH+LZFSH)/(LZTWM+LZFPM+LZFSM))
cvk  better to keep original definition of DEFR using total water
      DEFR=1.0-((LZTWC+LZFPC+LZFSC)/(LZTWM+LZFPM+LZFSM))

C     DEFR IS THE LOWER ZONE MOISTURE DEFICIENCY RATIO
c--      FR=1.0
C     FR IS THE CHANGE IN PERCOLATION WITHDRAWAL DUE TO FROZEN GROUND.
c--      FI=1.0
C     FI IS THE CHANGE IN INTERFLOW WITHDRAWAL DUE TO FROZEN GROUND.
c--      IF (IFRZE.EQ.0) GO TO 239
c--       UZDEFR=1.0-((UZTWC+UZFWC)/(UZTWM+UZFWM))
CVK
CVK     CALL FGFR1(DEFR,FR,UZDEFR,FI)
CVK      IF( IVERS .EQ. 1) THEN
CVK  IF IVERS=1, OLD VERSION; IF IVERS=2, NEW VERS. FROST INDEX,
CVK  BUT OLD VERS. OF PERCOLAT. AND INTERFLOW REDUCTION
c--      IF( IVERS .LE. 2) CALL FGFR1(DEFR,FR,UZDEFR,FI)

c--      IF(IVERS .EQ. 3 .AND. FGPM(5) .GT. 0.) THEN
CVK  OPTIONAL VERSION TO ACCOUNT FOR ADDITIONAL IMPERVIOUS
CVK  AREAS EFFECTS DUE FROZEN GROUND
c--       FR=1-SURFRZ1(FGCO(1),FGPM(6),FGPM(5))
c--       FI=FR
c--      ENDIF

c--  239 PERC=PERC*(1.0+ZPERC*(DEFR**REXP))*FR
  239 PERC=PERC*(1.0+SACPAR(7)*(DEFR**SACPAR(8)))
C     NOTE...PERCOLATION OCCURS FROM UZFWC BEFORE PAV IS ADDED.
CV.K      IF(PERC.LT.UZFWC) GO TO 241
      IF(PERC.LT.UZFWH) GO TO 241
C      PERCOLATION RATE EXCEEDS UZFWH.
CV.K      PERC=UZFWC
      PERC=UZFWH
C     PERCOLATION RATE IS LESS THAN UZFWH.
  241 UZFWC=UZFWC-PERC
CV.K  ADJUST UNFROZEN STORAGE ALSO
      UZFWH=UZFWH-PERC

C     CHECK TO SEE IF PERCOLATION EXCEEDS LOWER ZONE DEFICIENCY.
      CHECK=LZTWC+LZFPC+LZFSC+PERC-LZTWM-LZFPM-LZFSM
      IF(CHECK.LE.0.0) GO TO 242
      PERC=PERC-CHECK
      UZFWC=UZFWC+CHECK
CV.K  ADJUST UNFROZEN STARAGE ALSO
      UZFWH=UZFWH+CHECK

  242 SPERC=SPERC+PERC
C     SPERC IS THE TIME INTERVAL SUMMATION OF PERC
C
C     COMPUTE INTERFLOW AND KEEP TRACK OF TIME INTERVAL SUM.
C     NOTE...PINC HAS NOT YET BEEN ADDED
CV.K      DEL=UZFWC*DUZ*FI
CVK  INTERFLOW ALSO REDUCED DUE FROFEN GROUND (DUZ REDUCED BY RUZICE)
CVK  ADDITIONAL REDUCTION DUE IMPERVIOUS FROZEN AREAS (FI) IS OPTIONAL
CVK  IN THE NEW VERSION. BASIC OPTION IS FI=1
c--      DEL=UZFWH*DUZ*FI
      DEL=UZFWH*DUZ
      SIF=SIF+DEL
      UZFWC=UZFWC-DEL
CV.K  ADJUST UNFROZEN STORAGE ALSO
      UZFWH=UZFWH-DEL

C     DISTRIBE PERCOLATED WATER INTO THE LOWER ZONES
C     TENSION WATER MUST BE FILLED FIRST EXCEPT FOR THE PFREE AREA.
C     PERCT IS PERCOLATION TO TENSION WATER AND PERCF IS PERCOLATION
C         GOING TO FREE WATER.
      PERCT=PERC*(1.0-SACPAR(14))
      xx1=PERCT+LZTWC
      IF (xx1.GT.LZTWM) GO TO 243
      LZTWC=LZTWC+PERCT
CV.K  ADJUST UNFROZEN STORAGE ALSO
      LZTWH=LZTWH+PERCT
      PERCF=0.0
      GO TO 244
  243 PERCF=PERCT+LZTWC-LZTWM
CV.K  CHANGE UNFROZEN WATER STORAGE
      LZTWH=LZTWH+LZTWM-LZTWC
      LZTWC=LZTWM
C
C      DISTRIBUTE PERCOLATION IN EXCESS OF TENSION
C      REQUIREMENTS AMONG THE FREE WATER STORAGES.
  244 PERCF=PERCF+PERC*SACPAR(14)
      IF(PERCF.EQ.0.0) GO TO 245
      HPL=LZFPM/(LZFPM+LZFSM)
C     HPL IS THE RELATIVE SIZE OF THE PRIMARY STORAGE
C     AS COMPARED WITH TOTAL LOWER ZONE FREE WATER STORAGE.

c VK changed to account for ZERO MAX storage
      if(LZFPM .ne. 0.) then
       RATLP=LZFPC/LZFPM
      else
       RATLP = 1.
      endif
      if(LZFSM .ne. 0.) then
       RATLS=LZFSC/LZFSM
      else
       RATLS = 1.
      endif

C     RATLP AND RATLS ARE CONTENT TO CAPACITY RATIOS, OR
C     IN OTHER WORDS, THE RELATIVE FULLNESS OF EACH STORAGE
      FRACP=(HPL*2.0*(1.0-RATLP))/((1.0-RATLP)+(1.0-RATLS))
C     FRACP IS THE FRACTION GOING TO PRIMARY.
      IF (FRACP.GT.1.0) FRACP=1.0
      PERCP=PERCF*FRACP
      PERCS=PERCF-PERCP
C     PERCP AND PERCS ARE THE AMOUNT OF THE EXCESS
C     PERCOLATION GOING TO PRIMARY AND SUPPLEMENTAL
C      STORGES,RESPECTIVELY.
      LZFSC=LZFSC+PERCS
CV.K      IF(LZFSC.LE.LZFSM) GO TO 246
      IF(LZFSC.LE.LZFSM) THEN
       LZFSH=LZFSH+PERCS
       GO TO 246
      ENDIF

      PERCS=PERCS-LZFSC+LZFSM
CV.K  ADJUST UNFROZEN STORAGE ALSO
      LZFSH=LZFSH+PERCS

      LZFSC=LZFSM
  246 LZFPC=LZFPC+(PERCF-PERCS)
C     CHECK TO MAKE SURE LZFPC DOES NOT EXCEED LZFPM.
CV.K      IF (LZFPC.LE.LZFPM) GO TO 245
      IF (LZFPC.LE.LZFPM) THEN
       LZFPH=LZFPH+(PERCF-PERCS)
       GO TO 245
      ENDIF

      EXCESS=LZFPC-LZFPM
      LZTWC=LZTWC+EXCESS
CV.K  ADJUST UNFROZEN STORAGES ALSO
      LZTWH=LZTWH+EXCESS
      LZFPH=LZFPH+(PERCF-PERCS)-EXCESS
      LZFPC=LZFPM
C
C     DISTRIBUTE PINC BETWEEN UZFWC AND SURFACE RUNOFF.
  245 IF(PINC.EQ.0.0) GO TO 249
C     CHECK IF PINC EXCEEDS UZFWM
      xx1=PINC+UZFWC
      IF(xx1.GT.UZFWM) GO TO 248
C     NO SURFACE RUNOFF
      UZFWC=UZFWC+PINC
CV.K  ADJUST UNFROZEN STORAGE ALSO
      UZFWH=UZFWH+PINC
      GO TO 249
C
C     COMPUTE SURFACE RUNOFF (SUR) AND KEEP TRACK OF TIME INTERVAL SUM.
  248 SUR=PINC+UZFWC-UZFWM
      UZFWC=UZFWM
CV.K  ADJUST UNFROZEN STORAGE ALSO
      UZFWH=UZFWH+PINC-SUR
      SSUR=SSUR+SUR*PAREA
      ADSUR=SUR*(1.0-ADDRO/PINC)
C     ADSUR IS THE AMOUNT OF SURFACE RUNOFF WHICH COMES
C     FROM THAT PORTION OF ADIMP WHICH IS NOT
C     CURRENTLY GENERATING DIRECT RUNOFF.  ADDRO/PINC
C     IS THE FRACTION OF ADIMP CURRENTLY GENERATING
C     DIRECT RUNOFF.
      SSUR=SSUR+ADSUR*ADIMP
C
C     ADIMP AREA WATER BALANCE -- SDRO IS THE 6 HR SUM OF
C          DIRECT RUNOFF.
  249 ADIMC=ADIMC+PINC-ADDRO-ADSUR
      xx1=UZTWM+LZTWM
      IF (ADIMC.LE.xx1) GO TO 247
      ADDRO=ADDRO+ADIMC-xx1
      ADIMC=xx1
  247 SDRO=SDRO+ADDRO*ADIMP
      IF (ADIMC.LT.0.00001) ADIMC=0.0
  240 CONTINUE

C.......................................
C     END OF INCREMENTAL DO LOOP.
C.......................................

C     COMPUTE SUMS AND ADJUST RUNOFF AMOUNTS BY THE AREA OVER
C     WHICH THEY ARE GENERATED.
      EUSED=E1+E2+E3
C     EUSED IS THE ET FROM PAREA WHICH IS 1.0-ADIMP-PCTIM
      SIF=SIF*PAREA
C
C     SEPARATE CHANNEL COMPONENT OF BASEFLOW
C     FROM THE NON-CHANNEL COMPONENT
      TBF=SBF*PAREA
C     TBF IS TOTAL BASEFLOW
      BFCC=TBF*(1.0/(1.0+SACPAR(15)))
C     BFCC IS BASEFLOW, CHANNEL COMPONENT
      BFP=SPBF*PAREA/(1.0+SACPAR(15))
      BFS=BFCC-BFP
      IF(BFS.LT.0.0)BFS=0.0
      BFNCC=TBF-BFCC
C     BFNCC IS BASEFLOW,NON-CHANNEL COMPONENT
C
C     ADD TO MONTHLY SUMS.
c--      SINTFT=SINTFT+SIF
c--      SGWFP=SGWFP+BFP
c--      SGWFS=SGWFS+BFS
c--      SRECHT=SRECHT+BFNCC
c--      SROST=SROST+SSUR
c--      SRODT=SRODT+SDRO
C
C     COMPUTE TOTAL CHANNEL INFLOW FOR THE TIME INTERVAL.
      TCI=ROIMP+SDRO+SSUR+SIF+BFCC
        GRND = SIF + BFCC   ! interflow is part of ground flow
CC	GRND = BFCC         ! interflow is part of surface flow
	SURF = TCI - GRND
C
C     COMPUTE E4-ET FROM RIPARIAN VEGETATION.
	E4=(EDMND-EUSED)*SACPAR(6)
C
C     SUBTRACT E4 FROM CHANNEL INFLOW
	TCI=TCI-E4
	IF(TCI.GE.0.0) GO TO 250
	E4=E4+TCI
	TCI=0.0
cc  250 SROT=SROT+TCI
250	CONTINUE
	GRND = GRND - E4
	IF (GRND .LT. 0.) THEN
	   SURF = SURF + GRND
	   GRND = 0.
	 IF (SURF .LT. 0.) SURF = 0.
	END IF
C
C     COMPUTE TOTAL EVAPOTRANSPIRATION-TET
      EUSED=EUSED*PAREA
      TET=EUSED+E5+E4
c--      SETT=SETT+TET
c--      SE1=SE1+E1*PAREA
c--      SE3=SE3+E3*PAREA
c--      SE4=SE4+E4
c--      SE5=SE5+E5
C     CHECK THAT ADIMC.GE.UZTWC
      IF (ADIMC.LT.UZTWC) ADIMC=UZTWC
C
c  Return back SAC states
      SACST(1)=UZTWC
      SACST(2)=UZFWC
      SACST(3)=LZTWC
      SACST(4)=LZFSC
      SACST(5)=LZFPC
      SACST(6)=ADIMC

c new change: check negative states
      do i=1,6
       if(sacst(i) .lt. -1.0) then
        write(*,*) ' SAC state#',i,'<-1.',sacst(i)
C        stop
        ERROR = 1
        return
       endif
       if(sacst(i) .lt. 0.0) sacst(i)=0.0
      enddo
      if(uztwh .lt. 0.0) uztwh=0.0
      if(uzfwh .lt. 0.0) uzfwh=0.0
      if(lztwh .lt. 0.0) lztwh=0.0
      if(lzfsh .lt. 0.0) lzfsh=0.0
      if(lzfph .lt. 0.0) lzfph=0.0
      if(sacst(1) .lt. uztwh) uztwh=sacst(1)
      if(sacst(2) .lt. uzfwh) uzfwh=sacst(2)
      if(sacst(3) .lt. lztwh) lztwh=sacst(3)
      if(sacst(4) .lt. lzfsh) lzfsh=sacst(4)
      if(sacst(5) .lt. lzfph) lzfph=sacst(5)
c new change

c      WRITE(*,*) 'FLAND1: end, UZFWC = ', UZFWC
CVK  NEW VERSION OF FROST INDEX  ------------------------------
       IF (IVERS .NE. 0) THEN
        IF(FRZST(6) .LT. 0.) THEN
         FRZST(6)=FRZST(6)+UZTWH
        ELSE
         FRZST(6)=UZTWH
        ENDIF
        IF(FRZST(7) .LT. 0.) THEN
         FRZST(7)=FRZST(7)+UZFWH
        ELSE
         FRZST(7)=UZFWH
        ENDIF
        IF(FRZST(8) .LT. 0.) THEN
         FRZST(8)=FRZST(8)+LZTWH
        ELSE
         FRZST(8)=LZTWH
        ENDIF
        IF(FRZST(9) .LT. 0.) THEN
         FRZST(9)=FRZST(9)+LZFSH
        ELSE
         FRZST(9)=LZFSH
        ENDIF
        IF(FRZST(10) .LT. 0.) THEN
         FRZST(10)=FRZST(10)+LZFPH
        ELSE
         FRZST(10)=LZFPH
        ENDIF

        CALL FROST2_1(PXV,TA,WE,AESC,SH,FRZPAR,SACPAR,FRZST,SACST,
     +        SACST_PRV,SMC,SH2O,DTFRZ,IDTFRZ,NSOIL,NUPL,NSAC,IVERS,
     +        FRZDUP,FRZDBT,FROST, SMAX, SMCDRY, HRAPX, HRAPY)

CVK_02  NEW OPTION TO INTERPOLATE MODEL SOIL LAYER TEMP. INTO DESIRED LAYERS
      NMOD=NSOIL+1
      DSMOD(1)=0.
      DSMOD(1)=-0.5*FRZPAR(10)
      TSTMP(1)=FRZST(1)
      SWTMP(1)=SMC(2)
      SWHTMP(1)=SH2O(2)
      DSMOD(NMOD)=-FRZPAR(5)
      TSTMP(NMOD)=FRZPAR(2)-T0
      SWTMP(NMOD)=SMAX
      SWHTMP(NMOD)=SMAX
      do i=2,nmod-1
       DSMOD(I)=-0.5*(FRZPAR(I+8)+FRZPAR(I+9))
       TSTMP(I)=FRZST(I)
       SWTMP(I)=SMC(I)
       SWHTMP(I)=SH2O(I)
      ENDDO
cc-      do ii = 1, nmod
cc-        WRITE(*,*) SWTMP(ii), SWHTMP(ii), DSMOD(ii)
cc-      ENDDO

      CALL SOIL_INT1(TSTMP,NMOD,DSMOD,DSINT,NDSINT,TSINT)
      CALL SOIL_INT1(SWTMP,NMOD,DSMOD,DSINTW,NDINTW,SWINT)
      CALL SOIL_INT1(SWHTMP,NMOD,DSMOD,DSINTW,NDINTW,SWHINT)

cvk  1/2008 Option to generate normalized soil moisture content (SR)
      if(NORMALIZE .eq. 1) then
       DO I=1,NDINTW
        SWINT(I) =  (SWINT(I) -  FRZPAR(9))/(SMAX-FRZPAR(9))
        SWHINT(I) = (SWHINT(I) - FRZPAR(9))/(SMAX-FRZPAR(9))
       ENDDO
      endif
cvk  1/2008  end soil moisture normalization

C--      DO I=1,NINTW
C--       IF(I .EQ. 1) THEN
C--        SWINT(I)=SWINT(I)*DSINTW(I)*1000.
C--        SWHINT(I)=SWHINT(I)*DSINTW(I)*1000
C--       ELSE
C--        SWINT(I)=SWINT(I)*(DSINTW(I)-DSINTW(I-1))*1000.
C--        SWHINT(I)=SWHINT(I)*(DSINTW(I)-DSINTW(I-1))*1000
C--       ENDIF
C--      ENDDO

c        WRITE (*,905) (sacst(ii),ii=1,6),
c     1  SPERC,ROIMP,SDRO,SSUR,SIF,BFS,BFP,TCI,
c     2  EDMND,TET,PXV,(frzst(II),II=6,10),WE,SH,AESC,TA,frost,
c     3  (frzst(II),II=1,nsoil),frzdup,frzdbt
c        write(*,977) smax,(swint(ii),swhint(ii),tstmp(ii),ii=1,nintw)
       ELSE
cc        WRITE (*,977) (sacst(ii),ii=1,6),
cc     1  SPERC,ROIMP,SDRO,SSUR,SIF,BFS,BFP,TCI,
cc     2  EDMND,TET,PXV,UZTWH,UZFWH,LZTWH,LZFSH,LZFPH,WE,SH,AESC,TA
       ENDIF
  905 FORMAT (1H ,14x,F7.2,F7.3,F7.2,F7.3,2F7.2,7F7.3,2F8.3,F7.3,
     +        F9.4,5f7.2,2f6.1,F6.3,F9.4,f7.2,8f7.2,20F7.1)
  977 FORMAT (1H ,14x,f7.2,6(3F7.2))
c             ,F7.3,F7.2,F7.3,2F7.2,7F7.3,2F8.3,F7.3,
c     +        F9.4,5F7.2,2f6.1,F6.3,F9.4,8f7.2,20F7.1)

C.......................................

      RETURN
      END

max_stars_repo_name: KathyNie/LISF

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF


text:
program block_test
    implicit none
    integer, parameter :: n = 3
    integer :: i
    integer, dimension(n, n) :: a = reshape( [ (i, i = 1, n*n) ], [n, n] )

    print '(A)', 'block'
    do i = 1, n
        block
            integer, dimension(n) :: row
            row = a(:, i)
            print '(A25, 5I4)', 'row: ', row
            row(1) = 123
            print '(A25, 3I4)', 'row after update: ', row
            print '(A25, 3I4)', 'matrix after update: ', a(:, i)
        end block
    end do
    print '(A)', 'associate'
    do i = 1, n
        associate( row => a(:, i) )
            print '(A25, 5I4)', 'row: ', row
            row(1) = 123
            print '(A25, 3I4)', 'row after update: ', row
            print '(A25, 3I4)', 'matrix after update: ', a(:, i)
        end associate
    end do
end program block_test

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


module comoving_nd_module

  use amrex_error_module
  use amrex_fort_module, only : rt => amrex_real

  contains

! :::
! ::: ----------------------------------------------------------------
! :::

      real(rt) function invEz(H0, Om, a)
        real(rt), intent(in) :: H0, Om, a
        invEz = 1.0d0 / ( H0*dsqrt(Om/a + (1.0d0-Om)*a*a) )
      end function invEz

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_integrate_time_given_a(a0, a1, dt) &
         bind(C, name="fort_integrate_time_given_a")

        use fundamental_constants_module, only: Hubble_const
        use comoving_module             , only: comoving_h, comoving_OmM

        real(rt), intent(in   ) :: a0, a1
        real(rt), intent(  out) :: dt

        real(rt), parameter :: xacc = 1.0d-6
        real(rt) :: H0, Om, prev_soln, h
        integer :: iter, n, j

        H0 = comoving_h*Hubble_const
        Om = comoving_OmM

        prev_soln = -1.0d0
        ! trapezoidal integration
        do iter = 1, 20  ! max allowed iterations

          n  = 2**iter

          h = (a1-a0)/(n-1)

          if (a0 .lt. 1.0d-10) then  ! prevent division by zero in invEz
             dt = 0.5*invEz(H0, Om, a1)
          else
             dt = 0.5*(invEz(H0, Om, a0) + invEz(H0, Om, a1))
          endif

          do j = 1, n-2
             dt = dt + invEz(H0, Om, a0+j*h)
          enddo
          dt = dt*h

          if (iter .gt. 4) then
             if (abs(dt-prev_soln) .le. xacc*abs(prev_soln)) return
          endif
          prev_soln = dt
        enddo

      end subroutine fort_integrate_time_given_a

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_integrate_comoving_a(old_a,new_a,dt) &
         bind(C, name="fort_integrate_comoving_a")

        use fundamental_constants_module, only: Hubble_const
        use comoving_module             , only: comoving_h, comoving_OmM, comoving_type

        implicit none

        real(rt), intent(in   ) :: old_a, dt
        real(rt), intent(  out) :: new_a

        real(rt), parameter :: xacc = 1.0d-8
        real(rt) :: H_0, OmL
        real(rt) :: Delta_t, prev_soln
        real(rt) :: start_a, end_a, start_slope, end_slope
        integer  :: iter, j, nsteps

        if (comoving_h .eq. 0.0d0) then
          new_a = old_a
          return
        endif

        H_0 = comoving_h * Hubble_const
        OmL = 1.d0 - comoving_OmM

        prev_soln = 2.0d0 ! 0<a<1 so a=2 will do as "wrong" solution
        do iter = 1, 11  ! max allowed iterations
          nsteps = 2**(iter-1)
          Delta_t = dt/nsteps
          end_a = old_a

          do j = 1, nsteps
            ! This uses RK2 to integrate the ODE:
            !   da / dt = H_0 * sqrt(OmM/a + OmL*a^2)
            start_a = end_a

            ! Compute the slope at the old time
            if (comoving_type > 0) then
                start_slope = H_0*dsqrt(comoving_OmM / start_a + OmL*start_a**2)
            else
                start_slope = comoving_h
            end if

            ! Compute a provisional value of ln(a) at the new time
            end_a = start_a + start_slope * Delta_t

            ! Compute the slope at the new time
            if (comoving_type > 0) then
                end_slope = H_0*dsqrt(comoving_OmM / end_a + OmL*end_a**2)
            else
                end_slope = comoving_h
            end if

            ! Now recompute a at the new time using the average of the two slopes
            end_a = start_a + 0.5d0 * (start_slope + end_slope) * Delta_t
          enddo

          new_a  = end_a
          if (abs(1.0d0-new_a/prev_soln) .le. xacc) return
          prev_soln = new_a

        enddo

      end subroutine fort_integrate_comoving_a

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_integrate_comoving_a_to_z(old_a,z_value,dt) &
         bind(C, name="fort_integrate_comoving_a_to_z")

        use fundamental_constants_module, only: Hubble_const
        use comoving_module             , only: comoving_h, comoving_OmM, comoving_type

        implicit none

        real(rt), intent(in   ) :: old_a, z_value
        real(rt), intent(inout) :: dt

        real(rt), parameter :: xacc = 1.0d-8
        real(rt) :: H_0, OmL
        real(rt) :: Delta_t
        real(rt) :: start_a, end_a, start_slope, end_slope
        real(rt) :: a_value
        integer  :: j, nsteps

        if (comoving_h .eq. 0.0d0) &
          call amrex_error("fort_integrate_comoving_a_to_z: Shouldn't be setting plot_z_values if not evolving a")

        H_0 = comoving_h * Hubble_const
        OmL = 1.d0 - comoving_OmM

        ! Translate the target "z" into a target "a"
        a_value = 1.d0 / (1.d0 + z_value)

        ! Use lots of steps if we want to nail the z_value
        nsteps = 1024

        ! We integrate a, but stop when a = a_value (or close enough)
        Delta_t = dt/nsteps
        end_a = old_a
        do j = 1, nsteps
            ! This uses RK2 to integrate the ODE:
            !   da / dt = H_0 * sqrt(OmM/a + OmL*a^2)
            start_a = end_a

            ! Compute the slope at the old time
            if (comoving_type > 0) then
                start_slope = H_0*dsqrt(comoving_OmM / start_a + OmL*start_a**2)
            else
                start_slope = comoving_h
            end if

            ! Compute a provisional value of ln(a) at the new time
            end_a = start_a + start_slope * Delta_t

            ! Compute the slope at the new time
            if (comoving_type > 0) then
                end_slope = H_0*dsqrt(comoving_OmM / end_a + OmL*end_a**2)
            else
                end_slope = comoving_h
            end if

            ! Now recompute a at the new time using the average of the two slopes
            end_a = start_a + 0.5d0 * (start_slope + end_slope) * Delta_t

            ! We have crossed from a too small to a too big in this step
            if ( (end_a - a_value) * (start_a - a_value) < 0) then
                dt = ( (  end_a - a_value) * dble(j  ) + &
                       (a_value - start_a) * dble(j+1) ) / (end_a - start_a) * Delta_t
                exit
            end if
        end do

      end subroutine fort_integrate_comoving_a_to_z

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_integrate_comoving_a_to_a(old_a,a_value,dt) &
         bind(C, name="fort_integrate_comoving_a_to_a")

        use fundamental_constants_module, only: Hubble_const
        use comoving_module             , only: comoving_h, comoving_OmM, comoving_type

        implicit none

        real(rt), intent(in   ) :: old_a, a_value
        real(rt), intent(inout) :: dt

        real(rt), parameter :: xacc = 1.0d-8
        real(rt) :: H_0, OmL
        real(rt) :: Delta_t
        real(rt) :: start_a, end_a, start_slope, end_slope
        integer  :: j, nsteps

        real(rt) :: max_dt

        if (comoving_h .eq. 0.0d0) &
          call amrex_error("fort_integrate_comoving_a_to_z: Shouldn't be setting plot_z_values if not evolving a")

        H_0 = comoving_h * Hubble_const
        OmL = 1.d0 - comoving_OmM

        ! Use lots of steps if we want to nail the z_value
!        nsteps = 1024

        ! Use enough steps if we want to be close to the a_value
        nsteps = 2048

        ! We integrate a, but stop when a = a_value (or close enough)
        Delta_t = dt/nsteps
        end_a = old_a
        do j = 1, nsteps
            ! This uses RK2 to integrate the ODE:
            !   da / dt = H_0 * sqrt(OmM/a + OmL*a^2)
            start_a = end_a

            ! Compute the slope at the old time
            if (comoving_type > 0) then
                start_slope = H_0*dsqrt(comoving_OmM / start_a + OmL*start_a**2)
            else
                start_slope = comoving_h
            end if

            ! Compute a provisional value of ln(a) at the new time
            end_a = start_a + start_slope * Delta_t

            ! Compute the slope at the new time
            if (comoving_type > 0) then
                end_slope = H_0*dsqrt(comoving_OmM / end_a + OmL*end_a**2)
            else
                end_slope = comoving_h
            end if

            ! Now recompute a at the new time using the average of the two slopes
            end_a = start_a + 0.5d0 * (start_slope + end_slope) * Delta_t

            ! We have crossed from a too small to a too big in this step
            if ( (end_a - a_value) * (start_a - a_value) < 0) then
                dt = ( (  end_a - a_value) * dble(j  ) + &
                       (a_value - start_a) * dble(j+1) ) / (end_a - start_a) * Delta_t
                exit
            end if
        end do

      end subroutine fort_integrate_comoving_a_to_a
! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_est_maxdt_comoving_a(old_a,dt)

        use fundamental_constants_module, only: Hubble_const
        use comoving_module             , only: comoving_h, comoving_OmM, comoving_type

        implicit none

        real(rt), intent(in   ) :: old_a
        real(rt), intent(inout) :: dt

        real(rt) :: H_0, OmL
        real(rt) :: max_dt

        OmL = 1.d0 - comoving_OmM

        ! This subroutine computes dt based on not changing a by more than 5%
        ! if we use forward Euler integration
        !   d(ln(a)) / dt = H_0 * sqrt(OmM/a^3 + OmL)

        H_0 = comoving_h * Hubble_const

        if (H_0 .ne. 0.0d0) then
           if (comoving_type > 0) then
              max_dt = (0.05d0) / H_0 / dsqrt(comoving_OmM / old_a**3 + OmL)
           else
              max_dt = (0.05d0) / abs(comoving_h)
           end if
           dt = min(dt,max_dt)

        else

           ! dt is unchanged

        end if

      end subroutine fort_est_maxdt_comoving_a

! :::
! ::: ----------------------------------------------------------------
! :::

      ! This might only work for t=0=> a=.00625, although constant canceled
      subroutine fort_est_lindt_comoving_a(old_a,new_a,dt)

        use fundamental_constants_module, only: Hubble_const
        use comoving_module             , only: comoving_h, comoving_OmM, comoving_type

        implicit none

        real(rt), intent(in   ) :: old_a,new_a
        real(rt), intent(inout) :: dt

        real(rt) :: H_0, OmL
        real(rt) :: lin_dt

        OmL = 1.d0 - comoving_OmM

        ! This subroutine computes dt based on not changing a by more than 5%
        ! if we use forward Euler integration
        !   d(ln(a)) / dt = H_0 * sqrt(OmM/a^3 + OmL)

        H_0 = comoving_h * Hubble_const

        ! Could definately be optimized better
        if (H_0 .ne. 0.0d0) then

           lin_dt= ((new_a/(.75**(2/3)*(OmL+ comoving_OmM)**(1/3)))**(.75)  - &
                ((old_a/(.75**(2/3)*(OmL+ comoving_OmM)**(1/3)))**(.75) ) ) /H_0
           dt=lin_dt

        else

           ! dt is unchanged

        end if

      end subroutine fort_est_lindt_comoving_a

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_estdt_comoving_a(old_a,new_a,dt,change_allowed,fixed_da,final_a,dt_modified) &
         bind(C, name="fort_estdt_comoving_a")

        use comoving_module             , only: comoving_h

        implicit none

        real(rt), intent(in   ) :: old_a, change_allowed, fixed_da, final_a
        real(rt), intent(inout) :: dt
        real(rt), intent(  out) :: new_a
        integer , intent(  out) :: dt_modified
        real(rt) a_value
        real(rt) max_dt
        max_dt = dt

        if (comoving_h .ne. 0.0d0) then

           if( fixed_da .le. 0.0d0) then
              ! First call this to make sure dt that we send to integration routine isnt outrageous
              call fort_est_maxdt_comoving_a(old_a,dt)

              ! Initial call to see if existing dt will work
              call fort_integrate_comoving_a(old_a,new_a,dt)

              ! Make sure a isn't growing too fast
              call enforce_percent_change(old_a,new_a,dt,change_allowed)
           else
              ! First call this to make sure dt that we send to integration routine isnt outrageous
              new_a = (old_a +  fixed_da);
              call fort_est_lindt_comoving_a(old_a,new_a,dt)
              call fort_est_maxdt_comoving_a(old_a,dt)

              ! Then integrate old_a to a_value using dt as a guess for the maximum dt
              ! Output dt is based on a fraction of the input dt
              call fort_integrate_comoving_a_to_a(old_a,new_a,dt)
           endif

           ! Make sure we don't go past final_a (if final_a is set)
           if (final_a > 0.0d0) &
              call enforce_final_a(old_a,new_a,dt,final_a)

           dt_modified = 1

        else

           ! dt is unchanged by this call

           dt_modified = 0

        endif

      end subroutine fort_estdt_comoving_a

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine enforce_percent_change(old_a,new_a,dt,change_allowed)

        implicit none

        real(rt), intent(in   ) :: old_a, change_allowed
        real(rt), intent(inout) :: dt
        real(rt), intent(inout) :: new_a

        integer          :: i
        real(rt) :: factor

        factor = ( (new_a - old_a) / old_a ) / change_allowed

        ! Only go into this process if percent change exceeds change_allowed

        if (factor > 1.d0) then

           do i = 1, 100
              factor = ( (new_a - old_a) / old_a ) / change_allowed

              ! Note: 0.99 is just a fudge factor so we don't get bogged down.
              if (factor > 1.d0) then
                 dt = (1.d0 / factor) * dt * 0.99d0
                 call fort_integrate_comoving_a(old_a,new_a,dt)
              else if (i.lt.100) then
                 call fort_integrate_comoving_a(old_a,new_a,dt)
                 ! We're done
                 return
              else
                 call amrex_error("Too many iterations in enforce_percent_change")
              end if
           end do

        else
           ! We don't need to do anything
           return
        end if

      end subroutine enforce_percent_change

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine enforce_final_a(old_a,new_a,dt,final_a)

        implicit none

        real(rt), intent(in   ) :: old_a, final_a
        real(rt), intent(inout) :: dt
        real(rt), intent(inout) :: new_a

        integer  :: i
        real(rt) :: factor
        real(rt), parameter :: eps = 1.d-10

        if (old_a > final_a) then
           call amrex_error("Oops -- old_a > final_a")
        end if

        ! Only go into this process if new_a is past final_a
        if (new_a > final_a) then

           do i = 1, 100
              if ( (new_a > (final_a+eps)) .or. (new_a < final_a) ) then
                 factor = (final_a - old_a) / (new_a - old_a)
                 dt = dt * factor
                 call fort_integrate_comoving_a(old_a,new_a,dt)
              else if (i.lt.100) then
                 ! We're done
                 return
              else
                 call amrex_error("Too many iterations in enforce_final_a")
              end if
           end do

        else
           ! We don't need to do anything
           return
        end if

      end subroutine enforce_final_a

! ! :::
! ! ::: ----------------------------------------------------------------
! ! :::

!       subroutine fort_get_omb(frac) &
!          bind(C, name="fort_get_omb")

!         use comoving_module, only: comoving_OmB, comoving_OmM

!         real(rt) :: frac

!         frac = comoving_OmB / comoving_OmM

!       end subroutine fort_get_omb

! ! :::
! ! ::: ----------------------------------------------------------------
! ! :::

!       subroutine fort_get_omm(omm) &
!          bind(C, name="fort_get_omm")

!         use comoving_module, only: comoving_OmM

!         real(rt) :: omm

!         omm = comoving_OmM

!       end subroutine fort_get_omm

! ! :::
! ! ::: ----------------------------------------------------------------
! ! :::

        subroutine fort_get_hubble(hubble) &
           bind(C, name="fort_get_hubble")

          use comoving_module, only: comoving_h

          real(rt) :: hubble

          hubble = comoving_h

        end subroutine fort_get_hubble

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_set_omb(omb) &
         bind(C, name="fort_set_omb")

        use comoving_module, only: comoving_OmB

        real(rt), intent(in) :: omb

        comoving_OmB = omb

      end subroutine fort_set_omb

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_set_omm(omm) &
         bind(C, name="fort_set_omm")

        use comoving_module, only: comoving_OmM

        real(rt), intent(in) :: omm

        comoving_OmM = omm

      end subroutine fort_set_omm

! :::
! ::: ----------------------------------------------------------------
! :::

      subroutine fort_set_hubble(hubble) &
         bind(C, name="fort_set_hubble")

        use comoving_module, only: comoving_h

        real(rt), intent(in) :: hubble

        comoving_h = hubble

      end subroutine fort_set_hubble

end module comoving_nd_module


text:
*
*  -- LEN_TRIM is Fortran 95, so we use a replacement here
*
      FUNCTION LEN_TRIM(S)
*
      CHARACTER*(*)  S
      INTEGER LEN_TRIM
*
      INTRINSIC LEN
*
      DO LEN_TRIM = LEN(S), 1, -1
         IF (s(LEN_TRIM:LEN_TRIM) .NE. ' ') RETURN
      END DO
      END

max_stars_repo_name: globusgenomics/galaxy

max_issues_repo_name: globusgenomics/galaxy

max_forks_repo_name: globusgenomics/galaxy


text:
      module LOADBALANCER_COMMUNICATION
      use iso_fortran_env
      type t_lb_comm_entry
        integer(int32)                           :: proc
        integer(int32)                           :: length
        integer(int32)                           :: data_buf_start
        integer(int32)                           :: data_buf_end
!---- These entries are for communication the other way around
        integer(int32)                           :: result_buf_start
        integer(int32)                           :: result_buf_end
      end type

      type t_lb_comm
        integer(int32)                           :: exports_length
        integer(int32)                           :: imports_length
        type(t_lb_comm_entry), allocatable       :: exports(:)
        type(t_lb_comm_entry), allocatable       :: imports(:)
      end type

      type BYTEsliceptr
        BYTE, pointer,contiguous, dimension(:) :: p
      end type

      end module LOADBALANCER_COMMUNICATION

max_stars_repo_name: NLESC-JCER/QUICKLB

max_issues_repo_name: NLESC-JCER/QUICKLB

max_forks_repo_name: NLESC-JCER/QUICKLB


text:
!======================================================================!
  REAL FUNCTION Rnd (r)
!----------------------------------------------------------------------!
!   Computes random double precision number in the range 0.0-1.0.      !
!----------------------------------------------------------------------!
!-----------------------------[Parameters]-----------------------------!
  INTEGER :: r
!--------------------------------[CVS]---------------------------------!
!  $Id: Rnd.f90,v 1.4 2002/10/30 16:29:33 niceno Exp $
!  $Source: /home/muhamed/.CVSROOT/T-Rex/Library/Rnd.f90,v $
!======================================================================!

  call init_random_seed(r)
  call random_number(rnd)

  END FUNCTION Rnd

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text:
  program testffjet

  use class_four_vector
  use fluid_model_ffjet, only: initialize_ffjet_model, ffjet_vals, del_ffjet_data
  use kerr, only: calc_rms

  type (four_vector), dimension(:), allocatable :: x0,b,u
  real, dimension(:,:), allocatable :: r2d,z2d
  real, dimension(:), allocatable :: rho,p,r,z,rtemp,ztemp,bmag
  real :: a,rmin,rmax,zmin,zmax
  integer :: nr, i, k
  integer, dimension(:), allocatable :: indx

  nr=128; nz=128
  a=0.998; rmin=calc_rms(a); rmax=25.
  zmin=rmin; zmax=100.

  allocate(rtemp(nr)); allocate(ztemp(nz))
  allocate(z(nr*nz)); allocate(r(nr*nz)); allocate(rho(nr*nz)); allocate(p(nr*nz))
  allocate(r2d(nr,nz)); allocate(z2d(nr,nz)); allocate(bmag(nr*nz))

  allocate(x0(nr*nz)); allocate(b(nr*nz)); allocate(u(nr*nz))

  ztemp=(/(zmin*exp(real(k-1)/(nz-1)*log(zmax/zmin)),k=1,nz)/)
  rtemp=(/(rmin*exp(real(k-1)/(nr-1)*log(rmax/rmin)),k=1,nr)/)

  do i=1,nz
     r2d(:,i)=rtemp
  enddo
  do i=1,nr
     z2d(i,:)=ztemp
  enddo

  r=reshape(r2d,(/nr*nz/))
  z=reshape(z2d,(/nr*nz/))

!  write(6,*) 'ztemp: ',zmin,zmax,ztemp
!  write(6,*) 'rtemp: ',rmin,rmax,rtemp
!  write(6,*) 'after: ',size(r),size(z),size(r2d),size(z2d)

  call initialize_ffjet_model(dble(a))

  write(6,*) 'after init'

! Set up x four-vector:
  x0%data(1)=0.; x0%data(2)=sqrt(r*r+z*z)
  x0%data(3)=acos(z/sqrt(r*r+z*z)); x0%data(4)=0.

  call ffjet_vals(x0,a,rho,p,b,u,bmag)
! Make proper arrays of rho, v, b to check with BL09

!  write(6,*) 'after vals rho', rho
!  write(6,*) 'after vals b', b*b

! Output so I can plot in IDL

  open(unit=12,file='test_ffjet.out')
  write(12,*) nr*nz, a
  write(12,*) r, z
  write(12,*) rho, p
  write(12,*) b*b
  close(unit=12)

  write(6,*) 'unit test ffjet: ',maxval(abs(u*u+1.)),maxval(u*b),minval(rho)
  open(unit=12,file='unit_test_ffjet.out')
  write(12,*) minval(rho), maxval(abs(u*u+1.)), maxval(abs(u*b))
  close(unit=12)

  call del_ffjet_data()

  deallocate(rtemp); deallocate(ztemp); deallocate(z); deallocate(r); deallocate(rho)
  deallocate(p); deallocate(b); deallocate(x0); deallocate(u); deallocate(bmag)

  end program

max_stars_repo_name: HerculesJack/grtrans

max_issues_repo_name: RAnantua/grtrans

max_forks_repo_name: RAnantua/grtrans


text:
! { dg-do run { target fd_truncate } }
! Wide character I/O test 3, unformatted arrays
! Test case developed by Jerry DeLisle <jvdelisle@gcc.gnu.org>
program test1
  integer, parameter :: k4 = 4
  character(len=10,kind=4) :: wide
  character(len=10,kind=4), dimension(5,7) :: widearray
  wide = k4_"abcdefg"
  widearray = k4_"1234abcd"
  open(10, form="unformatted", status="scratch")
  write(10) wide
  rewind(10)
  wide = "wrong"
  read(10) wide
  if (wide /= k4_"abcdefg") call abort
  rewind(10)
  write(10) widearray(2:4,3:7)
  widearray(2:4,3:7)=""
  rewind(10)
  read(10) widearray(2:4,3:7)
  close(10)
  if (any(widearray.ne.k4_"1234abcd")) call abort
end program test1

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab