text:
module nonradiative
	use constants
	use ioutil
	implicit none

	type sysdata
		logical											:: do_rad, do_nonrad, do_write, skip_rate
		integer											:: ncut, nthresh, nlevels, nsamples, natoms, debug_level
		integer											:: maxnfix, minnfix, nrecords, lambda_n
		integer(bigint)									:: maxnoccs
		character(len=100)								:: bfile, gradfile, hessfile, algorithm, weighting
		character(len=100)								:: radfile, calctype, radunits, sortby, occprefix
		real(dbl)										:: k_ic, k_ic_ht, k_r, e_target, delta_e, gamma, tdm, memory, user_gamma
		real(dbl), dimension(:), allocatable			:: energies, hrfactors, masses, V_vq_j
		real(dbl), dimension(:, :), allocatable			:: fcfactors, V_jj, Bvqj
		integer, dimension(:, :), allocatable			:: cutoffs, bounds
		integer(sint), dimension(:), allocatable		:: energy_order
		integer(smallint), dimension(:), allocatable	:: minoccs, maxoccs
		type(memorymanager)								:: mm
	contains
		procedure	:: from_file => sysdata_from_file
		procedure	:: init => sysdata_init
		procedure	:: fc_compute => sysdata_fc_compute
		procedure	:: cutoff_compute => sysdata_cutoff_compute
		procedure	:: compute_zn => sysdata_compute_zn
		procedure	:: build_V => sysdata_build_V
		procedure	:: calculate_kic => sysdata_calc_kic
		procedure	:: calculate_from_lagrange => sysdata_calc_lagrange
		procedure	:: calculate_from_file => sysdata_calc_from_file
		procedure	:: calculate_gamma => sysdata_calc_gamma
		procedure	:: get_next_bound => sysdata_get_next_bound
		procedure	:: screened_ncombinations => sysdata_ncombinations
		procedure	:: free => sysdata_free
	end type sysdata

contains

	subroutine sysdata_from_file(sys, filename)
		class(sysdata), intent(inout)	:: sys
		character(len=*), intent(in)	:: filename

	    ! Input related variables
	    character(len=100) :: buffer, label, randostring
	    integer 	:: pos
	    integer 	:: ios = 0
	    integer 	:: line = 0
		integer 	:: ix
		real(dbl)	:: tmp_nsamples = 0d0

	    open(main_input_unit, file=filename)

	    ! ios is negative if an end of record condition is encountered or if
	    ! an endfile condition was detected.  It is positive if an error was
	    ! detected.  ios is zero otherwise.

		sys%debug_level = 0
		sys%do_write = .false.
		sys%skip_rate = .false.
		sys%memory = 0.5d0
		sys%nrecords = 0
		sys%maxnfix = -1
		sys%minnfix = -1
		sys%occprefix = 'occs'
		sys%hessfile  = 'none'
		sys%user_gamma=-1.0
	    do while (ios == 0)
	       read(main_input_unit, '(A)', iostat=ios) buffer
	       if (ios == 0) then
	          line = line + 1

	          ! Find the first instance of whitespace.  Split label and data.
	          pos = scan(buffer, ',')
	          label = trim(buffer(1:pos-1))
	          buffer = buffer(pos+1:)

	          select case (label)
			  case ('spectrum')
		  		 read(buffer, *, iostat=ios) sys%radfile
				 sys%radfile = trim(adjustl(sys%radfile))
			  case ('bfile')
			  	 read(buffer, *, iostat=ios) sys%bfile
				 sys%bfile = trim(adjustl(sys%bfile))
   			  case ('gradfile')
   			  	 read(buffer, *, iostat=ios) sys%gradfile
   				 sys%gradfile = trim(adjustl(sys%gradfile))
      		  case ('hessfile')
      		  	 read(buffer, *, iostat=ios) sys%hessfile
      			 sys%hessfile = trim(adjustl(sys%hessfile))
      		  case ('algorithm')
      			 read(buffer, *, iostat=ios) sys%algorithm
      			 sys%algorithm = trim(adjustl(sys%algorithm))
         	  case ('weighting')
         		 read(buffer, *, iostat=ios) sys%weighting
         		 sys%weighting = trim(adjustl(sys%weighting))
              case ('occprefix')
            	 read(buffer, *, iostat=ios) sys%occprefix
            	 sys%occprefix = trim(adjustl(sys%occprefix))
			  case ('calculation')
			  	 read(buffer, *, iostat=ios) sys%calctype
				 sys%calctype = trim(adjustl(sys%calctype))
   			  case ('radunits')
   			  	 read(buffer, *, iostat=ios) sys%radunits
   				 sys%radunits = trim(adjustl(sys%radunits))
      		  case ('sortby')
      			 read(buffer, *, iostat=ios) sys%sortby
      		     sys%sortby = trim(adjustl(sys%sortby))
   	          case ('debug')
   	             read(buffer, *, iostat=ios) sys%debug_level
			  case ('nrecords')
			  	 read(buffer, *, iostat=ios) sys%nrecords
			  case ('dump')
			  	 read(buffer, *, iostat=ios) sys%mm%dumplevel
	          case ('ncut')
	             read(buffer, *, iostat=ios) sys%ncut
			  case ('tofile')
			  	 read(buffer, *, iostat=ios) randostring
				 if (trim(adjustl(randostring)) .eq. 'true') then
					 sys%do_write = .true.
				 end if
   			  case ('norate')
   			  	 read(buffer, *, iostat=ios) randostring
   				 if (trim(adjustl(randostring)) .eq. 'true') then
   					 sys%skip_rate = .true.
   				 end if
   	          case ('natoms')
   	             read(buffer, *, iostat=ios) sys%natoms
	          case ('nthresh')
	             read(buffer, *, iostat=ios) sys%nthresh
			  case ('nsamples')
			  	 read(buffer, *, iostat=ios) tmp_nsamples
				 sys%nsamples = int(tmp_nsamples)
			  case ('nmax')
			  	 read(buffer, *, iostat=ios) sys%maxnfix
   			  case ('nmin')
   			  	 read(buffer, *, iostat=ios) sys%minnfix
      		  case ('memory')
      			 read(buffer, *, iostat=ios) sys%memory
   			  case ('energy')
   			  	 read(buffer, *, iostat=ios) sys%e_target
   			  case ('deltae')
   			  	 read(buffer, *, iostat=ios) sys%delta_e
			  case ('tdm')
			  	 read(buffer, *, iostat=ios) sys%tdm
      		  case ('gamma')
      			 read(buffer, *, iostat=ios) sys%user_gamma
			  case ('levels')
			  	 read(buffer, *, iostat=ios) sys%nlevels
				 allocate(sys%energies(sys%nlevels))
				 allocate(sys%hrfactors(sys%nlevels))
				 do ix=1,sys%nlevels
					 read(main_input_unit, '(A)', iostat=ios) buffer
					 pos = scan(buffer, ',')
					 label = buffer(1:pos-1)
					 buffer = buffer(pos+1:)
					 read(label, *, iostat=ios) sys%energies(ix)
					 read(buffer, *, iostat=ios) sys%hrfactors(ix)
				 end do
	          case default
	             print *, 'Skipping invalid label at line', line
	          end select
	       end if
	    end do
		close(main_input_unit)
	end subroutine sysdata_from_file

	subroutine sysdata_init(sys)
		class(sysdata), intent(inout)	:: sys

		integer	:: ix, tmp
		real(dbl)	:: nestimate
		real(dbl), dimension(sys%nlevels) :: unsorted_levels

		allocate(sys%energy_order(sys%nlevels))
		do ix=1,sys%nlevels
			sys%energy_order(ix) = ix
		end do

		if (sys%sortby .eq. 'energy') then
			unsorted_levels(:) = sys%energies(:)
			call sort_list(sys%nlevels, unsorted_levels, sys%energies, sys%energy_order)
			call reorder_list(sys%nlevels, sys%hrfactors, sys%energy_order)
		else
			! reverse sort by HR factor by default
			unsorted_levels(:) = -sys%hrfactors(:)
			call sort_list(sys%nlevels, unsorted_levels, sys%hrfactors, sys%energy_order)
			sys%hrfactors = -sys%hrfactors
			call reorder_list(sys%nlevels, sys%energies, sys%energy_order)
		end if

		! calculate maximum no. of occs and initialise memory manager
		if (sys%do_nonrad) then
			call sys%fc_compute
			call sys%cutoff_compute

			allocate(sys%minoccs(sys%nlevels))
			allocate(sys%maxoccs(sys%nlevels))
			sys%minoccs = 0
			sys%maxoccs = [(sys%cutoffs(ix, sys%nthresh+1), ix=1,sys%nlevels)]
			sys%maxnoccs = ncombinations(sys%nlevels, sys%maxoccs, sys%minoccs)

			if (sys%maxnfix .lt. 0) then
				sys%maxnfix = sum(sys%maxoccs)
				if (sys%algorithm .eq. 'fixedn') then
					! Estimate maxn
					nestimate = 10**(2*sys%nthresh/sys%nlevels)
					nestimate = nestimate * exp(-sum(sys%hrfactors))
					nestimate = nestimate ** (1.0/sys%nlevels)
					sys%maxnfix = ceiling(sys%nlevels * nestimate)
				end if
			end if

			if (sys%minnfix .lt. 0) then
				sys%minnfix = ceiling((sys%e_target*TO_EV - sys%delta_e)/(maxval(sys%energies)))
			end if

			if (sys%algorithm .eq. 'stochastic') then
				sys%mm%notunique = .true.
				sys%memory = 0.5 * sys%memory
			end if
			call sys%mm%initialise(sys%nlevels, sys%maxnoccs, sys%memory, sys%nthresh, sys%occprefix)
			call sys%mm%block_swap(sys%energies, sys%hrfactors, sys%mm%chunk_size, tmp)
		end if
	end subroutine sysdata_init

	subroutine sysdata_fc_compute(sys)
		class(sysdata), intent(inout)	:: sys

		integer		:: ix, nx
		real(dbl)	:: tmp, yk

		if (allocated(sys%fcfactors)) deallocate(sys%fcfactors)
		allocate(sys%fcfactors(sys%nlevels, sys%ncut+1))

		do ix=1, sys%nlevels
			yk = sys%hrfactors(ix)
			tmp = exp(-yk)
			do nx=0, sys%ncut
				sys%fcfactors(ix, nx+1) = sqrt(tmp)
				tmp = tmp*yk/real(nx+1)
			end do
		end do
	end subroutine sysdata_fc_compute

	subroutine sysdata_cutoff_compute(sys)
		class(sysdata), intent(inout)	:: sys

		integer		:: ix, nx, jx, tmpix
		real(dbl)	:: tmp
		logical		:: found

		if (allocated(sys%cutoffs)) deallocate(sys%cutoffs)
		allocate(sys%cutoffs(sys%nlevels, sys%nthresh+1))
		sys%cutoffs = 1

		do ix=1,sys%nlevels
			tmpix = 1
			tmp = log10(sys%fcfactors(ix, tmpix))
			do nx=0,sys%nthresh
				found = .false.
				do while (.not. found)
					if ((tmp .lt. -nx) .or. (tmpix .eq. sys%ncut)) then
						sys%cutoffs(ix, nx+1) = tmpix
						found = .true.
					else
						tmpix = tmpix+1
						tmp = log10(sys%fcfactors(ix, tmpix))
					end if
				end do
			end do
		end do

		if (allocated(sys%bounds)) deallocate(sys%bounds)
		allocate(sys%bounds(sys%nlevels, sys%ncut+1))
		sys%bounds = sys%nthresh+1
		do ix=1,sys%nlevels
			tmpix = 1
			inner: do nx=1,sys%nthresh+1
				do jx=tmpix,sys%ncut+1
					sys%bounds(ix, jx) = nx
				end do
				tmpix = sys%cutoffs(ix, nx)+1
				if (tmpix .eq. sys%ncut+1) exit inner
			end do inner
		end do

	end subroutine sysdata_cutoff_compute

	subroutine sysdata_free(sys)
		class(sysdata), intent(inout)	:: sys
		if (allocated(sys%energies)) deallocate(sys%energies)
		if (allocated(sys%hrfactors)) deallocate(sys%hrfactors)
		if (allocated(sys%masses)) deallocate(sys%masses)
		if (allocated(sys%fcfactors)) deallocate(sys%fcfactors)
		if (allocated(sys%V_vq_j)) deallocate(sys%V_vq_j)
		if (allocated(sys%Bvqj)) deallocate(sys%Bvqj)
		if (allocated(sys%V_jj)) deallocate(sys%V_jj)
		if (allocated(sys%cutoffs)) deallocate(sys%cutoffs)
		if (allocated(sys%bounds)) deallocate(sys%bounds)
		if (allocated(sys%energy_order)) deallocate(sys%energy_order)
		if (allocated(sys%maxoccs)) deallocate(sys%maxoccs)
		if (allocated(sys%minoccs)) deallocate(sys%minoccs)
	end subroutine sysdata_free

	real(dbl) function	compute_kfcn(sys, occs) result(kfcn)
		class(sysdata), intent(in) 					:: sys
		integer(smallint), dimension(sys%nlevels), intent(in)	:: occs

		integer :: ix
		kfcn = 1d0
		do ix = 1, sys%nlevels
			kfcn = kfcn * sys%fcfactors(ix, occs(ix)+1)
		end do
	end function compute_kfcn

	subroutine sysdata_compute_zn(sys, occs, res)
		class(sysdata), intent(inout)					:: sys
		integer(smallint), dimension(sys%nlevels), intent(in)		:: occs
		real(dbl), dimension(sys%nlevels), intent(out)	:: res

		integer		:: jx
		real(dbl)	:: kfcn, tmp

		kfcn = compute_kfcn(sys, occs)
		do jx=1, sys%nlevels
			res(jx) = 0.5d0 * sys%energies(jx) / sys%hrfactors(jx)
			tmp = real(occs(jx)) - sys%hrfactors(jx)
			res(jx) = sqrt(res(jx) * tmp * tmp)
			res(jx) = res(jx) * kfcn
		end do
	end subroutine sysdata_compute_zn

	subroutine sysdata_build_V(sys, grads, hessian)
		class(sysdata), intent(inout)					 	 	:: sys
		real(dbl), dimension(sys%natoms, 3), intent(inout)  	:: grads
		real(dbl), dimension(:, :), allocatable, intent(inout)  :: hessian

		integer								:: b_ios, g_ios, vx, qx, dummy1, dummy2
		character(len=1)					:: dummy_char1, dummy_char2
		real(dbl), dimension(sys%nlevels)	:: tmp
		real(dbl), dimension(3*sys%natoms, sys%nlevels) :: tmphess

		if (allocated(sys%masses)) deallocate(sys%masses)
		allocate(sys%masses(sys%natoms))
		if (allocated(sys%V_vq_j)) deallocate(sys%V_vq_j)
		allocate(sys%V_vq_j(sys%nlevels))
		if (allocated(sys%Bvqj)) deallocate(sys%Bvqj)
		allocate(sys%Bvqj(3*sys%natoms, sys%nlevels))

		if (sys%skip_rate) then
			sys%V_vq_j = 1.0/sqrt(sys%gamma)
		else
			open(gradfile_unit, file=sys%gradfile)
			readgrads: do vx=1,sys%natoms
				do qx=1,3
					read(gradfile_unit, *, iostat=g_ios) sys%masses(vx), grads(vx, qx)
					if (g_ios .ne. 0) then
						write(*, '(1x,a,1x,i4,1x,a,1x,i4)') 'Error reading grads, ierr=', g_ios, 'line=', vx
						exit readgrads
					end if
				end do
			end do readgrads
			close(gradfile_unit)
			write(*, '(1x,a,1x,a,1x,a,1x,i4,1x,a)') 'Read grads from', trim(adjustl(sys%gradfile)), 'with', vx-1, 'atoms'

			open(bfile_unit, file=sys%bfile)
			read(bfile_unit, *, iostat=b_ios)
			sys%V_vq_j = 0d0
			readbfile: do vx=1,sys%natoms
				do qx=1,3
					read(bfile_unit, *, iostat=b_ios) dummy_char1, dummy_char2, dummy1, dummy2, tmp(:)
					call reorder_list(sys%nlevels, tmp, sys%energy_order)
					sys%Bvqj(3*(vx-1)+qx, :) = tmp(:) / sqrt(sys%masses(vx))
					if (b_ios .ne. 0) then
						write(*, '(1x,a,1x,i4,1x,a,1x,i4)') 'Error reading B-vectors, ierr=', g_ios, 'line=', vx
						exit readbfile
					else
						sys%V_vq_j(:) = sys%V_vq_j(:) + grads(vx, qx) * tmp(:) / sys%masses(vx)
					end if
				end do
			end do readbfile
			close(bfile_unit)
			write(*, '(1x,a,1x,a,1x,a,1x,i4,1x,a)') 'Read B-vectors from', trim(adjustl(sys%bfile)), 'with', vx-1, 'atoms'

			if (sys%hessfile .ne. 'none') then
				if (allocated(hessian)) deallocate(hessian)
				allocate(hessian(3*sys%natoms, 3*sys%natoms))
				if (allocated(sys%V_jj)) deallocate(sys%V_jj)
				allocate(sys%V_jj(sys%nlevels, sys%nlevels))

				open(hessfile_unit, file=sys%hessfile)
				readhess: do vx=1,3*sys%natoms
					do qx=1,vx
						read(hessfile_unit, *, iostat=b_ios) hessian(vx, qx)
						if (b_ios .ne. 0) then
							write(*, '(1x,a,1x,i4,1x,a,1x,i4)') 'Error reading hessian, ierr=', b_ios, 'line=', vx
							exit readhess
						end if
						hessian(qx, vx) = hessian(vx, qx)
					end do
				end do readhess
				close(hessfile_unit)
				write(*, '(1x,a,1x,a,1x,a,1x,i4,1x,a)') 'Read Hessian from', trim(adjustl(sys%hessfile)), 'with', (vx-1)/3, 'atoms'

				tmphess = matmul(hessian, sys%Bvqj)
				sys%V_jj = matmul(transpose(sys%Bvqj), tmphess)

				sys%V_jj = sys%V_jj * V_CONVERT / TO_MD

			end if
		end if

		sys%V_vq_j = sys%V_vq_j * V_CONVERT / TO_S
		if (sys%debug_level .gt. 1) then
			write(*, *) 'Elements of V_vq_j (FC)'
			write(*, '(1x,a4,1x,a10,1x,a20)') 'J', 'ENERGY', 'V_VQ_J'
			do qx=1,sys%nlevels
				write(*, '(1x,i4,1x,f10.6,1x,e20.8)') qx, sys%energies(qx), sys%V_vq_j(qx)
			end do

			if ((sys%debug_level .gt. 2) .and. (sys%hessfile .ne. 'none')) then
				write(*, *) 'Elements of V_ij (HT)'
				write(*, '(1x,a4,1x,a4,1x,a20)') 'I', 'J', 'V_VQ_J'
				do qx=1,sys%nlevels
					do vx=1,qx
						write(*, '(1x,i4,1x,i4,1x,e20.8)') qx, vx, sys%V_jj(qx, vx)
					end do
				end do
			end if
		end if
	end subroutine sysdata_build_V

	subroutine sysdata_calc_gamma(sys)
		class(sysdata), intent(inout)	:: sys

		integer		:: ix
		real(dbl)	:: be, w, tmp
		sys%gamma = 0d0
		do ix=1, sys%nlevels
			w  = sys%energies(ix)
			be = exp(w / (KB*300d0)) - 1d0
			be = 1d0/be
			w  = w / PLANCK
			tmp = w * w * sys%hrfactors(ix) * (2d0 * be + 1)
			sys%gamma = sys%gamma + tmp
		end do

		sys%gamma = SQRT_8LN2 * sqrt(sys%gamma)
	end subroutine sysdata_calc_gamma

	subroutine sysdata_calc_from_file(sys, record, prefix, ntot)
		class(sysdata), intent(inout)	:: sys
		integer, intent(in)				:: record
		character(len=*), intent(in)	:: prefix
		integer, intent(out)			:: ntot

		integer	:: nrows, nchunk, ix, startrow, endrow

		call sys%mm%get_nrows(sys%nlevels, record, prefix, nrows)
		write(*, *) 'Reading from file, expecting ', nrows, ' rows'
		startrow = 1
		endrow = min(sys%mm%chunk_size, nrows)-1
		ntot = 0
		main: do while (ntot .lt. nrows)
			nchunk = endrow - startrow + 1
			call sys%mm%read_from_bin(sys%nlevels, record, sys%mm%current_block, sys%energies, sys%hrfactors, prefix, startrow, endrow)
			call sys%calculate_kic(nchunk, sys%mm%current_block, init=.false., stopix=nchunk)
			ntot = ntot + nchunk
			startrow = endrow + 1
			endrow = startrow + min(nrows - ntot, sys%mm%chunk_size) - 1
		end do main
	end subroutine sysdata_calc_from_file

	subroutine sysdata_calc_kic(sys, noccs, occs, init, stopix)
		class(sysdata), intent(inout)						:: sys
		integer(bigint), intent(in)							:: noccs
		integer(smallint), dimension(sys%nlevels, noccs), intent(in)	:: occs
		logical, intent(in)									:: init
		integer, intent(in)									:: stopix

		integer			:: nx, ix
		integer(bigint)	:: counter(20)
		real(dbl)		:: z(sys%nlevels), res(sys%nlevels), tmp
		real(dbl)		:: grads(sys%natoms, 3), sums(20), nsums(100)

		real(dbl), dimension(:, :), allocatable  :: hessian

		if (init) then
			counter = 0
			sums = 0d0
			nsums = 0d0
			sys%k_ic = 0d0
			sys%k_ic_ht = 0d0
			call sys%calculate_gamma
			if (.not. allocated(sys%V_vq_j)) call sys%build_V(grads, hessian)
		end if

		do nx=1,stopix
			call sys%compute_zn(occs(:, nx), res)
			tmp = dot_product(res, sys%V_vq_j)
			tmp = tmp**2
			sys%k_ic = sys%k_ic + tmp
			tmp = 4d0 * tmp / sys%gamma
			ix = int(min(max(1d0, log10(tmp)+15d0), 20d0))
			counter(ix) = counter(ix) + 1
			sums(ix) = sums(ix) + tmp
			ix = int(sum(occs(:, nx)))
			nsums(ix) = nsums(ix) + tmp
			if (sys%debug_level .gt. 3) then
				if (tmp > 1.0d0) then
					write(*, *) occs(:, nx), tmp
				end if
			end if
		end do

		if (sys%debug_level .gt. 2) then
			write(*, *) '\nLog10', 'Counter', 'Sum'
			do nx=1,20
				write(*, *) nx-15, counter(nx), sums(nx)
			end do
		end if

		if (sys%debug_level .gt. 1) then
			write(*, *) '\nNSUMS'
			do nx=1,100
				if (nsums(nx) .gt. 1e-12) then
					write(*, *) nx, nsums(nx)
				end if
			end do
		end if

		if (sys%hessfile .ne. 'none') then
			! Do Hertzberg-Teller calculation
			do nx=1,stopix
				call sys%compute_zn(occs(:, nx), z)
				res = matmul(sys%V_jj, z)
				do ix=1,sys%nlevels
					z(ix) = 0.5 * (occs(ix, nx) + sys%hrfactors(ix))**2
					z(ix) = z(ix) / (sys%energies(ix)*sys%hrfactors(ix))
					z(ix) = sqrt(z(ix))
				end do
				tmp = dot_product(res, z)
				tmp = tmp**2
				sys%k_ic_ht = sys%k_ic_ht + tmp
			end do
		end if

	end subroutine sysdata_calc_kic

	subroutine sysdata_calc_lagrange(sys, occs, klag)
		class(sysdata), intent(inout)					:: sys
		real(dbl), dimension(sys%nlevels), intent(in)	:: occs
		real(dbl), intent(out)							:: klag

		real(dbl) 	:: tmp, kfcn, z(sys%nlevels), res(sys%nlevels), yk, nfac
		integer		:: jx

		kfcn = 1d0
		do jx=1, sys%nlevels
			yk = sys%hrfactors(jx)
			nfac = factorial(nint(occs(jx)))
			tmp = (yk**occs(jx)) * exp(-yk) / nfac
			kfcn = kfcn * sqrt(tmp)
		end do

		do jx=1, sys%nlevels
			res(jx) = 0.5d0 * sys%energies(jx) / sys%hrfactors(jx)
			tmp = occs(jx) - sys%hrfactors(jx)
			res(jx) = sqrt(res(jx) * tmp * tmp)
			res(jx) = res(jx) * kfcn
		end do

		tmp = dot_product(res, sys%V_vq_j)
		tmp = tmp**2
		tmp = 4d0 * tmp / sys%gamma
		klag = tmp * 1.9652703588D+3

		write(*, '(1x,a,1x,e14.8,1x,a)') 'Lagrange k_ic_fc=', klag, 's-1'

		if (sys%hessfile .ne. 'none') then
			z = matmul(sys%V_jj, res)
			do jx=1, sys%nlevels
				res(jx) = 0.5 * (occs(jx) + sys%hrfactors(jx))**2
				res(jx) = res(jx) / (sys%energies(jx)*sys%hrfactors(jx))
				res(jx) = sqrt(res(jx))
			end do
			tmp = dot_product(res, z)
			tmp = tmp**2
			tmp = 4d0 * tmp / sys%gamma
			tmp = tmp * 1.9652703588D+3

			write(*, '(1x,a,1x,e14.8,1x,a)') 'Lagrange k_ic_ht=', tmp, 's-1'
		end if
		write(*, *)

	end subroutine sysdata_calc_lagrange

	subroutine sysdata_get_next_bound(sys, occ, ix, bnd)
		class(sysdata), intent(in)					:: sys
		integer(smallint), dimension(sys%nlevels), intent(in)	:: occ
		integer, intent(in)							:: ix
		integer, intent(out)						:: bnd

		integer	:: jx
		bnd = sys%nthresh + 2
		do jx=1,ix
			bnd = bnd - sys%bounds(jx, occ(jx)+1)
		end do
		bnd = max(1, bnd)

	end subroutine sysdata_get_next_bound

	subroutine sysdata_ncombinations(sys, maxnocc)
		class(sysdata), intent(in)		:: sys
		integer(bigint), intent(out)	:: maxnocc

		integer(smallint) :: occ(sys%nlevels)
		maxnocc = 0
		occ = 0
		occ(1) = sys%cutoffs(1, sys%nthresh+1)
		call next_n(sys, occ, 1, maxnocc)
	end subroutine

	recursive subroutine next_n(sys, occ, ix, maxnocc)
		type(sysdata), intent(in)						:: sys
		integer(bigint), intent(out) 					:: maxnocc
		integer, intent(in)								:: ix
		integer(smallint), dimension(sys%nlevels), intent(inout)	:: occ

		integer :: maxix, jx, bnd
		integer(bigint) :: tmpmax
		maxix = occ(ix)
		maxnocc = 1
		if (ix .eq. sys%nlevels) then
			maxnocc = maxix + 1
		else
			do jx=0,maxix
				occ(ix) = jx
				call sys%get_next_bound(occ, ix, bnd)
				occ(ix+1) = sys%cutoffs(ix+1, bnd)
				call next_n(sys, occ, ix+1, tmpmax)
				maxnocc = maxnocc + tmpmax
			end do
		end if
	end subroutine next_n
end module

max_stars_repo_name: robashaw/knapsack

max_issues_repo_name: robashaw/knapsack

max_forks_repo_name: robashaw/knapsack


text:
      subroutine sphere1(itype,icount,nrt,radius,xc,yc,zc,uradius)
C
C#######################################################################
C
C     PURPOSE -
C
C     THIS ROUTINE GENERATES A SPHERE OF A SPECIFIED RADIUS
C
C
C     FORMAT: SPHERE/ITYPE=1/NPT/RADIUS/XCEN/YCEN/ZCEN/VELOCITY
C
C
C     INPUT:
C
C           NRT - NUMBER OF POINTS ON RADIUS SPHERE (THIS WILL BE ROUNDED
C                 TO THE NEAREST SQRT(INT(NUMBER OF POINTS/6 SIDES))
C
C        RADIUS - RADIUS OF THE SPHERE
C
C            XC - X COORDINATE OF THE CENTER OF THE SPHERE
C            YC - Y COORDINATE OF THE CENTER OF THE SPHERE
C            ZC - Z COORDINATE OF THE CENTER OF THE SPHERE
C
C       URADIUS - VELOCITY (RELATIVE TO THE CENTER )
C
C
C     INPUT ARGUMENTS -
C
C
C
C     OUTPUT ARGUMENTS -
C
C
C
C     CHANGE HISTORY -
C
C        $Log: sphere1.f,v $
C        Revision 2.00  2007/11/09 20:04:03  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.1   Wed Apr 05 13:35:04 2000   nnc
CPVCS    Minor source modifications required by the Absoft compiler.
CPVCS
CPVCS       Rev 1.0   Mon Jan 31 11:14:18 2000   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.7   Mon Apr 14 17:02:04 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.6   Tue Mar 05 12:35:36 1996   het
CPVCS    Correct an error with the zero point.
CPVCS
CPVCS       Rev 1.5   Thu Feb 08 15:55:22 1996   dcg
CPVCS    replace uic,vic,wic with user defined velocity vector
CPVCS
CPVCS       Rev 1.4   11/16/95 15:22:46   dcg
CPVCS    replace character literals in calls
CPVCS
CPVCS       Rev 1.3   10/22/95 13:19:00   het
CPVCS    Fix the call to rz because the calling sequence changed
CPVCS
CPVCS       Rev 1.2   08/29/95 12:17:02   het
CPVCS    Add the cmowrk storage block for each CMO
CPVCS
CPVCS       Rev 1.1   08/15/95 18:19:50   het
CPVCS    Cleanup code and correct errors
CPVCS
CPVCS       Rev 1.0   01/17/95 16:39:34   pvcs
CPVCS     Original version
C
C
C#######################################################################
C
C
      implicit none
C
C#######################################################################
C
      include 'cmo.h'
      include 'chydro.h'
C
      character*4096 ibuff
      character*32 cvelnm,isubname
      pointer(ipout,out)
      real*8 out(*)
      pointer(ipvels,vels)
      real*8 vels(3,*)
      integer ipointj,length,icmotype,icscode,ierror,nrt,icts,nrt1,
     *  iout,lout,icount,itype,ier,i1,nrtx,nrty,nrtz,ntets,npoints,
     *  lenxic,ityp
      real*8 xc,yc,zc,rout,radius,uradius,xdiff,ydiff,zdiff,radxyz,
     * xradius,yradius,zradius
C
C#######################################################################

      isubname='shpere1'
C
      call cmo_get_name(cmo,ierror)
C
      call cmo_get_info('ipointj',cmo,
     *                ipointj,length,icmotype,icscode)
      if (icscode .ne. 0) call x3d_error(isubname,'ipointj')
C
      call cmo_get_info('xic',cmo,ipxic,lenxic,icmotype,ierror)
      length=ipointj+nrt
      if(length.gt.lenxic) then
         call cmo_set_info('nnodes',cmo,length,1,1,ierror)
         call cmo_newlen(cmo,ierror)
      endif
C
      call cmo_get_info('nnodes',cmo,npoints,length,icmotype,ierror)
      call cmo_get_info('nelements',cmo,ntets,length,icmotype,ierror)
      call cmo_get_info('mbndry',cmo,mbndry,length,icmotype,ierror)
      call cmo_get_info('xic',cmo,ipxic,length,icmotype,ierror)
      call cmo_get_info('yic',cmo,ipyic,length,icmotype,ierror)
      call cmo_get_info('zic',cmo,ipzic,length,icmotype,ierror)
C
C.......................................................................
C
      icts=ipointj+1
      nrt1=int(sqrt(float(nrt+1)/6.0))
      if(radius.le.1.0e-10) then
         ipointj=ipointj+1
         xic(ipointj)=0.0+xc
         yic(ipointj)=0.0+yc
         zic(ipointj)=0.0+zc
         call cmo_set_info('ipointj',cmo,
     *                   ipointj,1,1,icscode)
         if (icscode .ne. 0) call x3d_error(isubname,'set_ipointj')
         goto 9999
      endif
C
      nrtx=nrt1
      nrty=nrt1
      nrtz=nrt1
      xradius=radius
      yradius=radius
      zradius=radius
C
      ibuff=' '
      write(ibuff,9000) 'rz/xyz/',
     *     nrtx ,     nrty ,        1 ,
     *  -xradius , -yradius , -zradius ,
     *  +xradius , +yradius , -zradius
      call dotaskx3d(ibuff,ierror)
C
      ibuff=' '
      write(ibuff,9000) 'rz/xyz/',
     *      nrtx ,     nrty ,        1 ,
     *  -xradius , -yradius , +zradius ,
     *  +xradius , +yradius , +zradius
      call dotaskx3d(ibuff,ierror)
C
C
      ibuff=' '
      write(ibuff,9000) 'rz/xyz/',
     *      nrtx ,        1 ,     nrtz ,
     *  -xradius , -yradius , -zradius ,
     *  +xradius , -yradius , +zradius
      call dotaskx3d(ibuff,ierror)
C
      ibuff=' '
      write(ibuff,9000) 'rz/xyz/',
     *      nrtx ,        1 ,     nrtz ,
     *  -xradius , +yradius , -zradius ,
     *  +xradius , +yradius , +zradius
      call dotaskx3d(ibuff,ierror)
C
      ibuff=' '
      write(ibuff,9000) 'rz/xyz/',
     *         1 ,     nrty ,     nrtz ,
     *  -xradius , -yradius , -zradius ,
     *  -xradius , +yradius , +zradius
      call dotaskx3d(ibuff,ierror)
C
      ibuff=' '
      write(ibuff,9000) 'rz/xyz/',
     *         1 ,     nrty ,     nrtz ,
     *  +xradius , -yradius , -zradius ,
     *  +xradius , +yradius , +zradius
      call dotaskx3d(ibuff,ierror)
C
      call cmo_get_info('ipointj',cmo,
     *                icount,length,icmotype,icscode)
      if (icscode .ne. 0) call x3d_error(isubname,'get_ipointj')
C
      call cmo_get_info('nnodes',cmo,npoints,length,icmotype,ierror)
      call cmo_get_info('nelements',cmo,ntets,length,icmotype,ierror)
      call cmo_get_info('mbndry',cmo,mbndry,length,icmotype,ierror)
      call cmo_get_info('xic',cmo,ipxic,length,icmotype,ierror)
      call cmo_get_info('yic',cmo,ipyic,length,icmotype,ierror)
      call cmo_get_info('zic',cmo,ipzic,length,icmotype,ierror)
C
      do 100 i1=icts,icount
         radxyz=sqrt(xic(i1)*xic(i1)+yic(i1)*yic(i1)+zic(i1)*zic(i1))
         xic(i1)=xic(i1)*radius/radxyz+xc
         yic(i1)=yic(i1)*radius/radxyz+yc
         zic(i1)=zic(i1)*radius/radxyz+zc
 100  continue
C
      if(abs(uradius).lt.1.0e-15) goto 9999
      call cmo_get_attinfo('velname',cmo,iout,rout,cvelnm,
     *   ipout,lout,ityp,ierror)
      call cmo_get_info(cvelnm,cmo,ipvels,length,icmotype,ier)
      if(ier.ne.0) go to 9999
      do 200 i1=icts,icount
         xdiff=xic(i1)-xc
         ydiff=yic(i1)-yc
         zdiff=zic(i1)-zc
         radxyz=sqrt(xdiff*xdiff+ydiff*ydiff+zdiff*zdiff)
         vels(1,i1)=uradius*xdiff/radxyz
         vels(2,i1)=uradius*ydiff/radxyz
         vels(3,i1)=uradius*zdiff/radxyz
 200  continue
C

 9000 format(a7,3(i3,'/'),6(1pe15.7,'/'),' ; finish')
 9999 return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
module ointerface

  use userinterface
  use filters      , only: absfilter

  abstract interface
     subroutine optimization(n,y,l,u,me,mi,uevalf,uevalc,uevallc,uevalljac, &
          nf,alpha,ffilter,hfilter,filterTest,outiter,epsfeas,epsopt,       &
          verbose,delta,fy,hynorm,rho,flag)

       ! TODO: Comment this interface!

       ! SCALAR ARGUMENTS
       logical :: verbose
       integer :: flag,me,mi,n,nf,outiter
       real(8) :: alpha,delta,epsfeas,epsopt,fy,hynorm,rho

       ! ARRAY ARGUMENTS
       real(8) :: ffilter(nf),hfilter(nf),l(n),u(n),y(n)

       ! EXTERNAL SUBROUTINES
       procedure(evalf)     :: uevalf
       procedure(evalc)     :: uevallc,uevalc
       procedure(evaljac)   :: uevalljac
       procedure(absfilter) :: filterTest

       intent(in   ) :: alpha,epsfeas,epsopt,ffilter,hfilter,l,me,mi,n, &
            nf,outiter,u,verbose
       intent(out  ) :: hynorm,flag
       intent(inout) :: delta,fy,rho,y
     end subroutine optimization
  end interface

contains

end module ointerface

max_stars_repo_name: fsobral/fkss

max_issues_repo_name: fsobral/fkss

max_forks_repo_name: fsobral/fkss


text:
!======================================================================!
SUBROUTINE Split(sub, Npar)
    !----------------------------------------------------------------------!
    !   Splits the domain by a geometrical multisection technique.         !
    !----------------------------------------------------------------------!
    !------------------------------[Modules]-------------------------------!
    USE allp_mod, only : huge
    USE all_mod
    USE gen_mod
    USE div_mod
    USE par_mod
    !----------------------------------------------------------------------!
    IMPLICIT NONE
    !-----------------------------[Parameters]-----------------------------!
    INTEGER :: sub                           ! subdomain to be splitted
    INTEGER :: Npar                          ! number of new partitions
    !-------------------------------[Locals]-------------------------------!
    CHARACTER :: dir                         ! direction for splitting
    INTEGER   :: c, ic, j
    INTEGER   :: Nfiled
    REAL      :: xmax,ymax,zmax,xmin,ymin,zmin, delx,dely,delz,dxyz
    !--------------------------------[CVS]---------------------------------!
    !  $Id: Split.f90,v 1.8 2002/10/30 16:29:19 niceno Exp $
    !  $Source: /home/muhamed/.CVSROOT/T-Rex/Divide/Split.f90,v $
    !======================================================================!

    !--------------------------------------------!
    !     Find the smallest moment of inertia    !
    !--------------------------------------------!
    if(ALGOR == INERTIAL) then
        call Inertia(sub)
    end if

    !-------------------------------------------------!
    !     Find the largest dimension of the domain    !
    !-------------------------------------------------!
    if(ALGOR == COORDINATE) then
        xmax=-HUGE
        ymax=-HUGE
        zmax=-HUGE
        xmin=+HUGE
        ymin=+HUGE
        zmin=+HUGE
        do c=1,NC
            if(proces(c) == sub) then
                xmax=max(xmax, xc(c))
                ymax=max(ymax, yc(c))
                zmax=max(zmax, zc(c))
                xmin=min(xmin, xc(c))
                ymin=min(ymin, yc(c))
                zmin=min(zmin, zc(c))
            end if
        end do
        delx = xmax - xmin
        dely = ymax - ymin
        delz = zmax - zmin

        dxyz = max(delx,dely,delz)
        if(delz == dxyz) dir='z'
        if(dely == dxyz) dir='y'
        if(delx == dxyz) dir='x'
    end if

    do j=1,Npar-1

        Nfiled   = 0

        !----------------------------!
        !     Fill the subdomain     !
        !----------------------------!

        if(ALGOR==COORDINATE) then

            if(dir == 'x') then
                do c=1,NC
                    ic=ix(c)
                    if(proces(ic) == sub) then
                        proces(ic) = Nsub+j
                        Nfiled = Nfiled + 1
                        if(Nfiled >= subNC(sub)/(Npar-j+1)) goto 2
                    end if
                end do
            end if

            if(dir == 'y') then
                do c=1,NC
                    ic=iy(c)
                    if(proces(ic) == sub) then
                        proces(ic) = Nsub+j
                        Nfiled = Nfiled + 1
                        if(Nfiled >= subNC(sub)/(Npar-j+1)) goto 2
                    end if
                end do
            end if

            if(dir == 'z') then
                do c=1,NC
                    ic=iz(c)
                    if(proces(ic) == sub) then
                        proces(ic) = Nsub+j
                        Nfiled = Nfiled + 1
                        if(Nfiled >= subNC(sub)/(Npar-j+1)) goto 2
                    end if
                end do
            end if
        end if

        if(ALGOR==INERTIAL) then
            do c=1,NC
                ic=iin(c)
                if(proces(ic) == sub) then
                    proces(ic) = Nsub+j
                    Nfiled = Nfiled + 1
                    if(Nfiled >= subNC(sub)/(Npar-j+1)) goto 2
                end if
            end do
        end if

        !--------------------------------!
        !     Subdomain is filled up     !
        !--------------------------------!
2       subNC(Nsub+j) = Nfiled
        subNC(sub) = subNC(sub) - Nfiled

    end do  ! j

    Nsub = Nsub + Npar - 1

END SUBROUTINE Split

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


! { dg-do compile }
! { dg-options "-fcoarray=single" }
!
! PR fortran/44646
!
! DO CONCURRENT
!
implicit none
integer :: i, j

outer: do, concurrent ( i = 1 : 4)
  do j = 1, 5
    if (j == 1) cycle ! OK
    cycle outer ! OK: C821   FIXME
    exit outer ! { dg-error "EXIT statement at .1. leaves DO CONCURRENT construct" }
  end do
end do outer

do concurrent (j = 1:5)
  cycle ! OK
end do

outer2: do j = 1, 7
  do concurrent (j=1:5:2) ! cycle outer2 - bad: C821
    cycle outer2 ! { dg-error "leaves DO CONCURRENT construct" }
  end do
end do outer2

do concurrent ( i = 1 : 4)
  exit ! { dg-error "EXIT statement at .1. leaves DO CONCURRENT construct" }
end do
end

subroutine foo()
  do concurrent ( i = 1 : 4)
    return   ! { dg-error "Image control statement RETURN" }
    sync all ! { dg-error "Image control statement SYNC" }
    call test () ! { dg-error "Subroutine call to .test. in DO CONCURRENT block at .1. is not PURE" }
    stop ! { dg-error "Image control statement STOP" }
  end do
  do concurrent ( i = 1 : 4)
    critical ! { dg-error "Image control statement CRITICAL at .1. in DO CONCURRENT block" }
      print *, i
!    end critical
  end do

  critical
    do concurrent ( i = 1 : 4) ! OK
    end do
  end critical
end

subroutine caf()
  use iso_fortran_env
  implicit none
  type(lock_type), allocatable :: lock[:]
  integer :: i
  do, concurrent (i = 1:3)
    allocate (lock[*]) ! { dg-error "ALLOCATE of coarray at .1. in DO CONCURRENT block" }
    lock(lock) ! { dg-error "Image control statement LOCK" }
    unlock(lock) ! { dg-error "Image control statement UNLOCK" }
    deallocate (lock) ! { dg-error "DEALLOCATE of coarray at .1. in DO CONCURRENT block" }
  end do

  critical
    allocate (lock[*]) ! { dg-error "ALLOCATE of coarray at .1. in CRITICAL block" }
    lock(lock) ! { dg-error "Image control statement LOCK" }
    unlock(lock) ! { dg-error "Image control statement UNLOCK" }
    deallocate (lock) ! { dg-error "DEALLOCATE of coarray at .1. in CRITICAL block" }
  end critical
end subroutine caf


text:
	SUBROUTINE FL_FLUN  ( lun, iret )
C************************************************************************
C* FL_FLUN								*
C* 									*
C* This subroutine frees a logical unit number that was allocated by 	*
C* FL_GLUN.  A logical unit number should be freed when it is no 	*
C* longer needed.							*
C* 									*
C* FL_FLUN  ( LUN, IRET )						*
C* 									*
C* Input parameters:							*
C*	LUN		INTEGER		Logical unit number 		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					   0 = normal return		*
C**									*
C* Log:									*
C* G.C.Chatters/RDS	 3/84						*
C* M. desJardins/GSFC	 3/86	Changed comments			*
C* M. desJardins/GSFC	 7/90	Eliminate system service		*
C* M. desJardins/NMC	 8/94	Documentation				*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
C------------------------------------------------------------------------
C*	Free logical unit number.
C
	IF  ( ( lun .ge. 11 ) .and. ( lun .le. 20 ) )  THEN
	    lungem ( lun - 10 ) = 0
	END IF
C*
	iret = 0
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak


text:
! example of simple Fortran AMD GPU offloading
program main
  parameter (nsize=1000000)
  real a(nsize), b(nsize), c(nsize)
  integer i

  do i=1,nsize
    a(i)=0
    b(i) = i
    c(i) = 10
  end do

  call foo(a,b,c)

  write(6,*)"a(1)=", a(1), "    a(2)=", a(2)
  if (a(1).ne.11 .or. a(2).ne.22) then
    stop "ERROR: wrong answers"
  endif
  write(6,*)"Success: if a diagnostic line starting with DEVID was output"
  return
end
subroutine foo(a,b,c)
  parameter (nsize=1000000)
  real a(nsize), b(nsize), c(nsize)
  integer i
!$omp target map(from:a) map(to:b,c)
!$omp parallel do
  do i=1,nsize
    a(i) = b(i) * c(i) + i
  end do
!$omp end target
  return
end

max_stars_repo_name: illuhad/aomp

max_issues_repo_name: illuhad/aomp

max_forks_repo_name: illuhad/aomp


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:27 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE IS_LEAP_YEAR__genmod
          INTERFACE
            FUNCTION IS_LEAP_YEAR(YR)
              INTEGER(KIND=4) :: YR
              LOGICAL(KIND=4) :: IS_LEAP_YEAR
            END FUNCTION IS_LEAP_YEAR
          END INTERFACE
        END MODULE IS_LEAP_YEAR__genmod


text:
Michael Corbett served on the City Council from 6/20/1986 to 6/20/1990. During that time, he served as Mayor from 1988 to 1989. He and Judy Corbett are responsible for Village Homes. Michael and Judy separated, and Michael currently resides in downtown Davis. Michael was the principal designer of the Covell Village project. He is one of the authors of the Ahwahnee Principles.  He was responsible for expanding Central Park for the future expansion of the Farmers Market. He never swayed on his commitment to maintaining the integrity of Davis downtown keeping the theatres from going to a freeway location. He voted to get the city to build the theater in the core of downtown on city property. One of the many choices he made as Mayor to keep the Village concept strong. He also had the drainage ponds in North and west Davis converted to wildlife habitat.

Michael and former Davisite Users/JohnOliver John Oliver were very involved in Measure S in 1986.  That was a citizens referendum to reverse a city council decision to sell the ArdenMayfair site between 3rd and 4th streets for commercial development (to a company called Terranomics).  That site is now the piece of Central Park that is home to the Farmers Market and the US Bicycling Hall of Fame.

You can read about Measure S at http://www.city.davis.ca.us/pb/cultural/30years/chapt06.cfm the city website.

20070625 16:26:57 nbsp Micheal Corbett was the sole designer and developer of the world famous Village Homes in Davis, CA. Based on his lifelong interest in architecture and the environment, he developed one of the first solar and environmental communities in the world. He won acclaim not only from our President at the time Jimmy Carter, but from many other world leaders. Many concepts from his first book A Better Place to Live were used as he forged ahead to build a community many at the time thought was a little offbeat! He also designed the most recent environmental development, Covell Village which was voted down by the citizens of Davis. As Davis claims to be so concerned about growth, and environmental issues they were put to the test when this development failed. Mace Ranch and Wildhorse Davis first production sprawl developments must have caught on. Having one of the most acclaimed developments in the world was more than the citizens could handle. Users/Louis

I have a hard time with the form and substance of Louis posting; following the link to Covell Village will show anyone with an interest that there are many perspectives & arguments about the environmental impact of the failed/foiled development. Id say Mike was put to the test, and his environmental credentials couldnt cover up the traffic and financial costs the project would have had. But thats just another opinion. Users/DougWalter

20070626 09:16:09 nbsp While agree with a portion of Users/Louiss comments regarding the Covell Village project, and the citizens desire to not be at the forefront of sustainable development, the overall fawning tone makes me uncomfortable, like I should genuflect or somthing when Mr. Corbetts name is mentioned.   Users/RocksandDirt

Mike IS an influential thinker and developer. Hes conscious of it, but  then again  hes a very approachable person. Thats the combination that helped him succeed in politics. Theres a lot to debate about issues like sprawl, infill, development on agricultural land, and creating community; the character of individuals is bound to seep into that, but theres no reason that personalities need to become the focus. Users/DougWalter

  I agree that personalities should not be the focus in these debates and almost always get sucked into them.  It actually hurts the quality of the debate substantially.  I actually feel (and this is totally my own un substantiated opinion) that the fawning tone most people take towards Mr. Corbett and his ideas hurt the Covell Village project.  Davis (like most towns) has a built in vote against personality in these situations, and I think it hurt the project. Users/RocksandDirt

20070626 22:20:16 nbsp I was a supporter of Covell Village and have always admired Michel Corbett, but Louis seems to think he is the next Saint Francis of Assis. Nature is under good stewardship as long as Mike stays around.  Users/Jeffreymh

20080109 15:33:23 nbsp Ive known Michael Corbett for some years now; he is a man who loves Davis very much, and wanted to give something back to it, in the way he knew how:  by designing an environmentally sustainable community within Davis, called Covell Village.  The plans I saw were beautiful; he could see into the future at what that lifestyle within a neighborhood be like, and probably would have lived there himself if it had been approved.  It had the walking and bike paths, senior center and garden sites, school and parks, water/wetlands areas.  Davis is going to grow no matter what people want; but I guess what they want is more Mace Ranch monoliths that look like every other subdivision in California.  Can you say SOUND WALL? Users/Beahah

20080109 16:02:06 nbsp Mike was also involved in the planning of Manzanita, which would have been a small city of several thousand people just over the county line from Winters in Solano County. Touted for its environmentally friendly features, Manzanita triggered a countywide vote and was defeated by Solano County voters in 1985. The proposal was partly responsible for Solano Countys Measure A, which restricted growth to the existing cities.  Users/DonShor

20081214 22:04:32 nbsp I think the Covell Village proposal was just too big.  People had a hard time with the size.  Regardless of how trustworthy the developers were. Users/JeffShaw

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE GDTGRF  ( iytype, ratio, xstrt, xstop, ystrt, ystop,
     +			     rmargn, iret )
C************************************************************************
C* GDTGRF								*
C*									*
C* This subroutine sets the graph margins and sets up the graph in	*
C* GEMPLT.  The margins are set to 5, 4, 5, 3.				*
C*									*
C* GDTGRF  ( IYTYPE, RATIO, XSTRT, XSTOP, YSTRT, YSTOP, RMARGN, IRET )	*
C*									*
C* Input parameters:							*
C*	IYTYPE		INTEGER		Coordinate system		*
C*					  1 = linear			*
C*					  2 = log			*
C*					  3 = ** KAPPA			*
C*					  4 = log with skew x		*
C*	RATIO		REAL		Height to width plot ratio	*
C*	XSTRT		REAL		Left x				*
C*	XSTOP		REAL		Right x				*
C*	YSTRT		REAL		Bottom y			*
C*	YSTOP		REAL		Top y				*
C*	RMARGN (4)	REAL		Margins				*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -8 = error setting graph	*
C**									*
C* Log:									*
C* G. Huffman/GSC	 2/89	Adapted from GDPGRF			*
C* S. Schotz/GSC	 7/90	Update for rmargn array			*
C* K. Brill/NMC         01/92   Replace GERROR with ER_WMSG             *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	REAL		rmargn (*)
C----------------------------------------------------------------------
	iret = 0
C
C*	Set graph margins, with a default of 5, 4, 5, 3 (characters).
C
	IF  ( rmargn (1) .lt. 0 ) THEN
	    CALL GSGMGN  ( 5., 4., 5., 3., ier )
	  ELSE
	    CALL GSGMGN  ( rmargn (1), rmargn (2), rmargn (3),
     +			   rmargn (4), ier )
	END IF
C
C*	Set graph.
C
	jxtype = 1
	jytype = iytype
	IF  ( iytype .eq. 4 )  THEN
	    jxtype = 4
	    jytype = 2
	END IF
	CALL GSGRAF  ( jxtype, jytype, ratio, xstrt, ystrt, xstop,
     +		       ystop,  iret )
	IF  ( iret .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'GEMPLT', iret, ' ', ier )
	    iret = -8
	    CALL ER_WMSG  ( 'GDTSER', iret, ' ', ier )
	END IF
C*
	RETURN
	END

max_forks_repo_name: oxelson/gempak


text:
$ Declare
   SUBROUTINE FactorKr ( a, c, B, W, theta, Ym, Yb, Spm, Spb, Ssm, Ssb, K1C, Lr, Kr )
!
! Parameters
   IMPLICIT NONE
   DOUBLE PRECISION, INTENT(IN)    :: a            ! Crack depth
   DOUBLE PRECISION, INTENT(IN)    :: c            ! Crack width
   DOUBLE PRECISION, INTENT(IN)    :: B            ! Plate thickness
   DOUBLE PRECISION, INTENT(IN)    :: W            ! Plate width
   DOUBLE PRECISION, INTENT(IN)    :: theta        ! FT = 2: Angle in radians (0 => c, pi/2 => a)
   DOUBLE PRECISION, INTENT(IN)    :: Ym           ! Geometric Correction Factor for membrane stress
   DOUBLE PRECISION, INTENT(IN)    :: Yb           ! Geometric Correction Factor for bending stress
   DOUBLE PRECISION, INTENT(IN)    :: Spm          ! primary membrane stress
   DOUBLE PRECISION, INTENT(IN)    :: Spb          ! primary bending stress
   DOUBLE PRECISION, INTENT(IN)	   :: Ssm          ! secondary membrane stress
   DOUBLE PRECISION, INTENT(IN)	   :: Ssb          ! secondary bending stress
   DOUBLE PRECISION, INTENT(IN)	   :: K1C          ! fracture toughness
   DOUBLE PRECISION, INTENT(IN)    :: Lr		   ! plasticity factor
   DOUBLE PRECISION, INTENT(OUT)   :: Kr		   ! fracture factor
!
! Declarations
   DOUBLE PRECISION                :: KIp, KIs, KI, chi, rho1, rho
   DOUBLE PRECISION                :: Mm, Mb, Mmmax
!
! Constants
   DOUBLE PRECISION :: pi
   pi = 3.14159265358979D0
!
! Calculate Kr
	  CALL SIMFactorMm ( a, c, B, W, theta, Mm, Mmmax)
	  CALL SIMFactorMb ( a, c, B, W, theta, Mb)
   	  KIp = (Ym * Spm + Yb * Spb) * DSQRT(pi * a)
      KIs = (Mmmax * Ssm + Mb * Ssb) * DSQRT(pi * a)
	  KI = KIp + KIs
	  chi = (KIs * Lr / KIp)
	  rho1 = 0.1D0 * (chi)**0.714D0 - 0.007D0 * chi * chi + 3.D-5 * (chi)**5
	  IF (chi .GT. 4.D0) THEN
		 STOP 'Error in Subroutine FactorKr: KIs*Lr/KIp >4'
	  ELSE IF (Lr .LE. 0.8D0) THEN
		 rho = rho1
	  ELSE IF (Lr .LT. 1.05D0) THEN
		 rho = 4.D0 * rho1 * (1.05D0 - Lr)
	  ELSE
		 rho = 0.D0
	  END IF
	  Kr = KI/K1C + rho
!
! End subroutine
 RETURN
 END

max_stars_repo_name: monicanicoreac/ProRail-ravelingen

max_issues_repo_name: monicanicoreac/ProRail-ravelingen

max_forks_repo_name: monicanicoreac/ProRail-ravelingen


!==============================================================================
! CALCULO DE LOS COEFICIENTES DE VAN DER WAALS A PARTIR DE LOS ARCHIVOS DE
! LOS ARCHIVOS DE LA G(r) OBTENIDOS DE LA SIMULACION CON POTENCIAL SW
!
! AUTOR : MARTIN ALEJANDRO PAREDES SOSA
!==============================================================================

!DECLARACION DE VARIABLES
Program Waals
  !Use Basic
  Implicit None
  Integer,Parameter :: Part = 600                 ! PARA NOMBRE DE ARCHIVO
  Real, Parameter :: Dens = 0.5                ! CONCENTRACION
  Integer :: State                                ! ESTADO DE LECTURA
  Integer :: k, i, j                              ! CONTADOR
  Character (len=3), Parameter :: Start = "gdr"   ! NOMBRE DE ARCHIVO DE ENTRADA
  Character (len=4), Parameter :: En = ".dat"     ! EXTENSION ARCHIVO DE ENTRADA
  Character (len=10):: Filename, cons             ! NOMBRE DE ARCHIVO
  Real, Parameter :: PI = 4.0 * ATAN(1.0)         ! VALOR DE PI
  Real, Parameter :: Rc = 17.3205090              ! RADIO DE CORTE
  Real, Dimension(:),Allocatable ::  R , G        ! RADIO | DISTRIBUCION RADIAL
  Real :: a                                       ! TOTAL
  Real :: cte, delta                              ! PARAMETROS PARA CALCULO CTE
  Real :: Intg                                    ! ACUMULADOR PARA INTEGRACION

  Write(*,*) "========================"

  !TAMANO DEL ARCHIVO POR LEER
  Write(Cons,256) Part
  Filename = start//trim(Cons)//En
  Write(*,*) "Archivo: ",Filename

  Open( 1, File = Trim(Filename), action= "read", Status ="old" ) !ARCHIVO DE ENTRADA
  k = 0 !REINICIA CONTADOR
  Sizes: Do                         !BUSCANDO TAMANO DE ARCHIVO (RENGLONES QUE QUE TIENE)

     Read( 1,*, iostat = state  )
     k = k + 1
     If ( state .LT. 0 ) Exit

  End Do Sizes

  Write(*,*) "Tiene", k, "Renglones" !DEBUG LINE (SIZE OF FILE)

  Rewind 1 !REINICIAR ARCHIVO DE ENTRADA
  Allocate ( R(k), G(k) )

  !SAVING FILE DATA
  Saves : Do i = 1, k+1

     Read( 1,*, iostat = state  ) R(i), G(i)
     If ( state .LT. 0 ) Exit

  End Do Saves
  Close(1)

  Write(*,*) "DATOS GUARDADOS EN MEMORIA"

  !NUMERO DE PARTICULAS
  Delta = 0.05                  ! CAMBIAR SEGUN EL ARCHIVO
  cte = 8.0*Dens

  ! ALGORITMO DE INTEGRACION POR TRAPECIO
  i=1
  Intg = R(i)*G(i)*0.5
  Integrando: Do j = i+1 ,k

     If ( R(j+1) .GT. RC) Exit
     Intg = Intg + R(j) * G(j)

  End Do Integrando

  Intg = Intg + 0.5*R(j)*G(j)

  a = Intg*delta*cte         ! TOTAL DE INTEGRACION

  Write(*,*) "a*= ", a
  !Write(8,*) Dens*0.1, a  , b
  !Write(9,*) Dens*0.1, a74, b

  Deallocate(R,G)
  Write(*,*) "========================"

512 Format (I5.5)
256 Format (I3.3)
End Program


text:
      function init_sequence()
      include 'common.h'
      integer init_sequence
      integer run_a_while
      integer iphase
      integer res

c     this is done by amuse_input call because it overwrites
c     AMUSE settings, so only explicit use please..
c     call input

      if (istart.eq.1) then
          write(6,*)   '  calling start..'
          call flush(6)

          call start
          iphase = 1
      else
          write(6,*)   '  not calling start..'
          call flush(6)

          call mydump(2)
c         hmmmm...
c         call input
          iphase = 2
      endif

c     only in evolve...
c     res = run_a_while(iphase)
      init_sequence = res
      end function

      function run_a_while(iphase)
      include 'common.h'
      integer evolve_step
      integer run_a_while
      integer iphase
      integer res, amuse_output

      if (iphase.ne.1) then
          call zone
          call relaxt
          iphase = 1
          res = -20
      endif

      if (iphase.eq.1) then
c     call output
c     we use our own 'cos the original stops, amuse doesnt
c     like that:
          res = amuse_output()
c     if (res.lt.0) exit
c     let python script decide when to stop...
          call zone
          call relaxt
      endif

      run_a_while = res
      end function

      function parameter_test(times)
      include 'common.h'
      integer parameter_test
      double precision times

      times = TIME
      parameter_test = 0
      end function

      function evolve_src(time_end)
      include 'common.h'
      integer evolve_src
      integer iphase
      integer res, dumpres
      double precision tcrit_, time_end

      res = 0

   10 if (time_end.gt.timet) then
          write(6,*) 'call zone'
          call zone
          write(6,*) 'call relaxt', timet
          call relaxt
          goto 10
      endif
      time_end = timet
      evolve_src = 0

      end function

      subroutine amuse_set_mmc_data_directory(datadir_)
        character(len=200) datadir
        common /AMUSE/ datadir
        character(len=200), intent(in) :: datadir_
        datadir=datadir_
      end subroutine

      subroutine amuse_input_src()
        include 'common.h'
        call input()
      end subroutine

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text: #redirect Pothole Project


text:
! For Example:
! The change in number density of C12 is
! d(n12)/dt = - 2 * 1/2 (n12)**2 <sigma v>
!
! where <sigma v> is the average of the relative velocity times the cross
! section for the reaction, and the factor accounting for the total number
! of particle pairs has a 1/2 because we are considering a reaction involving
! identical particles (see Clayton p. 293).  Finally, the -2 means that for
! each reaction, we lose 2 carbon nuclei.
!
! The corresponding Mg24 change is
! d(n24)/dt = + 1/2 (n12)**2 <sigma v>
!
! note that no factor of 2 appears here, because we create only 1 Mg nuclei.
!
! Switching over to mass fractions, using n = rho X N_A/A, where N_A is
! Avagadro's number, and A is the mass number of the nucleon, we get
!
! d(X12)/dt = -2 *1/2 (X12)**2 rho N_A <sigma v> / A12
!
! d(X24)/dt = + 1/2 (X12)**2 rho N_A <sigma v> (A24/A12**2)
!
! these are equal and opposite.
!
! The quantity [N_A <sigma v>] is what is tabulated in Caughlin and Fowler.

! we will always refer to the species by integer indices that come from
! the network module -- this makes things robust to a shuffling of the
! species ordering

subroutine rhs(n, t, y, ydot, rpar, ipar)

  use bl_types
  use bl_constants_module
  use network
  use network_indices
  use rpar_indices

  implicit none
  !$acc routine seq
  !$acc routine(screenz) seq

  ! our convention is that y(1:nspec) are the species (in the same
  ! order as defined in network.f90, and y(nspec+1) is the temperature
  integer,         intent(in   ) :: n, ipar
  real(kind=dp_t), intent(in   ) :: y(n), t
  real(kind=dp_t), intent(  out) :: ydot(n)
  real(kind=dp_t), intent(inout) :: rpar(:)

  integer :: k
!  real(kind=dp_t) :: ymass(nspec)

  !real(kind=dp_t) :: dens
  !real(kind=dp_t) :: temp, T9, T9a, dT9dt, dT9adt

  real(kind=dp_t) :: capture_rate, decay_rate

  real(kind=dp_t), PARAMETER :: &
                     one_twelvth = 1.0d0/12.0d0, &
                     five_sixths = 5.0d0/ 6.0d0, &
                       one_third = 1.0d0/ 3.0d0, &
                      two_thirds = 2.0d0/ 3.0d0

  !dens = rpar(irp_dens)
  !temp = rpar(irp_temp)

  ! compute the molar fractions -- needed for the screening
  !ymass(ic12_) = y(1)/aion(ic12_)
  !ymass(io16_) = X_O16/aion(io16_)
  !ymass(img24_) = (ONE - y(1) - X_O16)/aion(img24_)

  ! compute some often used temperature constants
  !T9     = temp/1.e9_dp_t
  !dT9dt  = ONE/1.e9_dp_t
  !T9a    = T9/(1.0e0_dp_t + 0.0396e0_dp_t*T9)
  !dT9adt = (T9a / T9 - (T9a / (1.0e0_dp_t + 0.0396e0_dp_t*T9)) * 0.0396e0_dp_t) * dT9dt

  capture_rate = 0.1_dp_t
  decay_rate   = 0.2_dp_t

  ydot(ine23_) = -decay_rate*y(ine23_) + capture_rate*y(ina23_)
  ydot(ina23_) = -capture_rate*y(ina23_) + decay_rate*y(ine23_)

  return

end subroutine rhs

max_stars_repo_name: zingale/ode-opanacc

max_issues_repo_name: zingale/ode-opanacc

max_forks_repo_name: zingale/ode-opanacc


text:
C
C	$Id: wmseti.f,v 1.18 2010-01-05 03:52:14 fred Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE WMSETI (CNP,IVP)
C
      SAVE
C
      CHARACTER*(*) CNP
C
C  This subroutine is called to give an integer value to a specified
C  parameter.
C
C  CNP is the name of the parameter whose value is to be set.
C
C  IVP is an integer variable containing the desired value.
C
C
C  Declare a local character variable in which to form an error message.
C
      CHARACTER*80 CTM
C
      include 'wmcomn.h'
C
C  Do a call forcing a BLOCKDATA to be loaded from a binary library.
C
      CALL WMBLDA
C
C  Check for an uncleared prior error.
C
      IF (ICFELL('WMSETI - Uncleared prior error',1) .NE. 0) RETURN
C
C  Check for a parameter name that is too short.
C
      IF (LEN(CNP) .LT. 3) THEN
        CTM(1:36) = 'WMSETI - Parameter name too short - '
        CTM(37:36+LEN(CNP)) = CNP
        CALL SETER (CTM(1:36+LEN(CNP)), 1, 1)
        GO TO 120
      ENDIF
C
C  Do a crude check to see if IVP is within range (primarily to
C  make sure that a floating point number has not been supplied).
C
      IF (IVP .GT. 10000000) THEN
        CTM(1:42) = 'WMSETI - Parameter value out of range for '
        CTM(43:45) = CNP(1:3)
        CALL SETER (CTM(1:45), 2, 1)
        GO TO 120
      ENDIF
C
C  MXS - Maximum number of symbols to allow on a front line.
C
      IF (CNP(1:3).EQ.'MXS' .OR. CNP(1:3).EQ.'mxs' .OR.
     +    CNP(1:3).EQ.'Mxs') THEN
        IF (IVP .GT. ISDIM) THEN
          CTM(1:42) = 'WMSETI - Parameter value out of range for '
          CTM(43:45) = CNP(1:3)
          CALL SETER (CTM(1:45), 2, 1)
          GO TO 120
        ENDIF
        MAXSYM = IVP
        GO TO 120
C
C  NMS - force a certain number of symbols on a front line (if > 0).
C
      ELSE IF (CNP(1:3).EQ.'NMS' .OR. CNP(1:3).EQ.'nms' .OR.
     +    CNP(1:3).EQ.'Nms') THEN
        NUMSYO = IVP
        GO TO 120
C
C  PAI - current parameter array index.
C
      ELSE IF (CNP(1:3).EQ.'PAI' .OR. CNP(1:3).EQ.'pai' .OR.
     +    CNP(1:3).EQ.'Pai') THEN
        IF (IVP .GT. ISDIM) THEN
          CTM(1:42) = 'WMSETI - Parameter value out of range for '
          CTM(43:45) = CNP(1:3)
          CALL SETER (CTM(1:45), 2, 1)
          GO TO 120
        ENDIF
        IARNDX = IVP
        GO TO 120
C
C  WTY - flag to indicate how line widths are to be implemented
C        (0=use GKS linewidths; 1=use internal routines WMDRFL).
C
      ELSE IF (CNP(1:3).EQ.'WTY' .OR. CNP(1:3).EQ.'wty' .OR.
     +    CNP(1:3).EQ.'Wty') THEN
        IF (ABS(IVP) .GE. 1) THEN
          IWDTYP = 1
        ELSE
          IWDTYP = 0
        ENDIF
        GO TO 120
C
C  STY - symbol type (uses PAI to determine position on line).
C
      ELSE IF (CNP(1:3).EQ.'STY' .OR. CNP(1:3).EQ.'sty' .OR.
     +    CNP(1:3).EQ.'Sty') THEN
        IF (ABS(IVP).GT.2 .OR. IVP.EQ.0) THEN
          CTM(1:42) = 'WMSETI - Parameter value out of range for '
          CTM(43:45) = CNP(1:3)
          CALL SETER (CTM(1:45), 2, 1)
          GO TO 120
        ENDIF
        ISTYPE(IARNDX) = IVP
        GO TO 120
C
C  SLF - slope flag (0=use SL1 & SL2; 1=use SL1 only; 2=use SL2 only;
C        3=use neither SL1 or SL2).
C
      ELSE IF (CNP(1:3).EQ.'SLF' .OR. CNP(1:3).EQ.'slf' .OR.
     +    CNP(1:3).EQ.'Slf') THEN
        IF (ABS(IVP).GT.3 .OR. IVP.LT.0) THEN
          ISLFLG = 3
        ELSE
          ISLFLG = IVP
        ENDIF
        GO TO 120
C
C  ALO - flag for surface or aloft.
C
      ELSE IF (CNP(1:3).EQ.'ALO' .OR. CNP(1:3).EQ.'alo' .OR.
     +    CNP(1:3).EQ.'Alo') THEN
        IF (IVP .NE. 0) THEN
          IALOFT = 1
        ELSE
          IALOFT = 0
        ENDIF
        GO TO 120
C
C  COL - specify color index to be used for all objects that require
C        only a single color setting.
C
      ELSE IF (CNP(1:3).EQ.'COL' .OR. CNP(1:3).EQ.'col' .OR.
     +    CNP(1:3).EQ.'Col') THEN
        ICOLOR = IVP
        GO TO 120
C
C  CBC - specify color index to be used for backgrounds of city and
C        dalily high/low labels.
C
      ELSE IF (CNP(1:3).EQ.'CBC' .OR. CNP(1:3).EQ.'cbc' .OR.
     +    CNP(1:3).EQ.'Cbc') THEN
        IBGCTY = IVP
        GO TO 120
C
C  RFC - specify color index to be used for foreground of regional
C        temperature labels and cities.
C
      ELSE IF (CNP(1:3).EQ.'RFC' .OR. CNP(1:3).EQ.'rfc' .OR.
     +    CNP(1:3).EQ.'Rbc') THEN
        IFGTRG = IVP
        GO TO 120
C
C  RLS - shadow color index for regional temperature labels.
C
      ELSE IF (CNP(1:3).EQ.'RLS' .OR. CNP(1:3).EQ.'rls' .OR.
     +    CNP(1:3).EQ.'Rls') THEN
        IRLLSC = IVP
        GO TO 120
C
C  EZF - EZMAP flag.  Flags if the NCL entries WMBARB, WMDRFT, and
C        WMLABS are being used in conjunction with EZMAP, in which
C        case the input coordinates are assumed to be in degrees
C        latitude and longitude.
C
      ELSE IF (CNP(1:3).EQ.'EZF' .OR. CNP(1:3).EQ.'ezf' .OR.
     +    CNP(1:3).EQ.'Ezf') THEN
        IEZFLG = IVP
        GO TO 120
C
C  ROS - outline color index for regional temperature labels.
C
      ELSE IF (CNP(1:3).EQ.'ROS' .OR. CNP(1:3).EQ.'ros' .OR.
     +    CNP(1:3).EQ.'Ros') THEN
        IRLOUC = IVP
        GO TO 120
C
C  RBS -  color index for background box for regional temperature labels.
C
      ELSE IF (CNP(1:3).EQ.'RBS' .OR. CNP(1:3).EQ.'rbs' .OR.
     +    CNP(1:3).EQ.'Rbs') THEN
        IRLBKC = IVP
        GO TO 120
C
C  DTC - specify color index to be used for the dots marking cities.
C
      ELSE IF (CNP(1:3).EQ.'DTC' .OR. CNP(1:3).EQ.'dtc' .OR.
     +    CNP(1:3).EQ.'dtc') THEN
        IDOTCO = IVP
        GO TO 120
C
C  DBC - specify color background color for dot symbols.
C
      ELSE IF (CNP(1:3).EQ.'DBC' .OR. CNP(1:3).EQ.'dbc' .OR.
     +    CNP(1:3).EQ.'dbc') THEN
        IDOTBG = IVP
        GO TO 120
C
C  NBZ - number of points to use in Bezier curves for warm front
C        symbols.
C
      ELSE IF (CNP(1:3).EQ.'NBZ' .OR. CNP(1:3).EQ.'nbz' .OR.
     +    CNP(1:3).EQ.'Nbz') THEN
        NPTSBZ = MAX(IVP,7)
        GO TO 120
C
C  REV - Reverses the directions that the symbols point.
C
      ELSE IF (CNP(1:3).EQ.'REV' .OR. CNP(1:3).EQ.'rev' .OR.
     +    CNP(1:3).EQ.'Rev') THEN
        DO 10 I=1,MAXSYM
          ISTYPE(I) = -ISTYPE(I)
   10   CONTINUE
        GO TO 120
C
C  WBF - Flag indicating whether the base of a wind barb should
C        be drawn to allow for the sky cover circle at its base.
C
      ELSE IF (CNP(1:3).EQ.'WBF' .OR. CNP(1:3).EQ.'wbf' .OR.
     +    CNP(1:3).EQ.'Wbf') THEN
        IWBBAS = IVP
        GO TO 120
C
C  WDF - Flag indicating wind barb direction.  The default (=0)
C        means plot the barb shaft away from the wind (which is
C        is not the meteorological convention) and WDF non-zero
C        means plot the barb shaft in the direction of the wind
C        (the meteorological convention).
C
      ELSE IF (CNP(1:3).EQ.'WDF' .OR. CNP(1:3).EQ.'wdf' .OR.
     +    CNP(1:3).EQ.'Wdf') THEN
        IWBDIR = IVP
        GO TO 120
C
C  SC1 - color index for center of sun symbol.
C
      ELSE IF (CNP(1:3).EQ.'SC1' .OR. CNP(1:3).EQ.'sc1' .OR.
     +    CNP(1:3).EQ.'Sc1') THEN
        ISUNC1 = IVP
        GO TO 120
C
C  SC2 - color index for the points of the sun symbol.
C
      ELSE IF (CNP(1:3).EQ.'SC2' .OR. CNP(1:3).EQ.'sc2' .OR.
     +    CNP(1:3).EQ.'Sc2') THEN
        ISUNC2 = IVP
        GO TO 120
C
C  SC3 - color index for the outline of the sun symbol.
C
      ELSE IF (CNP(1:3).EQ.'SC3' .OR. CNP(1:3).EQ.'sc3' .OR.
     +    CNP(1:3).EQ.'Sc3') THEN
        ISUNC3 = IVP
        GO TO 120
C
C  SC4 - color index for shadow of the sun symbol.
C
      ELSE IF (CNP(1:3).EQ.'SC4' .OR. CNP(1:3).EQ.'sc4' .OR.
     +    CNP(1:3).EQ.'Sc4') THEN
        ISUNC4 = IVP
        GO TO 120
C
C  CC1 - color index for interior of the cloud symbol.
C
      ELSE IF (CNP(1:3).EQ.'CC1' .OR. CNP(1:3).EQ.'cc1' .OR.
     +    CNP(1:3).EQ.'Cc1') THEN
        ICLDC1 = IVP
        GO TO 120
C
C  CC2 - color index for the outline of the cloud symbol.
C
      ELSE IF (CNP(1:3).EQ.'CC2' .OR. CNP(1:3).EQ.'cc2' .OR.
     +    CNP(1:3).EQ.'Cc2') THEN
        ICLDC2 = IVP
        GO TO 120
C
C  CC3 - color index for shadow of the cloud symbol.
C
      ELSE IF (CNP(1:3).EQ.'CC3' .OR. CNP(1:3).EQ.'cc3' .OR.
     +    CNP(1:3).EQ.'Cc3') THEN
        ICLDC3 = IVP
        GO TO 120
C
C  LC1 - color index for the interior of the lightening bolt symbol.
C
      ELSE IF (CNP(1:3).EQ.'LC1' .OR. CNP(1:3).EQ.'lc1' .OR.
     +    CNP(1:3).EQ.'Lc1') THEN
        ILTNC1 = IVP
        GO TO 120
C
C  LC2 - color index for the outline of the lightening bolt symbol.
C
      ELSE IF (CNP(1:3).EQ.'LC2' .OR. CNP(1:3).EQ.'lc2' .OR.
     +    CNP(1:3).EQ.'Lc2') THEN
        ILTNC2 = IVP
        GO TO 120
C
C  LC3 - color index for the shadow of the lightening bolt symbol.
C
      ELSE IF (CNP(1:3).EQ.'LC3' .OR. CNP(1:3).EQ.'lc3' .OR.
     +    CNP(1:3).EQ.'Lc3') THEN
        ILTNC3 = IVP
        GO TO 120
C
C  WFC - color index for warm front symbols.
C
      ELSE IF (CNP(1:3).EQ.'WFC' .OR. CNP(1:3).EQ.'wfc' .OR.
     +    CNP(1:3).EQ.'Wfc') THEN
        IWARMC = IVP
        GO TO 120
C
C  CFC - color index for cold front symbols.
C
      ELSE IF (CNP(1:3).EQ.'CFC' .OR. CNP(1:3).EQ.'cfc' .OR.
     +    CNP(1:3).EQ.'Cfc') THEN
        ICOLDC = IVP
        GO TO 120
C
C  T1C - color index for alternate colors of tropical front.
C
      ELSE IF (CNP(1:3).EQ.'T1C' .OR. CNP(1:3).EQ.'t1c' .OR.
     +    CNP(1:3).EQ.'T1c') THEN
        ITRO1C = IVP
        GO TO 120
C
C  T2C - another color index for alternate colors of tropical front.
C
      ELSE IF (CNP(1:3).EQ.'T2C' .OR. CNP(1:3).EQ.'t2c' .OR.
     +    CNP(1:3).EQ.'T2c') THEN
        ITRO2C = IVP
        GO TO 120
C
C  HIS - shadow color for high symbols.
C
      ELSE IF (CNP(1:3).EQ.'HIS' .OR. CNP(1:3).EQ.'his' .OR.
     +    CNP(1:3).EQ.'His') THEN
        IHIGC1 = IVP
        GO TO 120
C
C  HIF - symbol color for high symbols.
C
      ELSE IF (CNP(1:3).EQ.'HIF' .OR. CNP(1:3).EQ.'hif' .OR.
     +    CNP(1:3).EQ.'Hif') THEN
        IHIGC3 = IVP
        GO TO 120
C
C  HIB - character background color for high symbols.
C
      ELSE IF (CNP(1:3).EQ.'HIB' .OR. CNP(1:3).EQ.'hib' .OR.
     +    CNP(1:3).EQ.'Hib') THEN
        IHIGC2 = IVP
        GO TO 120
C
C  HIC - color of circumscribed circle.
C
      ELSE IF (CNP(1:3).EQ.'HIC' .OR. CNP(1:3).EQ.'hic' .OR.
     +    CNP(1:3).EQ.'Hic') THEN
        IHIGC4 = IVP
        GO TO 120
C
C  LOS - shadow color for low symbols.
C
      ELSE IF (CNP(1:3).EQ.'LOS' .OR. CNP(1:3).EQ.'los' .OR.
     +    CNP(1:3).EQ.'Los') THEN
        ILOWC1 = IVP
        GO TO 120
C
C  LOF - symbol color for low symbols.
C
      ELSE IF (CNP(1:3).EQ.'LOF' .OR. CNP(1:3).EQ.'lof' .OR.
     +    CNP(1:3).EQ.'Lof') THEN
        ILOWC3 = IVP
        GO TO 120
C
C  LOB - character background color for low symbols.
C
      ELSE IF (CNP(1:3).EQ.'LOB' .OR. CNP(1:3).EQ.'lob' .OR.
     +    CNP(1:3).EQ.'Lob') THEN
        ILOWC2 = IVP
        GO TO 120
C
C  LOC - color of circumscribed circle.
C
      ELSE IF (CNP(1:3).EQ.'LOC' .OR. CNP(1:3).EQ.'loc' .OR.
     +    CNP(1:3).EQ.'Loc') THEN
        ILOWC4 = IVP
        GO TO 120
C
C  RC1 - box outline color for regional weather labels.
C
      ELSE IF (CNP(1:3).EQ.'RC1' .OR. CNP(1:3).EQ.'rc1' .OR.
     +    CNP(1:3).EQ.'Rc1') THEN
        IRGLC1 = IVP
        GO TO 120
C
C  RC2 - background color for regional weather labels.
C
      ELSE IF (CNP(1:3).EQ.'RC2' .OR. CNP(1:3).EQ.'rc2' .OR.
     +    CNP(1:3).EQ.'Rc2') THEN
        IRGLC2 = IVP
        GO TO 120
C
C  RC3 - shadow color for regional weather labels.
C
      ELSE IF (CNP(1:3).EQ.'RC3' .OR. CNP(1:3).EQ.'rc3' .OR.
     +    CNP(1:3).EQ.'Rc3') THEN
        IRGLC3 = IVP
        GO TO 120
C
C  RC4 - character color for regional weather labels.
C
      ELSE IF (CNP(1:3).EQ.'RC4' .OR. CNP(1:3).EQ.'rc4' .OR.
     +    CNP(1:3).EQ.'Rc4') THEN
        IRGLC4 = IVP
        GO TO 120
C
C  RC5 - character outline color for regional weather labels.
C
      ELSE IF (CNP(1:3).EQ.'RC5' .OR. CNP(1:3).EQ.'rc5' .OR.
     +    CNP(1:3).EQ.'Rc5') THEN
        IRGLC5 = IVP
        GO TO 120
C
C  AWC - arrow color.
C
      ELSE IF (CNP(1:3).EQ.'AWC' .OR. CNP(1:3).EQ.'awc' .OR.
     +    CNP(1:3).EQ.'Awc') THEN
        IAROWC = IVP
        GO TO 120
C
C  SMF - Flags whether interpolation (=0) or smoothing (=1) is done.
C
      ELSE IF (CNP(1:3).EQ.'SMF' .OR. CNP(1:3).EQ.'smf' .OR.
     +    CNP(1:3).EQ.'Smf') THEN
        ISMOTH = IVP
        GO TO 120
C
C  ASC - arrow shadow color.
C
      ELSE IF (CNP(1:3).EQ.'ASC' .OR. CNP(1:3).EQ.'asc' .OR.
     +    CNP(1:3).EQ.'Asc') THEN
        IARSHC = IVP
        GO TO 120
C
C  AOC - arrow outline color.
C
      ELSE IF (CNP(1:3).EQ.'AOC' .OR. CNP(1:3).EQ.'aoc' .OR.
     +    CNP(1:3).EQ.'Aoc') THEN
        IAROUC = IVP
        GO TO 120
C
C  UNT - flags whether to use imperial units (as described in the
C        NOAA chart), or metric units as per:
C          http://weather.unisys.com/wxp/Appendices/Formats/SYNOP.html
C
      ELSE IF (CNP(1:3).EQ.'UNT' .OR. CNP(1:3).EQ.'unt' .OR.
     +    CNP(1:3).EQ.'Unt') THEN
        IUNITS = IVP
        GO TO 120
C
C  VCC - color index for vectors.
C
      ELSE IF (CNP(1:3).EQ.'VCC' .OR. CNP(1:3).EQ.'vcc' .OR.
     +    CNP(1:3).EQ.'Vcc') THEN
        VCCOLR = IVP
        GO TO 120
C
C  VLF - Foreground color for vector label box.
C
      ELSE IF (CNP(1:3).EQ.'VLF' .OR. CNP(1:3).EQ.'vlf' .OR.
     +         CNP(1:3).EQ.'Vlf') THEN
        VLBLFC = IVP
        GO TO 120
C
C  VLB - Background color for vector label box.
C
      ELSE IF (CNP(1:3).EQ.'VLB' .OR. CNP(1:3).EQ.'vlb' .OR.
     +         CNP(1:3).EQ.'Vlb') THEN
        VLBLBC = IVP
        GO TO 120
C
C  VVC - Flags whether to plot the visibility parameter VV as the
C        raw two-character SYNOP codes.
C
      ELSE IF (CNP(1:3).EQ.'VVC' .OR. CNP(1:3).EQ.'vvc' .OR.
     +         CNP(1:3).EQ.'Vvc') THEN
        IVVCOD = IVP
        GO TO 120
      ELSE
        CTM(1:36) = 'WMSETI - Parameter name not known - '
        CTM(37:39) = CNP(1:3)
        CALL SETER (CTM(1:39), 23, 1)
        GO TO 120
      ENDIF
C
  120 CONTINUE
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
C
C $Id: gasetc.f,v 1.6 2008-07-27 00:17:13 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GASETC (PNAM,CVAL)
C
        CHARACTER*(*) PNAM,CVAL
C
C The subroutine GASETC may be used to set GRIDAL parameters which have
C values of type CHARACTER.
C
C Declare the common block containing character parameters.
C
        COMMON /GACHAR/ FNLX,FNLY
        CHARACTER*10    FNLX,FNLY
        SAVE   /GACHAR/
C
C Check for an uncleared prior error.
C
        IF (ICFELL('GASETC - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Set the specified parameter.
C
        IF      (PNAM(1:3).EQ.'XLF'.OR.PNAM(1:3).EQ.'xlf') THEN
          FNLX=CVAL
        ELSE IF (PNAM(1:3).EQ.'YLF'.OR.PNAM(1:3).EQ.'ylf') THEN
          FNLY=CVAL
        ELSE
          CALL SETER ('GASETC - UNRECOGNIZED PARAMETER NAME',2,1)
        END IF
C
C Done.
C
        RETURN
C
      END


!> @brief Collection of subroutines to read in data.
!! @details
!! @author Mark Lomas
!! @date July 2016

module doly

use real_precision
use phenology_methods
use productivity_methods
use misc_values
use site_parameters
use tuning_parameters

implicit none

contains

!**********************************************************************!
!                                                                      !
!                          dolyday :: doly                             !
!                          ---------------                             !
!                                                                      !
!           (D)ynamic gl(O)ba(L) phytogeograph(Y) model                !
! SUBROUTINE DOLYDAY(tmp,prc,hum,ca,soilc,soiln,minn,adp,sfc,sw,sswc,  !
! awl,kd,kx,daygpp,resp_l,rlai,evap,tran,roff,interc,evbs,flow1,flow2, !
! year,mnth,day,pet,ht,thty_dys,ft,lmor_sc,nleaf,leaflitter,hrs,q,     !
! qdirect,qdiff,fpr,canga,gsn,rn,check_closure)                        !
!                                                                      !
!----------------------------------------------------------------------!
!> @brief Read internal parameters from "param.dat" file, and io
!! parameters from "misc_params.dat".
!! @details
!! @author Mark Lomas
!! @date Feb 2006
!----------------------------------------------------------------------!
subroutine dolyday(tmp,prc,hum,ca,soilc,soiln,minn,adp,sfc,sw,sswc, &
 awl,kd,kx,daygpp,resp_l,rlai,evap,tran,roff,interc,evbs,flow1,flow2, &
 year,mnth,day,pet,ht,thty_dys,ft,lmor_sc,nleaf,leaflitter,hrs,q, &
 qdirect,qdiff,fpr,canga,gsn,rn,check_closure)
!**********************************************************************!
real(dp), parameter :: oi = 21000.0
real(dp) :: bb0,bbmax,bblim,sslim,tmp,prc,hum,ca,maxevap, &
 lairat,interc,evbs,t,suma,sum,laiinc, &
 soilc,soiln,rh,tk,rd,rn,pet,petmm,petv,q,qdiff,sla, &
 qdirect,hrs,vpd,canga,lam,rho,s,gam,amx,gsn,dp2,pet2,wtfc, &
 etwt,etmm,ftagh,soilw,maxc,soil2g,can2a,daygpp,daynpp, &
 gsum,sapresp,nppstorx,leafmol,ht,tf(13), &
 eemm,amax,respref,can2g,cangs,svp,rlai, &
 wtwp,p,nppstor2,canres,et,ee,flow1,flow2,tran,rem, &
 adp(4),evap,f2,sfc(4),sw(4),awl(4),kd, &
 kx,f3,minn,roff,windspeed,leaflit,resp_l, &
 npp_eff,eff_dec,eff_age,eff,nppstore,stemnpp, &
 rootnpp,leafnpp,tfscale,sswc(4),yield,yld,resp,fpr, &
 s1in,stemfr,lmor_sc(3600),nleaf,leaflitter,old_total_carbon, &
 total_carbon
integer leafls,stemls,rootls,bbm,ssm,sss,ftphen,c3c4,thty_dys,ft, &
 mnth,i,iter,ndsum(12),lai,day,year,bb,bbgs,ftdth,ss,dsbb, &
 chill,dschill,co
logical veg,check_closure
!----------------------------------------------------------------------!

co = ssp%cohort

!----------------------------------------------------------------------!
! Check carbon closure.
!----------------------------------------------------------------------!
if ((check_closure).and.(pft(co)%sla > 0.0)) then
  old_total_carbon = ssv(co)%lai%tot*12.0/pft(co)%sla/18.0 + &
 ssv(co)%nppstore(1) + ssv(co)%stem%tot + ssv(co)%root%tot + ssv(co)%bio(1) + &
 ssv(co)%bio(2)
endif
!----------------------------------------------------------------------!

wtwp     = ssp%wilt
wtfc     = ssp%field
sla      = pft(ft)%sla
c3c4     = pft(ft)%c3c4
ftphen   = pft(ft)%phen
ftagh    = pft(ft)%crop
ftdth    = pft(ft)%d2h
leafls   = pft(ft)%lls
stemls   = pft(ft)%sls
rootls   = pft(ft)%rls
bbm      = pft(ft)%bbmem
bb0      = pft(ft)%bb0
bbmax    = pft(ft)%bbmax
bblim    = pft(ft)%bblim
ssm      = pft(ft)%senm
sss      = pft(ft)%sens
sslim    = pft(ft)%senlim
lairat   = pft(ft)%lrat

leaflitter = 0.0

!----------------------------------------------------------------------!
! SET parameter for THROUGHFALL.                                       !
!----------------------------------------------------------------------!
tf(1) = 1.0
tf(2) = 0.95
tf(3) = 0.935
tf(4) = 0.92
tf(5) = 0.905
tf(6) = 0.89
tf(7) = 0.875
tf(8) = 0.86
tf(9) = 0.845
tf(10)= 0.83
tf(11)= 0.815
tf(12)= 0.80
tf(13)= 0.785

tfscale = 2.75
do i=1,13
  tf(i) = 1.0 - tfscale*(1.0 - tf(i))
enddo

!----------------------------------------------------------------------!
! Data input for number of days for each month of the year 'ndays',    !
! the mid julian day of each month 'ndsum'.                            !
!----------------------------------------------------------------------!
data ndsum/16,46,75,106,136,167,197,228,259,289,320,350/

p = 101325.0

rlai = ssv(ssp%cohort)%lai%tot
rem = rlai - int(rlai)
lai = int(rlai) + 1

!----------------------------------------------------------------------!
! INITIALISATIONS.                                                     !
!----------------------------------------------------------------------!
! Night-time respiration ? CHECK
rd = 0.82e-6
! maxc, moisture response parameter
maxc = 690.0/622.6
! Leaf molecular weight.
if (sla>0.0) then
  leafmol = 1.0/(sla*18.0)
else
  leafmol =0.0
endif

! initialisation a faire a chaque fois
if (pft(ft)%phen == 0) then
   veg = .false.
else
   veg = .true.
endif

!----------------------------------------------------------------------!
! Water model variables and rain generater.                            !
!----------------------------------------------------------------------!
iter = 30

! unsused anymore... so not important if set to this value every day.
! Better: should be remove from nppcalc
amx = -1000.0
amax = -1000.0

!----------------------------------------------------------------------!
soilw = ssv(ft)%soil_h2o(1) + ssv(ft)%soil_h2o(2) + &
 ssv(ft)%soil_h2o(3) + ssv(ft)%soil_h2o(4)
soil2g = soilw/(ssp%soil_depth*10.0)

if (soil2g>1.0) soil2g = 1.0

!----------------------------------------------------------------------!
! Bare ground IF statement                                             !
!----------------------------------------------------------------------!
if (.not.(veg)) rlai = 0.0
!----------------------------------------------------------------------!

t = tmp
tk = 273.0 + t
rh = hum
!      abmin = tmp(mnth,day)*1.29772 - 19.5362
if (rh>95.0)  rh=95.0
if (rh<30.0)  rh=30.0

! Sapwood respiration N.B. 8.3144 J/K/mol= Universal gas constant.
! sapresp = exp(21.6 - 5.367e4/(8.3144*tk))
sapresp = 0.4*exp(0.06*t)*0.1
sapresp = 0.4*exp(0.02*t)*0.35

! Converted to a monthly figure (N.B. 24 hour day).
! sapresp = (sapresp*3600.0*24.0*30.0)/1000000.0
! Use daily value
sapresp = (sapresp*3600.0*24.0)/1000000.0

! Totalled for the year.
if (soil2g>ssp%wilt) then
  respref = sapresp*tgp%p_resp*((soil2g - ssp%wilt)/&
 min(1.0,(wtfc - ssp%wilt)))**tgp%p_kgw
else
  respref = 0.0
endif

if (t<0.0) then
  respref = 0.0
endif

t = tmp
lam = 2500.0 - (2.367*t)
s = 48.7*exp(0.0532*t)
rn = 0.96*(q*1000000.0/4.0 + 208.0 + 6.0*t)
rn = rn*0.52
pet = (1.26*s*rn)/(s + 66.0)*tgp%p_pet
petmm = (pet*3600.0)/(lam*1000.0)*10.0
if (petmm>0.0) then
  petv = petmm/1.0
else
  petv = 0.0
endif
pet = petv

!----------------------------------------------------------------------!
! canga=k^2 u / (log[(z-d)/z0])^2
! k=von Karman constant. k=0.41
! z=reference height
! d=zero plane displacement
! z0=roughness length
!----------------------------------------------------------------------!
windspeed= 5.0 ! in m/s
canga = 0.168*windspeed/log((200.0 - 0.7*ht)/(0.1*ht))**2

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

npp_eff = 0.0
sum = 0.0
eff_dec = 0.75
eff_age = 90.0
!DO i=1,leafls
!  IF (i.LE.eff_age) THEN
!    eff = 1.0
!  ELSE
!    eff = eff_dec + real(leafls - i)*(1.0 - eff_dec)/(real(leafls) - &
! eff_age + 1.0)
!  ENDIF
!  npp_eff = npp_eff +  leafv(i)*eff
!  sum = sum + leafv(i)
!ENDDO

!IF (abs(sum-rlai).gt.0.001) THEN
!  WRITE(*,*) 'leafv not = rlai ',SUM,RLAI,mnth,day,nppstore
!  STOP
!ENDIF
!IF (sum.GT.0.0) THEN
!  npp_eff = npp_eff/sum
!ELSE
!  npp_eff = 0.0
!ENDIF
npp_eff = 1.0

if (veg) then
  call NPPCALC(npp_eff,c3c4,maxc,soilc,soiln,minn,soil2g,wtwp,wtfc, &
 rd,rlai,t,rh,ca,oi,rn,qdirect,qdiff,can2a,can2g,canres,suma, &
 amx,amax,gsum,hrs,canga/1.3,p,day,nleaf,fpr)

!      write(*,*) 'o',canga,can2a,can2g,canres,suma,gsum

!----------------------------------------------------------------------!
! Calculate the canopy conductance.                                    !
!----------------------------------------------------------------------!
  cangs = 8.3144*tk/p*can2g
  cangs = 1.6*cangs             ! convert CO2 gs into H2O gs
  gsn = cangs
!----------------------------------------------------------------------!
else
  can2a = 0.0
  can2g = 0.0
  canres = 0.
  cangs = 0.0
  gsn = 0.0
  suma = 0.0
endif

call SUMA_ADD(suma)

dp2 = prc
pet2 = petv

pet = eemm
pet2 = pet

!----------------------------------------------------------------------!
! Set switch 'stseas' to find when a new growing season has started.   *
! 'stseas = 1' for the first day of a new growing season.              !
!----------------------------------------------------------------------!

soil2g = soilw/(ssp%soil_depth*10.0)

!----------------------------------------------------------------------!
!                         DAILY npp                                    !
!----------------------------------------------------------------------!
daynpp = can2a*3600.0*hrs/1000000.0 - canres*3600.0*(24.0-hrs)/1000000.0
daynpp = daynpp*12.0
resp_l = canres*3600.0*(24.0-hrs)*12.0/1000000.0
daygpp= can2a*3600.0*hrs*12.0/1000000.0

!----------------------------------------------------------------------!
! End of the DAILY LOOP.                                               !
!----------------------------------------------------------------------!

!----------------------------------------------------------------------!
! Check carbon closure.
!----------------------------------------------------------------------!
if ((check_closure).and.(pft(co)%sla > 0.0)) then
  total_carbon = ssv(co)%lai%tot*12.0/pft(co)%sla/18.0 + &
 ssv(co)%nppstore(1) + ssv(co)%stem%tot + ssv(co)%root%tot + ssv(co)%bio(1) + &
 ssv(co)%bio(2)
  if (abs(total_carbon-old_total_carbon) >  1.0e-3) then
    write(*,*) 'Breach of carbon closure in DOLYDAY:', &
 total_carbon-old_total_carbon,' g/m^2.'
    stop
  endif
endif

end subroutine dolyday

!**********************************************************************!
!                                                                      !
!                     evapotranspiration :: doly                       !
!                     --------------------------                       !
!                                                                      !
! SUBROUTINE evapotranspiration(t,rh,rn,canga,gsn,hrs,eemm,etmm)       !
!                                                                      !
!----------------------------------------------------------------------!
!> @brief Read internal parameters from "param.dat" file, and io
!! parameters from "misc_params.dat".
!! @details
!! @author Mark Lomas
!! @date Feb 2006
!----------------------------------------------------------------------!
subroutine evapotranspiration(t,rh,rn,canga,gsn,hrs,eemm,etmm)
!**********************************************************************!
real(dp) :: t,rh,rn,canga,gsn,eemm,etmm,et,svp,vpd,lam,rho,s,gam,hrs, &
 maxevap,etwt,ee

!----------------------------------------------------------------------!
! Penman-Monteith equation for evapotranspiration.                     !
! Units of ET = W/m2  (CHECK) N.B. Conductances in moles.              !
!----------------------------------------------------------------------!
svp = 6.108*exp((17.269*t)/(237.3 + t))*100.0
vpd = (1.0 - rh/100.0)*svp
lam = 2500.0 - 2.367*t
rho = 1288.4 - 4.103*t
s = 48.7*exp(0.0532*t)
gam = 101325.0*1.012/(0.622*lam)

if ((ssv(ssp%cohort)%lai%tot>0.1).and.(msv%mv_soil2g>ssp%wilt)) then
  et = (s*rn + rho*1.012*canga*vpd)/(s + gam*(1.0 + canga/gsn))*tgp%p_et
! watch dog ajoute par ghislain
  if (et<0.0) then
    et=0.0
  endif

  etwt = (et*3600.0*hrs)/lam
  etmm = etwt/1000.0
  if (etmm>0.0) then
    if (etmm/hrs>maxevap)  maxevap = etmm/hrs
  endif
else
  et = 0.0
  etwt = 0.0
  etmm = 0.0
endif

ee = (s*rn + rho*1.012*canga*vpd)/(s + gam)
eemm = (ee*3600.0*hrs)/(lam*1000.0)

! added by Ghislain 20/10/03
if (ee<0.0) then
! WRITE(*,*) 'ee is negativ ee=',ee,rn
  ee=0.0
  eemm=0.0
endif

end subroutine evapotranspiration

end module doly


text:
><Image(Surprised.JPG)

><James was premed at UC Davis, and graduated with a B.S. in College of Biological Sciences Biological Sciences in Spring 2006.

><James has red hair and chicken legs.  He is surprisingly good at Softball softball despite these features.  He plays for Heaven Help Us and likes to eat ice cream at Baskin Robbins after the games.

><He worked at BJs Restaurant and Brewery BJs (until very recently) and thinks their pizza is awesome. He is too lazy to make his own page, which is why he gave Users/JanelleAlvstadMattson Janelle full authorization to make this page for him.

20070709 22:24:55 nbsp Hi James, welcome to the wiki and thanks for your contributions. Please remember to preview your edits. Thanks! Users/ArlenAbraham


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Log: graexp.f,v $
C Revision 1.2  2009/03/25 12:36:44  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.1  1994/04/07 20:01:58  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.2  1990/12/14  08:51:19  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE GRAEXP (AXTYP, NNUM, RNUM, ATIC)
C=======================================================================

C   --*** GRAEXP *** (GRPLIB) Choose axis exponent (PLT)
C   --   Written by Amy Gilkey - revised 02/20/87
C   --
C   --GRAEXP chooses an exponent from the axes limits and sets this
C   --value and the number of decimal digits for the given axis.
C   --It will convert to engineering notation if possible.
C   --
C   --Parameters:
C   --   AXTYP - IN - the axis (' ' for both, 'X' for X, 'Y' for Y)
C   --   NNUM - IN - the number of numbers in the array set;
C   --      0 for no exponent
C   --   RNUM - IN - the array of axes limits
C   --   ATIC - IN - the axis tick mark interval

C   --Routines Called:
C   --   PLTSTG - (PLTLIB) Set graph parameter
C   --      13, 14 = (KXEXP, KYEXP) X, Y axis exponent
C   --      19, 20 = (KXNDIG, KYNDIG) X, Y axis number of decimal digits

      PARAMETER (KXEXP=13, KYEXP=14, KXNDIG=19, KYNDIG=20)

      CHARACTER AXTYP
      INTEGER NNUM
      REAL RNUM(*)
      REAL ATIC

      LOGICAL LDUM, PLTSTG
      CHARACTER*12 SCRSTR

      IF (NNUM .GT. 0) THEN

C      --Convert all to E notation and find the minimum and maximum exponent
C      --   MINE and MAXE are the minimum and maximum exponents

         MINE = 999
         MAXE = -999
         DO 100 I = 1, NNUM
            IF (RNUM(I) .NE. 0.0) THEN
               WRITE (SCRSTR(1:8), '(0PE8.1)') RNUM(I)
               READ (SCRSTR(6:8), '(I3)') IE
               MINE = MIN (MINE, IE)
               MAXE = MAX (MAXE, IE)
            END IF
  100    CONTINUE

         IF (MAXE .GT. 0) THEN
            IEXP = INT ((MAXE - 1) / 3) * 3
         ELSE
            IEXP = INT ((MAXE - 2) / 3) * 3
         END IF

         WRITE (SCRSTR(1:12), '(0PE12.5)') ATIC
         READ (SCRSTR(10:12), '(I3)') ITICE
         DO 110 I = 8, 3, -1
            IF (SCRSTR(I:I) .NE. '0') GOTO 120
  110    CONTINUE
  120    CONTINUE
         NSIG = MAX (2, I-3 + MAXE - ITICE)
         NDIG = MAX (0, IEXP - MAXE + NSIG)

      ELSE

C      --Set exponent and number of digits to zero

         IEXP = 0
         NDIG = 0
      END IF

      IF (AXTYP .NE. 'Y') THEN
         LDUM = PLTSTG (KXEXP, REAL(IEXP))
         LDUM = PLTSTG (KXNDIG, REAL(NDIG))
      END IF
      IF (AXTYP .NE. 'X') THEN
         LDUM = PLTSTG (KYEXP, REAL(IEXP))
         LDUM = PLTSTG (KYNDIG, REAL(NDIG))
      END IF

      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
function increment_value(a) result(b)
  integer :: a
  integer :: b

  b = a + 1
end function

program main
   integer :: d = 2

   print *, increment_value(d)
end program

max_stars_repo_name: ouankou/rose

max_issues_repo_name: ouankou/rose

max_forks_repo_name: ouankou/rose