text:
      subroutine gg_ZZ(p,msqgg)
      implicit none
c--- Author: J. M. Campbell, September 2013
c--- Matrix element squared for the process gg->ZZ
c--- Effects of massive quarks in the third generation may be included
c--- (default: included)
      include 'constants.f'
      include 'ewcouple.f'
      include 'qcdcouple.f'
      include 'docheck.f'
      include 'interference.f'
      logical includegens1and2,includebottom,includetop
      integer h1,h2,h34,h56
      double precision p(mxpart,4),pswap(mxpart,4),msqgg,fac,oprat
      double complex
     & Mloop_uptype(2,2,2,2),Mloop_dntype(2,2,2,2),
     & Mloop_bquark(2,2,2,2),Mloop_tquark(2,2,2,2),
     & Sloop_uptype(2,2,2,2),Sloop_dntype(2,2,2,2),
     & Sloop_bquark(2,2,2,2),Sloop_tquark(2,2,2,2),
     & Mamp,Samp

c--- set this to true to include generations 1 and 2 of (light) quarks
      includegens1and2=.true.
c--- set this to true to include massive bottom quark
      includebottom=.true.
c--- set this to true to include massive top quark
      includetop=.true.

c--- if set, performs check against numerical results at specific PS point
      docheck=.false.

c--- compute all gg->ZZ amplitudes
      call getggZZamps(p,includegens1and2,includebottom,includetop,
     & Mloop_uptype,Mloop_dntype,Mloop_bquark,Mloop_tquark)

      if (interference) then
c--- for interference, compute amplitudes after 4<->6 swap
       pswap(1,:)=p(1,:)
       pswap(2,:)=p(2,:)
       pswap(3,:)=p(3,:)
       pswap(4,:)=p(6,:)
       pswap(5,:)=p(5,:)
       pswap(6,:)=p(4,:)
       call getggZZamps(pswap,includegens1and2,includebottom,includetop,
     &  Sloop_uptype,Sloop_dntype,Sloop_bquark,Sloop_tquark)
      endif

      msqgg=0d0
      do h1=1,2
      do h2=1,2
      do h34=1,2
      do h56=1,2

      Mamp=
     &  +2d0*Mloop_uptype(h1,h2,h34,h56)
     &  +2d0*Mloop_dntype(h1,h2,h34,h56)
     &      +Mloop_bquark(h1,h2,h34,h56)
     &      +Mloop_tquark(h1,h2,h34,h56)

      if (interference .eqv. .false.) then
c--- normal case
        msqgg=msqgg+cdabs(Mamp)**2
      else
c--- with interference
        Samp=
     &  +2d0*Sloop_uptype(h1,h2,h34,h56)
     &  +2d0*Sloop_dntype(h1,h2,h34,h56)
     &      +Sloop_bquark(h1,h2,h34,h56)
     &      +Sloop_tquark(h1,h2,h34,h56)
        if (h34 .eq. h56) then
          oprat=1d0-2d0*dble(dconjg(Mamp)*Samp)
     &                 /(cdabs(Mamp)**2+cdabs(Samp)**2)
        else
          oprat=1d0
        endif

        msqgg=msqgg+cdabs(Mamp)**2*oprat
     &  +cdabs(Samp)**2*oprat

      endif

      enddo
      enddo
      enddo
      enddo

c--- overall factor extracted (c.f. getggZZamps.f)
      fac=avegg*V*(4d0*esq*gsq/(16d0*pisq)*esq)**2

      msqgg=msqgg*fac*vsymfact

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
module marquardt
use vars
use SEE
use rt_coef
use synth
use svd
use maths
use io
use l_bfgs_b, only : setulb
implicit none
contains

!*************************************************************
!*************************************************************
! MARQUARDT METHOD
!*************************************************************
!*************************************************************

!------------------------------------------------------------
! Calculate the chi^2 function
!------------------------------------------------------------
	function compute_chisq(in_observation,in_inversion)
	type(type_inversion) :: in_inversion
	type(type_observation) :: in_observation
	real(kind=8) :: compute_chisq, weight
	integer :: nf, i, j

		compute_chisq = 0.d0

! Compute the chisq using the appropriate weights for each cycle
		do i = 0, 3
			weight = in_inversion%stokes_weights(i,in_inversion%loop_cycle)
			compute_chisq = compute_chisq + weight * &
				sum((in_observation%stokes(i,:)-in_inversion%stokes_unperturbed(i,:))**2/in_observation%sigma(i,:)**2) / &
				(4.d0*in_observation%n)
		enddo

	end function compute_chisq

!------------------------------------------------------------
! Calculate the maximum relative change in the parameters
!------------------------------------------------------------
	function compute_params_relative_change(params,trial)
	real(kind=8) :: compute_params_relative_change
	type(variable_parameters) :: params, trial
	real(kind=8), allocatable :: rel_change(:)
	integer :: i, j

		allocate(rel_change(params%n_total))

		if (params%bgauss /= 0.d0) then
			rel_change(1) = abs((params%bgauss-trial%bgauss) / params%bgauss)
		else
			rel_change(1) = 0.d0
		endif

		if (params%thetabd /= 0.d0) then
			rel_change(2) = abs((params%thetabd-trial%thetabd) / params%thetabd)
		else
			rel_change(2) = 0.d0
		endif

		if (params%chibd /= 0.d0) then
			rel_change(3) = abs((params%chibd-trial%chibd) / params%chibd)
		else
			rel_change(3) = 0.d0
		endif

		if (params%vdopp /= 0.d0) then
			rel_change(4) = abs((params%vdopp-trial%vdopp) / params%vdopp)
		else
			rel_change(4) = 0.d0
		endif

		if (params%dtau /= 0.d0) then
			rel_change(5) = abs((params%dtau-trial%dtau) / params%dtau)
		else
			rel_change(5) = 0.d0
		endif

		if (params%delta_collision /= 0.d0) then
			rel_change(6) = abs((params%delta_collision-trial%delta_collision) / params%delta_collision)
		else
			rel_change(6) = 0.d0
		endif

		if (params%vmacro /= 0.d0) then
			rel_change(7) = abs((params%vmacro-trial%vmacro) / params%vmacro)
		else
			rel_change(7) = 0.d0
		endif

		if (params%damping /= 0.d0) then
			rel_change(8) = abs((params%damping-trial%damping) / params%damping)
		else
			rel_change(8) = 0.d0
		endif

		if (params%beta /= 0.d0) then
			rel_change(9) = abs((params%beta-trial%beta) / params%beta)
		else
			rel_change(9) = 0.d0
		endif

		if (params%height /= 0.d0) then
			rel_change(10) = abs((params%height-trial%height) / params%height)
		else
			rel_change(10) = 0.d0
		endif

		if (params%dtau2 /= 0.d0) then
			rel_change(11) = abs((params%dtau2-trial%dtau2) / params%dtau2)
		else
			rel_change(11) = 0.d0
		endif

		if (params%vmacro2 /= 0.d0) then
			rel_change(12) = abs((params%vmacro2-trial%vmacro2) / params%vmacro2)
		else
			rel_change(12) = 0.d0
		endif

		if (params%bgauss2 /= 0.d0) then
			rel_change(13) = abs((params%bgauss2-trial%bgauss2) / params%bgauss2)
		else
			rel_change(13) = 0.d0
		endif

		if (params%thetabd2 /= 0.d0) then
			rel_change(14) = abs((params%thetabd2-trial%thetabd2) / params%thetabd2)
		else
			rel_change(14) = 0.d0
		endif

		if (params%chibd2 /= 0.d0) then
			rel_change(15) = abs((params%chibd2-trial%chibd2) / params%chibd2)
		else
			rel_change(15) = 0.d0
		endif

		if (params%vdopp2 /= 0.d0) then
			rel_change(16) = abs((params%vdopp2-trial%vdopp2) / params%vdopp2)
		else
			rel_change(16) = 0.d0
		endif

		if (params%ff /= 0.d0) then
			rel_change(17) = abs((params%ff-trial%ff) / params%ff)
		else
			rel_change(17) = 0.d0
		endif

		compute_params_relative_change = maxval(rel_change)*100.d0

		deallocate(rel_change)

	end function compute_params_relative_change

!------------------------------------------------------------
! Check if the values are inside physical boundaries
! If not inside the boundaries, use the previous value
!------------------------------------------------------------
	subroutine check_boundaries(original,trial,correct)
	type(variable_parameters) :: original, trial
	logical :: correct
	real(kind=8) :: upper(17), lower(17)
	integer :: i

! Read ranges
		open(unit=24,file=direct_ranges,action='read',status='old')
		call lb(24,12)
		do i = 1, 17
			call lb(24,2)
			read(24,*) lower(i), upper(i)
		enddo
		close(24)

! Verify if the parameters are inside the boundaries. If any of them is outside, use the previous value
! of this parameter.

		if (trial%bgauss < lower(1) .or. trial%bgauss > upper(1)) then
			trial%bgauss = original%bgauss
		endif
		if (trial%thetabd < lower(2) .or. trial%thetabd > upper(2)) then
			trial%thetabd = original%thetabd
		endif
		if (trial%chibd < lower(3) .or. trial%chibd > upper(3)) then
			trial%chibd = original%chibd
		endif
		if (trial%vdopp < lower(4) .or. trial%vdopp > upper(4)) then
			trial%vdopp = original%vdopp
		endif
		if (trial%dtau < lower(5) .or. trial%dtau > upper(5)) then
			trial%dtau = original%dtau
		endif
		if (original%nslabs /= 1) then
			if (trial%dtau2 < lower(11) .or. trial%dtau2 > upper(11)) then
				trial%dtau2 = original%dtau2
			endif
		endif
		if (trial%delta_collision < lower(6) .or. trial%delta_collision > upper(6)) then
			trial%delta_collision = original%delta_collision
		endif
		if (trial%vmacro < lower(7) .or. trial%vmacro > upper(7)) then
			trial%vmacro = original%vmacro
		endif
		if (original%nslabs /= 1) then
			if (trial%vmacro2 < lower(12) .or. trial%vmacro2 > upper(12)) then
				trial%vmacro2 = original%vmacro2
			endif
		endif
		if (trial%damping < lower(8) .or. trial%damping > upper(8)) then
			trial%damping = original%damping
		endif
		if (trial%beta < lower(9) .or. trial%beta > upper(9)) then
			trial%beta = original%beta
		endif
		if (trial%height < lower(10) .or. trial%height > upper(10)) then
			trial%height = original%height
		endif

! If two components with two different fields
		if (original%nslabs == 3 .or. original%nslabs == -2) then
			if (trial%bgauss2 < lower(13) .or. trial%bgauss2 > upper(13)) then
				trial%bgauss2 = original%bgauss2
			endif
			if (trial%thetabd2 < lower(14) .or. trial%thetabd2 > upper(14)) then
				trial%thetabd2 = original%thetabd2
			endif
			if (trial%chibd2 < lower(15) .or. trial%chibd2 > upper(15)) then
				trial%chibd2 = original%chibd2
			endif
			if (trial%vdopp2 < lower(16) .or. trial%vdopp2 > upper(16)) then
				trial%vdopp2 = original%vdopp2
			endif
		endif

! Two slabs in the same pixel
		if (original%nslabs == -2) then
			if (trial%ff < lower(17) .or. trial%ff > upper(17)) then
				trial%ff = original%ff
			endif
		endif

		correct = .TRUE.

	end subroutine check_boundaries

!------------------------------------------------------------
! Perturb a given parameter
!------------------------------------------------------------
	subroutine perturb_parameter(in_params,out_params,which,perturbation,delta)
	type(variable_parameters) :: in_params, out_params
	integer :: which
	real(kind=8) :: perturbation, delta

		select case(which)
			case(1)
				if (in_params%bgauss > 1.d-2) then
					out_params%bgauss = in_params%bgauss * (1.d0+perturbation)
				else
					out_params%bgauss = in_params%bgauss + perturbation
				endif
				delta = out_params%bgauss - in_params%bgauss

			case(2)
				if (abs(in_params%thetabd) > 1.d-2) then
					out_params%thetabd = in_params%thetabd * (1.d0+perturbation)
				else
					out_params%thetabd = in_params%thetabd + perturbation
				endif
				delta = out_params%thetabd - in_params%thetabd

			case(3)
				if (abs(in_params%chibd) > 1.d-2) then
					out_params%chibd = in_params%chibd * (1.d0+perturbation)
				else
					out_params%chibd = in_params%chibd + perturbation
				endif
				delta = out_params%chibd - in_params%chibd

			case(4)
				if (in_params%vdopp > 1.d-2) then
					out_params%vdopp = in_params%vdopp * (1.d0+perturbation)
				else
					out_params%vdopp = in_params%vdopp + perturbation
				endif
				delta = out_params%vdopp - in_params%vdopp

			case(5)
				if (in_params%dtau > 1.d-2) then
					out_params%dtau = in_params%dtau * (1.d0+perturbation)
				else
					out_params%dtau = in_params%dtau + perturbation
				endif
				delta = out_params%dtau - in_params%dtau

			case(6)
				if (in_params%delta_collision > 1.d-2) then
					out_params%delta_collision = in_params%delta_collision * (1.d0+perturbation)
				else
					out_params%delta_collision = in_params%delta_collision + perturbation
				endif
				delta = out_params%delta_collision - in_params%delta_collision

			case(7)
				if (abs(in_params%vmacro) > 1.d-2) then
					out_params%vmacro = in_params%vmacro * (1.d0+perturbation)
				else
					out_params%vmacro = in_params%vmacro + perturbation
				endif
				delta = out_params%vmacro - in_params%vmacro

			case(8)
				if (in_params%damping > 1.d-2) then
					out_params%damping = in_params%damping * (1.d0+perturbation)
				else
					out_params%damping = in_params%damping + perturbation
				endif
				delta = out_params%damping - in_params%damping

			case(9)
				if (abs(in_params%beta) > 1.d-2) then
					out_params%beta = in_params%beta * (1.d0+perturbation)
				else
					out_params%beta = in_params%beta + perturbation
				endif
				delta = out_params%beta - in_params%beta

			case(10)
				if (in_params%height > 1.d-2) then
					out_params%height = in_params%height * (1.d0+perturbation)
				else
					out_params%height = in_params%height + perturbation
				endif
				delta = out_params%height - in_params%height

			case(11)
				if (in_params%dtau2 > 1.d-2) then
					out_params%dtau2 = in_params%dtau2 * (1.d0+perturbation)
				else
					out_params%dtau2 = in_params%dtau2 + perturbation
				endif
				delta = out_params%dtau2 - in_params%dtau2

			case(12)
				if (abs(in_params%vmacro2) > 1.d-2) then
					out_params%vmacro2 = in_params%vmacro2 * (1.d0+perturbation)
				else
					out_params%vmacro2 = in_params%vmacro2 + perturbation
				endif
				delta = out_params%vmacro2 - in_params%vmacro2

			case(13)
				if (in_params%bgauss2 > 1.d-2) then
					out_params%bgauss2 = in_params%bgauss2 * (1.d0+perturbation)
				else
					out_params%bgauss2 = in_params%bgauss2 + perturbation
				endif
				delta = out_params%bgauss2 - in_params%bgauss2

			case(14)
				if (abs(in_params%thetabd2) > 1.d-2) then
					out_params%thetabd2 = in_params%thetabd2 * (1.d0+perturbation)
				else
					out_params%thetabd2 = in_params%thetabd2 + perturbation
				endif
				delta = out_params%thetabd2 - in_params%thetabd2

			case(15)
				if (abs(in_params%chibd2) > 1.d-2) then
					out_params%chibd2 = in_params%chibd2 * (1.d0+perturbation)
				else
					out_params%chibd2 = in_params%chibd2 + perturbation
				endif
				delta = out_params%chibd2 - in_params%chibd2

			case(16)
				if (in_params%vdopp2 > 1.d-2) then
					out_params%vdopp2 = in_params%vdopp2 * (1.d0+perturbation)
				else
					out_params%vdopp2 = in_params%vdopp2 + perturbation
				endif
				delta = out_params%vdopp2 - in_params%vdopp2

			case(17)
				if (in_params%ff > 1.d-2) then
					out_params%ff = in_params%ff * (1.d0+perturbation)
				else
					out_params%ff = in_params%ff + perturbation
				endif
				delta = out_params%ff - in_params%ff

		end select

	end subroutine perturb_parameter

!------------------------------------------------------------
! Perturb a given parameter
!------------------------------------------------------------
	subroutine add_to_parameter(in_params,out_params,which,new_value)
	type(variable_parameters) :: in_params, out_params
	integer :: which
	real(kind=8) :: new_value, rel_change

! Limit the maximum relative change to 50%
		if (new_value > 0.5d0) then
			write(*,FMT='(A,A)') 'Clipping ', trim(adjustl(parameters_name(which)))
			new_value = 0.5d0
		endif
		select case(which)
			case(1)
				out_params%bgauss = in_params%bgauss + new_value
			case(2)
				out_params%thetabd = in_params%thetabd + new_value
			case(3)
				out_params%chibd = in_params%chibd + new_value
			case(4)
				out_params%vdopp = in_params%vdopp + new_value
			case(5)
				out_params%dtau = in_params%dtau + new_value
			case(6)
				out_params%delta_collision = in_params%delta_collision + new_value
			case(7)
				out_params%vmacro = in_params%vmacro + new_value
			case(8)
				out_params%damping = in_params%damping + new_value
			case(9)
				out_params%beta = in_params%beta + new_value
			case(10)
				out_params%height = in_params%height + new_value
			case(11)
				out_params%dtau2 = in_params%dtau2 + new_value
			case(12)
				out_params%vmacro2 = in_params%vmacro2 + new_value
			case(13)
				out_params%bgauss2 = in_params%bgauss2 + new_value
			case(14)
				out_params%thetabd2 = in_params%thetabd2 + new_value
			case(15)
				out_params%chibd2 = in_params%chibd2 + new_value
			case(16)
				out_params%vdopp2 = in_params%vdopp2 + new_value
			case(17)
				out_params%ff = in_params%ff + new_value
		end select

	end subroutine add_to_parameter

!------------------------------------------------------------
! Calculate the derivatives of the Stokes profiles
!------------------------------------------------------------
	subroutine compute_dydx(in_params,in_fixed,in_inversion,in_observation)
	type(variable_parameters) :: in_params, params_scaled, params_scaled_modified, params_modified
	type(fixed_parameters) :: in_fixed
	type(type_inversion) :: in_inversion
	type(type_observation) :: in_observation
	integer :: i
	real(kind=8) :: delta

! Do the synthesis without perturbing the magnetic field
! 		call do_synthesis(in_params, in_fixed, in_observation, in_inversion%stokes_unperturbed)

		params_scaled = in_params
		params_modified = in_params

		call compress_parameters(in_params, params_scaled)
		params_scaled_modified = params_scaled

		do i = 1, in_params%n_total

			if (in_params%inverted(i) == 1) then
				write(*,FMT='(A,I2,A,A)') 'Perturbing parameter ', i, '  --- ', trim(adjustl(parameters_name(i)))

! Perturb the given parameter and do another synthesis
				call perturb_parameter(params_scaled,params_scaled_modified,i,0.0001d0,delta)

				call expand_parameters(params_scaled_modified, params_modified)

				call do_synthesis(params_modified, in_fixed, in_observation, in_inversion%stokes_perturbed)

				in_inversion%dydx(:,:,i) = (in_inversion%stokes_perturbed - in_inversion%stokes_unperturbed) / delta
			endif
		enddo

	end subroutine compute_dydx

!------------------------------------------------------------
! Calculate the derivatives of the Stokes profiles
!------------------------------------------------------------
	subroutine compute_trial_params(in_params,in_fixed,in_inversion,in_observation,in_trial)
	type(variable_parameters) :: in_params, in_trial, in_scaled, in_temp
	type(fixed_parameters) :: in_fixed
	type(type_inversion) :: in_inversion
	type(type_observation) :: in_observation
	real(kind=8), allocatable :: alpha(:,:), beta(:), w(:), v(:,:), x(:)
	real(kind=8) :: obs, syn, sig, weight, wmax
	integer :: i, j, k, l, np

		in_trial = in_params
		call compress_parameters(in_params,in_scaled)
		in_temp = in_scaled

		np = in_params%n_total

		allocate(alpha(np,np))
		allocate(beta(np))
		allocate(v(np,np))
		allocate(w(np))
		allocate(x(np))

		alpha = 0.d0
		beta = 0.d0

		do i = 1, in_fixed%no
			do j = 0, 3
				weight = in_inversion%stokes_weights(j,in_inversion%loop_cycle)
				do k = 1, np

					obs = in_observation%stokes(j,i)
					syn = in_inversion%stokes_unperturbed(j,i)
					sig = in_observation%sigma(j,i)

					if (in_params%inverted(k) == 1) then
						beta(k) = beta(k) - 2.d0 * weight * (obs-syn) / sig**2 * in_inversion%dydx(j,i,k) / (4.d0*in_fixed%no)
						do l = 1, np
							if (in_params%inverted(l) == 1) then
								alpha(k,l) = alpha(k,l) + 2.d0 * weight * &
									in_inversion%dydx(j,i,k) * in_inversion%dydx(j,i,l) / sig**2 / (4.d0*in_fixed%no)
							endif
						enddo
					endif
				enddo
			enddo
		enddo

		beta = -0.5d0 * beta
		alpha = 0.5d0 * alpha

! Apply the Lambda parameter of the Levenverb-Marquardt
		do i = 1, np
			alpha(i,i) = (1.d0+inversion%lambda) * alpha(i,i)
		enddo

		np = in_params%n_total
		call svdcmp(alpha,np,np,np,np,w,v)

! Regularize the linear system by setting all singular values below a certain
! threshold equal to zero
		wmax = maxval(w)
		do i = 1, np
			if (w(i) < wmax * 1.d-6) then
				w(i) = 0.d0
			endif
		enddo

		call svbksb(alpha,w,v,np,np,np,np,beta,x)

		do i = 1, np
			if (in_params%inverted(i) == 1) then
				call add_to_parameter(in_scaled,in_temp,i,x(i)/1.d0)
			endif
		enddo

		call expand_parameters(in_temp,in_trial)

		deallocate(alpha)
		deallocate(beta)
		deallocate(v)
		deallocate(w)
		deallocate(x)

	end subroutine compute_trial_params

!------------------------------------------------------------
! Re-scale the parameters to make them of the order of unity
!------------------------------------------------------------
	subroutine compress_parameters(in_params,in_params_scaled)
	type(variable_parameters) :: in_params, in_params_scaled
		in_params_scaled%bgauss = in_params%bgauss / 1000.d0
		in_params_scaled%thetabd = in_params%thetabd / 10.d0
		in_params_scaled%chibd = in_params%chibd / 10.d0
		in_params_scaled%vdopp = in_params%vdopp / 5.d0
		in_params_scaled%dtau = in_params%dtau / 1.5d0
		in_params_scaled%delta_collision = in_params%delta_collision
		in_params_scaled%vmacro = in_params%vmacro / 5.d0
		in_params_scaled%damping = in_params%damping
		in_params_scaled%beta = in_params%beta / 3.d0
		in_params_scaled%height = in_params%height / 30.d0

! Two components with same field and Doppler width
		if (in_params%nslabs /= 1) then
			in_params_scaled%dtau2 = in_params%dtau2 / 1.5d0
			in_params_scaled%vmacro2 = in_params%vmacro2 / 5.d0
		endif

! Two components with different field and Doppler width
		if (in_params%nslabs == 3 .or. in_params%nslabs == -2) then
			in_params_scaled%bgauss2 = in_params%bgauss2 / 1000.d0
			in_params_scaled%thetabd2 = in_params%thetabd2 / 10.d0
			in_params_scaled%chibd2 = in_params%chibd2 / 10.d0
			in_params_scaled%vdopp2 = in_params%vdopp2 / 5.d0
		endif

! Two components in the same pixel
		if (in_params%nslabs == -2) then
			in_params_scaled%ff = in_params%ff
		endif

	end subroutine compress_parameters

!------------------------------------------------------------
! Re-scale the parameters to their original values
!------------------------------------------------------------
	subroutine expand_parameters(in_params_scaled,in_params)
	type(variable_parameters) :: in_params, in_params_scaled
		in_params%bgauss = in_params_scaled%bgauss * 1000.d0
		in_params%thetabd = in_params_scaled%thetabd * 10.d0
		in_params%chibd = in_params_scaled%chibd * 10.d0
		in_params%vdopp = in_params_scaled%vdopp * 5.d0
		in_params%dtau = in_params_scaled%dtau * 1.5d0
		in_params%delta_collision = in_params_scaled%delta_collision
		in_params%vmacro = in_params_scaled%vmacro * 5.d0
		in_params%damping = in_params_scaled%damping
		in_params%beta = in_params_scaled%beta * 3.d0
		in_params%height = in_params_scaled%height * 30.d0

! Two components with same field and Doppler width
		if (in_params%nslabs /= 1) then
			in_params%dtau2 = in_params_scaled%dtau2 * 1.5d0
			in_params%vmacro2 = in_params_scaled%vmacro2 * 5.d0
		endif

! Two components with different field and Doppler width
		if (in_params%nslabs == 3 .or. in_params%nslabs == -2) then
			in_params%bgauss2 = in_params_scaled%bgauss2 * 1000.d0
			in_params%thetabd2 = in_params_scaled%thetabd2 * 10.d0
			in_params%chibd2 = in_params_scaled%chibd2 * 10.d0
			in_params%vdopp2 = in_params_scaled%vdopp2 * 5.d0
		endif

! Two components in the same pixel
		if (in_params%nslabs == -2) then
			in_params%ff = in_params_scaled%ff
		endif

	end subroutine expand_parameters

!*************************************************************
!*************************************************************
! DIRECT METHOD
!*************************************************************
!*************************************************************

!------------------------------------------------------------
! Invert some parameters with the DIRECT method
!------------------------------------------------------------
	subroutine invert_with_direct(out_params, synth_option)
	real(kind=8) :: DIReps, DIRf
	integer :: ndim, DIRmaxf, DIRmaxT, synth_option
	integer :: DIRalg
	integer :: IError, logfile
	real(kind=8) :: fglobal, fglper, volper, sigmaper
	real(kind=8), allocatable :: u(:), l(:), upper(:), lower(:)
	integer :: n, i, j
	real(kind=8), allocatable :: DIRx(:)

	integer, parameter :: iisize=300, idsize=300, icsize=30
	integer :: iidata(iisize)
	real(kind=8) :: ddata(idsize)
	character(len=40) :: cdata(icsize)
	character(len=120) :: location_file_direct

	type(variable_parameters) :: out_params

		open(unit=23,file='logfile',action='write',status='replace')

		logfile = 23
		fglobal = -1.d100
		fglper = 0.d0
		sigmaper  = -1.d0

		DIReps = 1.d-6
		DIRmaxT = 6000

		ndim = sum(params%inverted)
		allocate(u(ndim))
		allocate(l(ndim))
		allocate(DIRx(ndim))

		allocate(upper(17))
		allocate(lower(17))
		open(unit=24,file=direct_ranges,action='read',status='old')
		call lb(24,3)

		call lb(24,2)
		read(24,*) location_file_direct

		open(unit=25,file=location_file_direct,action='write',status='replace')

		call lb(24,2)
		read(24,*) DIRmaxf
		if (DIRmaxf < 0) DIRmaxf = 10000
		print *, 'Maximum number of function evaluations : ', DIRmaxf

		call lb(24,2)
		read(24,*) volper
		if (volper < 0) volper = 1.d-20
		print *, 'Stopping when the hypervolume with respect to the original is less than ', volper

		do i = 1, 17
			call lb(24,2)
			read(24,*) lower(i), upper(i)
		enddo
		close(24)

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				l(j) = lower(i)
				u(j) = upper(i)
				print *, 'Inverting ', trim(adjustl(parameters_name(i)))
				j = j + 1
			endif
		enddo

		deallocate(upper)
		deallocate(lower)

		DIRalg = 1

		if (synth_option == 2) then
			call DIRECT(fcn,DIRx,ndim,DIReps,DIRmaxf,DIRmaxT, DIRf, l, u, DIRalg, Ierror, logfile, &
				fglobal, fglper, volper, sigmaper, iidata, iisize, ddata, idsize, cdata, icsize)
		else
			call DIRECT(fcn_simplified_StokesI,DIRx,ndim,DIReps,DIRmaxf,DIRmaxT, DIRf, l, u, DIRalg, Ierror, logfile, &
				fglobal, fglper, volper, sigmaper, iidata, iisize, ddata, idsize, cdata, icsize)
		endif

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				select case(i)
					case(1)
						out_params%bgauss = DIRx(j)
					case(2)
						out_params%thetabd = DIRx(j)
					case(3)
						out_params%chibd = DIRx(j)
					case(4)
						out_params%vdopp = DIRx(j)
					case(5)
						out_params%dtau = DIRx(j)
					case(6)
						out_params%delta_collision = DIRx(j)
					case(7)
						out_params%vmacro = DIRx(j)
					case(8)
						out_params%damping = DIRx(j)
					case(9)
						out_params%beta = DIRx(j)
					case(10)
						out_params%height = DIRx(j)
					case(11)
						out_params%dtau2 = DIRx(j)
					case(12)
						out_params%vmacro2 = DIRx(j)
					case(13)
						out_params%bgauss2 = DIRx(j)
					case(14)
						out_params%thetabd2 = DIRx(j)
					case(15)
						out_params%chibd2 = DIRx(j)
					case(16)
						out_params%vdopp2 = DIRx(j)
					case(17)
						out_params%ff = DIRx(j)
				end select
				j = j + 1
			endif
		enddo

		close(23)
		close(25)

! K-means clustering
! Read the sampling carried out by DIRECT
! 		open(unit=25,file=location_file_direct,action='read',status='old')
!
! 		close(25)
! 		call cluster_initialize_5 ( 2, point_num, cluster_num, point, &
! 			cluster_center )
!
! 		call kmeans_01 ( dim_num, point_num, cluster_num, it_max, it, point, &
! 			cluster, cluster_center, cluster_population, cluster_energy )

	end subroutine invert_with_direct

!------------------------------------------------------------
! Function that returns the chi^2 for the DIRECT problem. This solves the full RT problem
!------------------------------------------------------------
	subroutine fcn(n, x, f, flag, iidata, iisize, ddata, idsize, cdata, icsize)
	integer :: n,flag
   real(kind=8) :: x(n)
	real(kind=8) :: f
	integer :: iisize, idsize, icsize, i, j
	integer :: iidata(iisize)
	real(kind=8) :: ddata(idsize)
	character(len=40) :: cdata(icsize)
	type(variable_parameters) :: trial
	character(len=120) :: temporal_file

		temporal_file = 'temporal.prof'

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				select case(i)
					case(1)
						params%bgauss = x(j)
					case(2)
						params%thetabd = x(j)
					case(3)
						params%chibd = x(j)
					case(4)
						params%vdopp = x(j)
					case(5)
						params%dtau = x(j)
					case(6)
						params%delta_collision = x(j)
					case(7)
						params%vmacro = x(j)
					case(8)
						params%damping = x(j)
					case(9)
						params%beta = x(j)
					case(10)
						params%height = x(j)
					case(11)
						params%dtau2 = x(j)
					case(12)
						params%vmacro2 = x(j)
					case(13)
						params%bgauss2 = x(j)
					case(14)
						params%thetabd2 = x(j)
					case(15)
						params%chibd2 = x(j)
					case(16)
						params%vdopp2 = x(j)
					case(17)
						params%ff = x(j)
				end select
				j = j + 1
			endif
		enddo

		call do_synthesis(params, fixed, observation, inversion%stokes_unperturbed)
		call write_final_profiles(temporal_file, observation, inversion)

		f = compute_chisq(observation,inversion)

		call print_parameters(params,'      -Parameters : ',.TRUE.)

 		print *, 'chi^2 : ', f
		write(25,FMT='(11(E15.5,2X))') x, f

		flag = 0

	end subroutine fcn

!------------------------------------------------------------
! Function that returns the chi^2 for the DIRECT problem. This solves an approximate RT problem for Stokes I
!------------------------------------------------------------
	subroutine fcn_simplified_StokesI(n, x, f, flag, iidata, iisize, ddata, idsize, cdata, icsize)
	integer :: n,flag
   real(kind=8) :: x(n)
	real(kind=8) :: f
	integer :: iisize, idsize, icsize, i, j
	integer :: iidata(iisize)
	real(kind=8) :: ddata(idsize)
	character(len=40) :: cdata(icsize)
	type(variable_parameters) :: trial
	character(len=120) :: temporal_file
	real(kind=8), allocatable :: onum(:), prof(:), prof2(:)
	real(kind=8) :: va, dnum, adamp, onum0, onum1, onum2

		temporal_file = 'temporal.prof'

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				select case(i)
					case(1)
						params%bgauss = x(j)
					case(2)
						params%thetabd = x(j)
					case(3)
						params%chibd = x(j)
					case(4)
						params%vdopp = x(j)
					case(5)
						params%dtau = x(j)
					case(6)
						params%delta_collision = x(j)
					case(7)
						params%vmacro = x(j)
					case(8)
						params%damping = x(j)
					case(9)
						params%beta = x(j)
					case(10)
						params%height = x(j)
					case(11)
						params%dtau2 = x(j)
					case(12)
						params%vmacro2 = x(j)
					case(13)
						params%bgauss2 = x(j)
					case(14)
						params%thetabd2 = x(j)
					case(15)
						params%chibd2 = x(j)
					case(16)
						params%vdopp2 = x(j)
					case(17)
						params%ff = x(j)
				end select
				j = j + 1
			endif
		enddo

		allocate(onum(fixed%no))
		allocate(prof(fixed%no))

		onum = -1.d8 * observation%wl / fixed%wl**2

		inversion%stokes_unperturbed = 0.d0

! Only one slab
		if (params%nslabs == 1) then
			dnum = params%vdopp*1.d5 / (fixed%wl*1.d-8*PC)
			va = params%vmacro*1.d5 / (fixed%wl*1.d-8*PC)
			adamp = params%damping

 			if (fixed%damping_treatment == 1) then
 				adamp = fixed%wl * 1.d-8 / (params%vdopp*1.d5) * aesto(fixed%nemiss) * abs(params%damping)
 			endif

			onum0 = -1.26d0*1.d-8 / (fixed%wl*1.d-8)**2
			onum1 = (-1.26d0+0.09d0)*1.d-8 / (fixed%wl*1.d-8)**2
			onum2 = 0.d0
			inversion%stokes_unperturbed(0,:) = profile(adamp,(onum0-onum-va)/dnum) + &
				0.6*profile(adamp,(onum1-onum-va)/dnum) + &
				0.2*profile(adamp,(onum2-onum-va)/dnum)

			if (fixed%Stokes_incident(0) == 0) then
				inversion%stokes_unperturbed(0,:) = 1.d0 - exp(-params%dtau / 1.56d0 * inversion%stokes_unperturbed(0,:))
			else
				inversion%stokes_unperturbed(0,:) = exp(-params%dtau / 1.56d0 * inversion%stokes_unperturbed(0,:))
			endif
		endif

! Two slabs
 		if (params%nslabs == 2 .or. params%nslabs == 3) then
! First slab
			dnum = params%vdopp*1.d5 / (fixed%wl*1.d-8*PC)
			va = params%vmacro*1.d5 / (fixed%wl*1.d-8*PC)
			adamp = params%damping

 			if (fixed%damping_treatment == 1) then
 				adamp = fixed%wl * 1.d-8 / (params%vdopp*1.d5) * aesto(fixed%nemiss) * abs(params%damping)
 			endif

			onum0 = -1.26d0*1.d-8 / (fixed%wl*1.d-8)**2
			onum1 = (-1.26d0+0.09d0)*1.d-8 / (fixed%wl*1.d-8)**2
			onum2 = 0.d0
			prof = profile(adamp,(onum0-onum-va)/dnum) + &
				0.6*profile(adamp,(onum1-onum-va)/dnum) + &
				0.2*profile(adamp,(onum2-onum-va)/dnum)

			if (fixed%Stokes_incident(0) == 0) then
				inversion%stokes_unperturbed(0,:) = 1.d0 - exp(-params%dtau / 1.56d0 * prof)
			else
				inversion%stokes_unperturbed(0,:) = exp(-params%dtau / 1.56d0 * prof)
			endif

! Second slab
! If nslabs=3, then we use a different Doppler width. If not, we use the same for both components
			if (params%nslabs == 3) then
				dnum = params%vdopp2*1.d5 / (fixed%wl*1.d-8*PC)
				if (fixed%damping_treatment == 1) then
 				adamp = fixed%wl * 1.d-8 / (params%vdopp2*1.d5) * aesto(fixed%nemiss) * abs(params%damping)
 			endif

			endif
			va = params%vmacro2*1.d5 / (fixed%wl*1.d-8*PC)
			adamp = params%damping

			onum0 = -1.26d0*1.d-8 / (fixed%wl*1.d-8)**2
			onum1 = (-1.26d0+0.09d0)*1.d-8 / (fixed%wl*1.d-8)**2
			onum2 = 0.d0
			prof = profile(adamp,(onum0-onum-va)/dnum) + &
				0.6*profile(adamp,(onum1-onum-va)/dnum) + &
				0.2*profile(adamp,(onum2-onum-va)/dnum)

			inversion%stokes_unperturbed(0,:) = exp(-params%dtau2 / 1.56d0 * prof) * inversion%stokes_unperturbed(0,:)

 		endif

! Two slabs on the same pixel
		if (params%nslabs == -2) then

			allocate(prof2(fixed%no))

! First slab
			dnum = params%vdopp*1.d5 / (fixed%wl*1.d-8*PC)
			va = params%vmacro*1.d5 / (fixed%wl*1.d-8*PC)
			adamp = params%damping

 			if (fixed%damping_treatment == 1) then
 				adamp = fixed%wl * 1.d-8 / (params%vdopp*1.d5) * aesto(fixed%nemiss) * abs(params%damping)
 			endif

			onum0 = -1.26d0*1.d-8 / (fixed%wl*1.d-8)**2
			onum1 = (-1.26d0+0.09d0)*1.d-8 / (fixed%wl*1.d-8)**2
			onum2 = 0.d0
			prof = profile(adamp,(onum0-onum-va)/dnum) + &
				0.6*profile(adamp,(onum1-onum-va)/dnum) + &
				0.2*profile(adamp,(onum2-onum-va)/dnum)

! Second slab
			dnum = params%vdopp2*1.d5 / (fixed%wl*1.d-8*PC)
			va = params%vmacro2*1.d5 / (fixed%wl*1.d-8*PC)
			adamp = params%damping

 			if (fixed%damping_treatment == 1) then
 				adamp = fixed%wl * 1.d-8 / (params%vdopp2*1.d5) * aesto(fixed%nemiss) * abs(params%damping)
 			endif

			onum0 = -1.26d0*1.d-8 / (fixed%wl*1.d-8)**2
			onum1 = (-1.26d0+0.09d0)*1.d-8 / (fixed%wl*1.d-8)**2
			onum2 = 0.d0
			prof2 = profile(adamp,(onum0-onum-va)/dnum) + &
				0.6*profile(adamp,(onum1-onum-va)/dnum) + &
				0.2*profile(adamp,(onum2-onum-va)/dnum)

			if (fixed%Stokes_incident(0) == 0) then
				inversion%stokes_unperturbed(0,:) = params%ff * (1.d0 - exp(-params%dtau / 1.56d0 * prof)) + &
					(1.d0-params%ff) * (1.d0 - exp(-params%dtau2 / 1.56d0 * prof2))
			else
				inversion%stokes_unperturbed(0,:) = params%ff * exp(-params%dtau / 1.56d0 * prof) + &
					(1.d0-params%ff) * exp(-params%dtau2 / 1.56d0 * prof2)
			endif

			deallocate(prof2)
		endif

		f = compute_chisq(observation,inversion)

		deallocate(onum)
		deallocate(prof)

 		call write_final_profiles(temporal_file, observation, inversion)

		call print_parameters(params,'      -Parameters : ',.TRUE.)

 		print *, 'chi^2 : ', f
		write(25,FMT='(11(E15.5,2X))') x, f

		flag = 0

	end subroutine fcn_simplified_StokesI

!*************************************************************
!*************************************************************
! PIKAIA METHOD (GENETIC ALGORITHM)
!*************************************************************
!*************************************************************

!------------------------------------------------------------
! Invert some parameters with the GENETIC ALGORITHM (PIKAIA)
!------------------------------------------------------------
	subroutine invert_with_pikaia(out_params)
	integer :: n, seed, status, i, j
	real(kind=4) :: ctrl(12), f
	real(kind=4), allocatable :: x(:)
	type(variable_parameters) :: out_params

		n = sum(params%inverted)

! Use default values
		ctrl(1) = 100
		ctrl(2) = 5
		ctrl(3) = 5
		ctrl(4) = 0.85
		ctrl(5) = 2
		ctrl(6) = 0.005
		ctrl(7) = 0.0005
		ctrl(8) = 0.25
		ctrl(9) = 1.0
		ctrl(10) = 1
		ctrl(11) = 1
		ctrl(12) = 1

		seed = 1234
		call rninit(seed)

		allocate(x(n))

		call pikaia(chisq_pikaia,n,ctrl,x,f,status)

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				select case(i)
					case(1)
						out_params%bgauss = x(j)
					case(2)
						out_params%thetabd = x(j)
					case(3)
						out_params%chibd = x(j)
					case(4)
						out_params%vdopp = x(j)
					case(5)
						out_params%dtau = x(j)
					case(6)
						out_params%delta_collision = x(j)
					case(7)
						out_params%vmacro = x(j)
					case(8)
						out_params%damping = x(j)
					case(9)
						out_params%beta = x(j)
					case(10)
						out_params%height = x(j)
					case(11)
						out_params%dtau2 = x(j)
					case(12)
						out_params%vmacro2 = x(j)
				end select
				j = j + 1
			endif
		enddo

		deallocate(x)

	end subroutine invert_with_pikaia

!------------------------------------------------------------
! Function that returns the chisq for the PIKAIA algorithm
!------------------------------------------------------------
	function chisq_pikaia(n, in_x)
	real(kind=4) :: chisq_pikaia
	integer :: n, i, j
   real(kind=4) :: x(n), in_x(n), f

		j = 1

		do i = 1, params%n_total
			if (params%inverted(i) == 1) then

! Undo the change of variables that puts everything in the range [0,1]
				x(j) = minim_pikaia(i) - (minim_pikaia(i) - maxim_pikaia(i)) * in_x(i)

				select case(i)
					case(1)
						params%bgauss = x(j)
					case(2)
						params%thetabd = x(j)
					case(3)
						params%chibd = x(j)
					case(4)
						params%vdopp = x(j)
					case(5)
						params%dtau = x(j)
					case(6)
						params%delta_collision = x(j)
					case(7)
						params%vmacro = x(j)
					case(8)
						params%damping = x(j)
					case(9)
						params%beta = x(j)
					case(10)
						params%height = x(j)
					case(11)
						params%dtau2 = x(j)
					case(12)
						params%vmacro2 = x(j)
				end select
				j = j + 1
			endif
		enddo

		call do_synthesis(params, fixed, observation, inversion%stokes_unperturbed)

		f = compute_chisq(observation,inversion)

		call print_parameters(params,'      -Parameters : ',.TRUE.)

		chisq_pikaia = 1.d0 / f
		print *, 'chi^2 : ', f, chisq_pikaia

	end function chisq_pikaia

!*************************************************************
!*************************************************************
! SIMPLIFIED RADIATIVE TRANSFER FOR STOKES I
!*************************************************************
!*************************************************************

!------------------------------------------------------------
! Invert some parameters with the DIRECT method
!------------------------------------------------------------
	subroutine invert_with_simplified_transfer(out_params)
	real(kind=8), allocatable :: u(:), l(:), upper(:), lower(:)
	integer :: n, i, j, ndim, iprint, m
	real(kind=8), allocatable :: x(:), g(:), xold(:), wa(:)
	real(kind=8) :: f, factr, pgtol, fnew
	type(variable_parameters) :: out_params
	character(len=60)      :: task, csave
	logical                :: lsave(4)
	integer                :: isave(44)
	real(kind=8)               :: dsave(29)
	integer,  allocatable  :: nbd(:), iwa(:)

		factr = 1.d7
		pgtol = 1.d-5
		iprint = 0
		m = 5

		ndim = sum(params%inverted)
		allocate(u(ndim))
		allocate(l(ndim))
		allocate(x(ndim))
		allocate(xold(ndim))

		allocate(nbd(ndim))
		allocate(iwa(3*ndim))
		allocate(g(ndim))
		allocate(wa(2*m*n + 5*n + 11*m*m + 8*m))

		nbd = 2

		allocate(upper(17))
		allocate(lower(17))
		open(unit=24,file=direct_ranges,action='read',status='old')
		call lb(24,3)

		call lb(24,2)
		read(24,*)

		call lb(24,2)
		read(24,*)
		call lb(24,2)
		read(24,*)

		do i = 1, 17
			call lb(24,2)
			read(24,*) lower(i), upper(i)
		enddo
		close(24)

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				l(j) = lower(i)
				u(j) = upper(i)
				print *, 'Inverting ', trim(adjustl(parameters_name(i)))
				j = j + 1
			endif
		enddo

		deallocate(upper)
		deallocate(lower)

! Invert the Stokes I profile with the simplified RT method

		task = 'START'

		do while(task(1:2) == 'FG' .or. task == 'NEW_X' .or. task == 'START')

!     This is the call to the L-BFGS-B code.
         call setulb ( ndim, m, x, l, u, nbd, f, g, factr, pgtol, &
                       wa, iwa, task, iprint,&
                       csave, lsave, isave, dsave )

! Evaluate the forward model and the gradient
         if (task(1:2) == 'FG') then

				call simplified_transfer(ndim, x, f)

				xold = x

! Calculate gradient numerically
				do i = 1, ndim
					x(i) = x(i) + 1.d-3
					call simplified_transfer(ndim, x, fnew)
					g(i) = (fnew - f) / 1.d-3
					x(i) = xold(i)
				enddo

				print *, 'f,g= ', f, g

         endif
      enddo

      stop

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				select case(i)
					case(1)
						out_params%bgauss = x(j)
					case(2)
						out_params%thetabd = x(j)
					case(3)
						out_params%chibd = x(j)
					case(4)
						out_params%vdopp = x(j)
					case(5)
						out_params%dtau = x(j)
					case(6)
						out_params%delta_collision = x(j)
					case(7)
						out_params%vmacro = x(j)
					case(8)
						out_params%damping = x(j)
					case(9)
						out_params%beta = x(j)
					case(10)
						out_params%height = x(j)
					case(11)
						out_params%dtau2 = x(j)
					case(12)
						out_params%vmacro2 = x(j)
					case(13)
						out_params%bgauss2 = x(j)
					case(14)
						out_params%thetabd2 = x(j)
					case(15)
						out_params%chibd2 = x(j)
					case(16)
						out_params%vdopp2 = x(j)
					case(17)
						out_params%ff = x(j)
				end select
				j = j + 1
			endif
		enddo

	end subroutine invert_with_simplified_transfer

!------------------------------------------------------------
! Function that returns the the DIRECT
!------------------------------------------------------------
	subroutine simplified_transfer(n, x, f)
	integer :: n, i, j
   real(kind=8) :: x(n)
	real(kind=8) :: f, g(n)
	type(variable_parameters) :: trial
	character(len=120) :: temporal_file
	real(kind=8), allocatable :: onum(:)
	real(kind=8) :: va, dnum, adamp, onum0, onum1, onum2

		j = 1
		do i = 1, params%n_total
			if (params%inverted(i) == 1) then
				select case(i)
					case(1)
						params%bgauss = x(j)
					case(2)
						params%thetabd = x(j)
					case(3)
						params%chibd = x(j)
					case(4)
						params%vdopp = x(j)
					case(5)
						params%dtau = x(j)
					case(6)
						params%delta_collision = x(j)
					case(7)
						params%vmacro = x(j)
					case(8)
						params%damping = x(j)
					case(9)
						params%beta = x(j)
					case(10)
						params%height = x(j)
					case(11)
						params%dtau2 = x(j)
					case(12)
						params%vmacro2 = x(j)
					case(13)
						params%bgauss2 = x(j)
					case(14)
						params%thetabd2 = x(j)
					case(15)
						params%chibd2 = x(j)
					case(16)
						params%vdopp2 = x(j)
					case(17)
						params%ff = x(j)
				end select
				j = j + 1
			endif
		enddo

		allocate(onum(fixed%no))

		onum = -1.d8 * observation%wl / fixed%wl**2

		if (params%nslabs == 1) then
			dnum = params%vdopp*1.d5 / (fixed%wl*1.d-8*PC)
			va = params%vmacro*1.d5 / (fixed%wl*1.d-8*PC)
			adamp = params%damping

			onum0 = 0.d0
			onum1 = -0.09d0*1.d-8 / (fixed%wl*1.d-8)**2
			onum2 = -1.26d0*1.d-8 / (fixed%wl*1.d-8)**2
			inversion%stokes_unperturbed(0,:) = profile(adamp,(onum0-onum-va)/dnum) + &
				0.6*profile(adamp,(onum1-onum-va)/dnum) + &
				0.2*profile(adamp,(onum2-onum-va)/dnum)

			if (fixed%Stokes_incident(0) == 0) then
				inversion%stokes_unperturbed(0,:) = 1.d0 - exp(-params%dtau * inversion%stokes_unperturbed(0,:))
			else
				inversion%stokes_unperturbed(0,:) = exp(-params%dtau * inversion%stokes_unperturbed(0,:))
			endif

			print *, 'x,f = ' ,x, f

			f = compute_chisq(observation,inversion)

		endif

		deallocate(onum)

		call print_parameters(params,'      -Parameters : ',.TRUE.)

 		print *, 'chi^2 : ', f
		write(25,FMT='(11(E15.5,2X))') x, f

	end subroutine simplified_transfer

end module marquardt

max_stars_repo_name: aasensio/pyiacsun

max_issues_repo_name: aasensio/pyiacsun

max_forks_repo_name: aasensio/pyiacsun


text:
Lake Spafford isnt a ponds lake so much as a place south of Mrak Hall in the Arboretum where Putah Creek becomes obese and then thins down again. It is named after Ed Spafford, who was ASUCD President in 1953 and Vice Chancellor for many years(????).  Also known as the location of the Picnic Day Battle of the Bands, and a good place to find kids feeding Ducks ducks.

On Sunday afternoon, August 17th, 2014, Eric Jensen and his wife saw something they had never seen beforeotters! They were a group of about 68 river otters that swam into Lake Spafford from the west end of the Lake. The sight became even more remarkable and a bit gruesome when they gathered on the southern side of the lake to attack and ultimately kill an adult duck. Several ducks flew around the lake in panic. It was too late, and minutes after killing the duck the group of otters began to tear apart the ducks body before a number of horrified visitors to the lake. What a sight!

 Plants found around the Lake
  RicePaper Plant

 Benches, Memorials, Trees

  Allean Burton Service Tree
  Rosie Gonzales Bench
  Gerry Taylor Bench
  Janet C. Hamilton Service Bench
  Warren Roberts Birthday Tree

 Photos

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!==============================================================================!
  subroutine Multiphase_Mod_Compute_Vof(mult, sol, dt, n)
!------------------------------------------------------------------------------!
!   Solves Volume Fraction equation using UPWIND ADVECTION and CICSAM          !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Work_Mod, only: beta_f => r_face_01, c_d => r_cell_01
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Multiphase_Type), target :: mult
  type(Solver_Type),     target :: sol
  real                          :: dt
  integer                       :: n    ! current temporal iteration
!-----------------------------------[Locals]-----------------------------------!
  type(Field_Type),  pointer :: flow
  type(Grid_Type),   pointer :: grid
  type(Var_Type),    pointer :: vof
  type(Face_Type),   pointer :: m_flux
  real, contiguous,  pointer :: vof_f(:)
  type(Matrix_Type), pointer :: a
  real, contiguous,  pointer :: b(:)
  integer                    :: s, c, c1, c2
  integer                    :: donor, accept, corr_num
  integer                    :: i_sub, n_sub
  real                       :: fs
  real                       :: courant_max, epsloc
  real,              pointer :: courant_max_param
  integer,           pointer :: n_sub_param, corr_num_max
  character(len=80)          :: solver
!==============================================================================!

  call Cpu_Timer_Mod_Start('Compute_Multiphase (without solvers)')

  ! Take aliases
  flow   => mult % pnt_flow
  grid   => flow % pnt_grid
  m_flux => flow % m_flux
  vof    => mult % vof
  vof_f  => mult % vof_f
  courant_max_param => mult % courant_max_param
  n_sub_param       => mult % n_sub_param
  corr_num_max      => mult % corr_num_max

  a => sol % a
  b => sol % b % val

  epsloc = epsilon(epsloc)

  if (vof % adv_scheme .eq. CICSAM .or. &
      vof % adv_scheme .eq. STACS) then

    ! Compute courant Number close to the interface:
    call Vof_Max_Courant_Number(mult, dt, c_d, 1, courant_max)

    n_sub = min(max(ceiling(courant_max / courant_max_param),1),n_sub_param)

  else
    ! Old volume fraction:
    vof % o(:) = vof % n(:)
  end if

  if (vof % adv_scheme .eq. UPWIND) then
    !-------------------------!
    !   Matrix Coefficients   !
    !-------------------------!

    call Multiphase_Mod_Vof_Coefficients(mult, a, b, dt, beta_f)

    ! Solve System
    call Multiphase_Mod_Vof_Solve_System(mult, sol, b)

    call Grid_Mod_Exchange_Real(grid, vof % n)

    !-----------------------------!
    !   Correct Volume Fraction   !
    !-----------------------------!

    do c = 1, grid % n_cells
      vof % n(c) = max(min(vof % n(c),1.0),0.0)
    end do

  else if (vof % adv_scheme .eq. CICSAM .or. &
           vof % adv_scheme .eq. STACS) then

    do i_sub = 1, n_sub

      ! Courant number full domain:
      call Vof_Max_Courant_Number(mult, dt / real(n_sub),    &
                                  c_d, 0, courant_max)

      ! Old volume fraction:
      vof % o(:) = vof % n(:)

      !---------------------------!
      !   Predict Beta at faces   !
      !---------------------------!
      ! Compute Gradient:
      call Field_Mod_Grad_Variable(flow, vof)

      call Multiphase_Mod_Vof_Predict_Beta(vof,                        &
                                           m_flux % n,                 &
                                           vof % x, vof % y, vof % z,  &
                                           grid % dx,                  &
                                           grid % dy,                  &
                                           grid % dz,                  &
                                           beta_f,                     &
                                           c_d)

      loop_corr:  do corr_num = 1, corr_num_max
        !-------------------------!
        !   Matrix Coefficients   !
        !-------------------------!

        call Multiphase_Mod_Vof_Coefficients(mult, a, b,         &
                                             dt / real(n_sub),   &
                                             beta_f)

        ! Solve System
        call Multiphase_Mod_Vof_Solve_System(mult, sol, b)

        call Grid_Mod_Exchange_Real(grid, vof % n)

        !---------------------------!
        !   Correct Beta at faces   !
        !---------------------------!
        call Multiphase_Mod_Vof_Correct_Beta(mult,               &
                                             vof,                &
                                             m_flux % n,         &
                                             beta_f,             &
                                             dt / real(n_sub))

        !------------------------!
        !   Correct Boundaries   !
        !------------------------!
        do s = 1, grid % n_faces

          c1 = grid % faces_c(1,s)
          c2 = grid % faces_c(2,s)

          if (c2 < 0) then
            if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. OUTFLOW) then
              vof % n(c2) = max(min(vof % n(c1),1.0),0.0)
            end if
            if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. WALL) then
              vof % n(c2) = max(min(vof % n(c1),1.0),0.0)
            end if
          end if

        end do

        !---------------------!
        !   Noise Reduction   !
        !---------------------!
        do c = 1, grid % n_cells
          if(vof % n(c) < epsloc) then
            vof % n(c) = 0.0
          end if

          if(1.0 - vof % n(c) < epsloc) then
            vof % n(c) = 1.0
          end if
        end do

      end do loop_corr

      !-----------------------------!
      !   Correct Volume Fraction   !
      !-----------------------------!
      do c = 1, grid % n_cells
        vof % n(c) = max(min(vof % n(c),1.0),0.0)
      end do

      call Grid_Mod_Exchange_Real(grid, vof % n)
    end do

  end if

  !------------------------------!
  !   Volume fraction at faces   !
  !------------------------------!

  if (vof % adv_scheme .eq. UPWIND) then

    do s = 1, grid % n_faces
      c1 = grid % faces_c(1,s)
      c2 = grid % faces_c(2,s)
      fs = grid % f(s)

      ! Face is inside the domain
      if(c2 > 0) then

        if (m_flux % n(s)>=0.0) then
          vof_f(s) = vof % n(c1)
        else
          vof_f(s) = vof % n(c2)
        end if

      ! Side is on the boundary
      else ! (c2 < 0)

        vof_f(s) = vof % n(c1)

      end if

    end do

  else if (vof % adv_scheme .eq. CICSAM .or. vof % adv_scheme .eq. STACS) then

    do s = 1, grid % n_faces
      c1 = grid % faces_c(1,s)
      c2 = grid % faces_c(2,s)
      fs = grid % f(s)

      if(c2 > 0) then
        if(m_flux % n(s) >= 0) then
          donor = c1
          accept = c2
        else
          donor = c2
          accept = c1
        end if

        vof_f(s) = 0.5 * ((1.0 - beta_f(s)) * ( vof % n(donor)      &
                                              + vof % o(donor) )    &
                               + beta_f(s)  * ( vof % n(accept)     &
                                              + vof % o(accept) ))

      else
        if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. OUTFLOW) then
          vof_f(s) = vof % n(c1)
        else if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. INFLOW) then
          vof_f(s) = vof % n(c2)
        else
          vof_f(s) = vof % n(c1)
        end if

      end if
    end do

  end if

  !-----------------------!
  !   Update properties   !
  !-----------------------!
  call Multiphase_Mod_Update_Physical_Properties(mult)

  call Field_Mod_Grad_Variable(flow, vof)

  !----------------------------------------!
  !   Surface Tension Force Contribution   !
  !----------------------------------------!

  ! If distance function is calculated
  if (mult % d_func) then
    call Multiphase_Mod_Compute_Distance_Function(mult, sol, flow % dt, n)
  end if

  if (mult % surface_tension > TINY ) then
    call Multiphase_Mod_Vof_Surface_Tension_Contribution(mult)
  end if

  call Cpu_Timer_Mod_Stop('Compute_Multiphase (without solvers)')

  end subroutine

max_stars_repo_name: EdinSmartLab/T-Flows

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                          .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... file islapes.f
c
c     this file includes documentation and code for
c     subroutine islapes         i
c
c ... files which must be loaded with islapes.f
c
c     sphcom.f, hrfft.f, shaes.f, shses.f
c
c     subroutine islapes(nlat,nlon,isym,nt,xlmbda,sf,ids,jds,a,b,
c    +mdab,ndab,wshses,lshses,work,lwork,pertrb,ierror)
c
c     islapes inverts the laplace or helmholz operator on an equally
c     spaced latitudinal grid using o(n**3) storage. given the
c     spherical harmonic coefficients a(m,n) and b(m,n) of the right
c     hand side slap(i,j), islapes computes a solution sf(i,j) to
c     the following helmhotz equation :
c
c           2                2
c     [d(sf(i,j))/dlambda /sint + d(sint*d(sf(i,j))/dtheta)/dtheta]/sint
c
c                   - xlmbda * sf(i,j) = slap(i,j)
c
c      where sf(i,j) is computed at colatitude
c
c                 theta(i) = (i-1)*pi/(nlat-1)
c
c            and longitude
c
c                 lambda(j) = (j-1)*2*pi/nlon
c
c            for i=1,...,nlat and j=1,...,nlon.
c
c
c     input parameters
c
c     nlat   the number of colatitudes on the full sphere including the
c            poles. for example, nlat = 37 for a five degree grid.
c            nlat determines the grid increment in colatitude as
c            pi/(nlat-1).  if nlat is odd the equator is located at
c            grid point i=(nlat+1)/2. if nlat is even the equator is
c            located half way between points i=nlat/2 and i=nlat/2+1.
c            nlat must be at least 3. note: on the half sphere, the
c            number of grid points in the colatitudinal direction is
c            nlat/2 if nlat is even or (nlat+1)/2 if nlat is odd.
c
c     nlon   the number of distinct longitude points.  nlon determines
c            the grid increment in longitude as 2*pi/nlon. for example
c            nlon = 72 for a five degree grid. nlon must be greater
c            than zero. the axisymmetric case corresponds to nlon=1.
c            the efficiency of the computation is improved when nlon
c            is a product of small prime numbers.
c
C*PL*ERROR* Comment line too long
c     isym   this parameter should have the same value input to subroutine
C*PL*ERROR* Comment line too long
c            shaes to compute the coefficients a and b for the scalar field
c            slap.  isym is set as follows:
c
c            = 0  no symmetries exist in slap about the equator. scalar
c                 synthesis is used to compute sf on the entire sphere.
c                 i.e., in the array sf(i,j) for i=1,...,nlat and
c                 j=1,...,nlon.
c
c           = 1  sf and slap are antisymmetric about the equator. the
c                synthesis used to compute sf is performed on the
c                northern hemisphere only.  if nlat is odd, sf(i,j) is
c                computed for i=1,...,(nlat+1)/2 and j=1,...,nlon.  if
c                nlat is even, sf(i,j) is computed for i=1,...,nlat/2
c                and j=1,...,nlon.
c
c
c           = 2  sf and slap are symmetric about the equator. the
c                synthesis used to compute sf is performed on the
c                northern hemisphere only.  if nlat is odd, sf(i,j) is
c                computed for i=1,...,(nlat+1)/2 and j=1,...,nlon.  if
c                nlat is even, sf(i,j) is computed for i=1,...,nlat/2
c                and j=1,...,nlon.
c
c
c   nt       the number of solutions. in the program that calls islapes
c            the arrays sf,a, and b can be three dimensional in which
c            case multiple solutions are computed. the third index
c            is the solution index with values k=1,...,nt.
c            for a single solution set nt=1. the description of the
c            remaining parameters is simplified by assuming that nt=1
c            and sf,a,b are two dimensional.
c
c   xlmbda   a one dimensional array with nt elements. if xlmbda is
c            is identically zero islapes solves poisson's equation.
c            if xlmbda > 0.0 islapes solves the helmholtz equation.
c            if xlmbda < 0.0 the nonfatal error flag ierror=-1 is
c            returned. negative xlambda could result in a division
c            by zero.
c
c   ids      the first dimension of the array sf as it appears in the
C*PL*ERROR* Comment line too long
c            program that calls islapes.  if isym = 0 then ids must be at
c            least nlat.  if isym > 0 and nlat is even then ids must be
c            at least nlat/2. if isym > 0 and nlat is odd then ids must
c            be at least (nlat+1)/2.
c
c   jds      the second dimension of the array sf as it appears in the
c            program that calls islapes. jds must be at least nlon.
c
c
c   a,b      two or three dimensional arrays (see input parameter nt)
c            that contain scalar spherical harmonic coefficients
c            of the scalar field slap. a,b must be computed by shaes
c            prior to calling islapes.
c
c
c   mdab     the first dimension of the arrays a and b as it appears
c            in the program that calls islapes.  mdab must be at
c            least min0(nlat,(nlon+2)/2) if nlon is even or at least
c            min0(nlat,(nlon+1)/2) if nlon is odd.
c
c   ndab     the second dimension of the arrays a and b as it appears
c            in the program that calls islapes. ndbc must be at least
c            least nlat.
c
c            mdab,ndab should have the same values input to shaes to
c            compute the coefficients a and b.
c
c
c   wshses   an array which must be initialized by subroutine shsesi.
c            once initialized, wshses can be used repeatedly by
c            islapes as long as nlat and nlon  remain unchanged.
c            wshses must not be altered between calls of islapes.
c
c    lshses  the dimension of the array wshses as it appears in the
c            program that calls islapes.  let
c
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            and
c
c               l2 = nlat/2        if nlat is even or
c               l2 = (nlat+1)/2    if nlat is odd
c
c            then lshses must be at least
c
c               (l1*l2*(nlat+nlat-l1+1))/2+nlon+15
c
c     work   a work array that does not have to be saved.
c
c     lwork  the dimension of the array work as it appears in the
c            program that calls islapes. define
c
c               l2 = nlat/2                    if nlat is even or
c               l2 = (nlat+1)/2                if nlat is odd
c               l1 = min0(nlat,(nlon+2)/2) if nlon is even or
c               l1 = min0(nlat,(nlon+1)/2) if nlon is odd
c
c            if isym is zero then lwork must be at least
c
c               (nt+1)*nlat*nlon + nlat*(2*nt*l1+1)
c
c            if isym is nonzero lwork must be at least
c
c               (nt+1)*l2*nlon + nlat*(2*nt*l1+1)
c
c
c     **************************************************************
c
c     output parameters
c
c
C*PL*ERROR* Comment line too long
c    sf      a two or three dimensional arrays (see input parameter nt) that
C*PL*ERROR* Comment line too long
c            inverts the scalar laplacian in slap - pertrb.  sf(i,j) is given
c            at the colatitude
c
c                 theta(i) = (i-1)*pi/(nlat-1)
c
c            and longitude
c
c                 lambda(j) = (j-1)*2*pi/nlon
c
c            for i=1,...,nlat and j=1,...,nlon.
c
c   pertrb  a one dimensional array with nt elements (see input
c           parameter nt). in the discription that follows we assume
c           that nt=1. if xlmbda > 0.0 then pertrb=0.0 is always
c           returned because the helmholtz operator is invertible.
c           if xlmbda = 0.0 then a solution exists only if a(1,1)
c           is zero. islapec sets a(1,1) to zero. the resulting
c           solution sf(i,j) solves poisson's equation with
c           pertrb = a(1,1)/(2.*sqrt(2.)) subtracted from the
c           right side slap(i,j).
c
C*PL*ERROR* Comment line too long
c  ierror    a parameter which flags errors in input parameters as follows:
c
c            = 0  no errors detected
c
c            = 1  error in the specification of nlat
c
c            = 2  error in the specification of nlon
c
c            = 3  error in the specification of ityp
c
c            = 4  error in the specification of nt
c
c            = 5  error in the specification of ids
c
c            = 6  error in the specification of jds
c
c            = 7  error in the specification of mdbc
c
c            = 8  error in the specification of ndbc
c
c            = 9  error in the specification of lshses
c
c            = 10 error in the specification of lwork
c
c
c **********************************************************************
c
c     end of documentation for islapes
c
c **********************************************************************
c
      SUBROUTINE DISLAPES(NLAT,NLON,ISYM,NT,XLMBDA,SF,IDS,JDS,A,B,MDAB,
     +                   NDAB,WSHSES,LSHSES,WORK,LWORK,PERTRB,IERROR)
      DOUBLE PRECISION XLMBDA
      DOUBLE PRECISION SF
      DOUBLE PRECISION A
      DOUBLE PRECISION B
      DOUBLE PRECISION WSHSES
      DOUBLE PRECISION WORK
      DOUBLE PRECISION PERTRB
      DIMENSION SF(IDS,JDS,NT),A(MDAB,NDAB,NT),B(MDAB,NDAB,NT)
      DIMENSION WSHSES(LSHSES),WORK(LWORK),XLMBDA(NT),PERTRB(NT)
c
c     check input parameters
c
      IERROR = 1
      IF (NLAT.LT.3) RETURN
      IERROR = 2
      IF (NLON.LT.4) RETURN
      IERROR = 3
      IF (ISYM.LT.0 .OR. ISYM.GT.2) RETURN
      IERROR = 4
      IF (NT.LT.0) RETURN
      IERROR = 5
      IMID = (NLAT+1)/2
      IF ((ISYM.EQ.0.AND.IDS.LT.NLAT) .OR.
     +    (ISYM.GT.0.AND.IDS.LT.IMID)) RETURN
      IERROR = 6
      IF (JDS.LT.NLON) RETURN
      IERROR = 7
      MMAX = MIN0(NLAT,NLON/2+1)
      IF (MDAB.LT.MMAX) RETURN
      IERROR = 8
      IF (NDAB.LT.NLAT) RETURN
      IERROR = 9
c
c     set and verify saved work space length
c
      IMID = (NLAT+1)/2
      LPIMN = (IMID*MMAX* (NLAT+NLAT-MMAX+1))/2
      IF (LSHSES.LT.LPIMN+NLON+15) RETURN
      IERROR = 10
c
c     set and verify unsaved work space length
c
      LS = NLAT
      IF (ISYM.GT.0) LS = IMID
      NLN = NT*LS*NLON
      MN = MMAX*NLAT*NT
c     lwkmin = nln+ls*nlon+2*mn+nlat
c     if (lwork .lt. lwkmin) return
      L2 = (NLAT+1)/2
      L1 = MIN0(NLAT,NLON/2+1)
      IF (ISYM.EQ.0) THEN
          LWKMIN = (NT+1)*NLAT*NLON + NLAT* (2*NT*L1+1)
      ELSE
          LWKMIN = (NT+1)*L2*NLON + NLAT* (2*NT*L1+1)
      END IF
      IF (LWORK.LT.LWKMIN) RETURN
      IERROR = 0
c
c     check sign of xlmbda
c
      DO K = 1,NT
          IF (XLMBDA(K).LT.0.0D0) THEN
              IERROR = -1
          END IF
      END DO
c
c     set work space pointers
c
      IA = 1
      IB = IA + MN
      IFN = IB + MN
      IWK = IFN + NLAT
      LWK = LWORK - 2*MN - NLAT
      CALL DISLPES1(NLAT,NLON,ISYM,NT,XLMBDA,SF,IDS,JDS,A,B,MDAB,NDAB,
     +             WORK(IA),WORK(IB),MMAX,WORK(IFN),WSHSES,LSHSES,
     +             WORK(IWK),LWK,PERTRB,IERROR)
      RETURN
      END

      SUBROUTINE DISLPES1(NLAT,NLON,ISYM,NT,XLMBDA,SF,IDS,JDS,A,B,MDAB,
     +                   NDAB,AS,BS,MMAX,FNN,WSHSES,LSHSES,WK,LWK,
     +                   PERTRB,IERROR)
      DOUBLE PRECISION XLMBDA
      DOUBLE PRECISION SF
      DOUBLE PRECISION A
      DOUBLE PRECISION B
      DOUBLE PRECISION AS
      DOUBLE PRECISION BS
      DOUBLE PRECISION FNN
      DOUBLE PRECISION WSHSES
      DOUBLE PRECISION WK
      DOUBLE PRECISION PERTRB
      DOUBLE PRECISION FN
      DIMENSION SF(IDS,JDS,NT),A(MDAB,NDAB,NT),B(MDAB,NDAB,NT)
      DIMENSION AS(MMAX,NLAT,NT),BS(MMAX,NLAT,NT),FNN(NLAT)
      DIMENSION WSHSES(LSHSES),WK(LWK),PERTRB(NT),XLMBDA(NT)
c
c     set multipliers and preset synthesis coefficients to zero
c
      DO N = 1,NLAT
          FN = DBLE(N-1)
          FNN(N) = FN* (FN+1.0D0)
          DO M = 1,MMAX
              DO K = 1,NT
                  AS(M,N,K) = 0.0D0
                  BS(M,N,K) = 0.0D0
              END DO
          END DO
      END DO
      DO K = 1,NT
c
c     compute synthesis coefficients for xlmbda zero or nonzero
c
          IF (XLMBDA(K).EQ.0.0D0) THEN
              DO N = 2,NLAT
                  AS(1,N,K) = -A(1,N,K)/FNN(N)
                  BS(1,N,K) = -B(1,N,K)/FNN(N)
              END DO
              DO M = 2,MMAX
                  DO N = M,NLAT
                      AS(M,N,K) = -A(M,N,K)/FNN(N)
                      BS(M,N,K) = -B(M,N,K)/FNN(N)
                  END DO
              END DO
          ELSE
c
c     xlmbda nonzero so operator invertible unless
c     -n*(n-1) = xlmbda(k) < 0.0  for some n
c
              PERTRB(K) = 0.0D0
              DO N = 1,NLAT
                  AS(1,N,K) = -A(1,N,K)/ (FNN(N)+XLMBDA(K))
                  BS(1,N,K) = -B(1,N,K)/ (FNN(N)+XLMBDA(K))
              END DO
              DO M = 2,MMAX
                  DO N = M,NLAT
                      AS(M,N,K) = -A(M,N,K)/ (FNN(N)+XLMBDA(K))
                      BS(M,N,K) = -B(M,N,K)/ (FNN(N)+XLMBDA(K))
                  END DO
              END DO
          END IF
      END DO
c
c     synthesize as,bs into sf
c
      CALL DSHSES(NLAT,NLON,ISYM,NT,SF,IDS,JDS,AS,BS,MMAX,NLAT,WSHSES,
     +           LSHSES,WK,LWK,IERROR)
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! From pwsfC:
! array-valued functions whose size depends on the size of a
! an argument whose correspoinding actual is a module allocatable
! array.
	module fn16
	real, allocatable :: zzz(:,:)
	contains
	    function ff1(xx)
	    real :: xx(:,:)
	    real :: ff1(size(xx))
	    ff1 = 11
	    endfunction
	    function ff2(xx)
	    real :: xx(:,:)
	    real :: ff2(size(xx,1))
	    ff2 = 13
	    endfunction
	    function ff3(xx,i)
	    real :: xx(:,:)
	    real :: ff3(size(xx,i))
	    ff3 = 17
	    endfunction
	endmodule
	subroutine sub(ii)
	use fn16
	integer :: result(3)
	integer :: expect(3)=(/66,26,51/)
	allocate(zzz(2,3))
	result(1) = sum(ff1(zzz))
	result(2) = sum(ff2(zzz))
	result(3) = sum(ff3(zzz,ii))
	call check(result, expect, 3)
	end
	call sub(2)
	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
c $Id: polfit.f,v 1.1 1996/04/19 17:12:51 ty7777 Exp $
c $Log: polfit.f,v $
c Revision 1.1  1996/04/19  17:12:51  ty7777
c Initial revision
c
c
      SUBROUTINE POLFIT(X,Y,NPTS,A,ARRAY,NTERM,CHISQR)
C     BEVINGTON p. 141
      character*256 rcsid /'$Header: /data/petsun4/src_solaris/librms/RCS/polfit.f,v 1.1 1996/04/19 17:12:51 ty7777 Exp $'/
      REAL*4 X(NPTS),Y(NPTS)
      REAL*4 A(NTERM),ARRAY(NTERM,NTERM)
      REAL*4 SUMX(19),SUMY(10)
      IF(NTERM.GT.10)STOP 'POLFIT NTERM TOO LARGE'
      NMAX=2*NTERM-1
      DO 13 N=1,NMAX
   13 SUMX(N)=0.
      DO 15 J=1,NTERM
   15 SUMY(J)=0.
      CHISQR=0.
      DO 50 I=1,NPTS
      XTERM=1.
      DO 44 N=1,NMAX
      SUMX(N)=SUMX(N)+XTERM
   44 XTERM=XTERM*X(I)
      YTERM=Y(I)
      DO 48 N=1,NTERM
      SUMY(N)=SUMY(N)+YTERM
   48 YTERM=YTERM*X(I)
      CHISQR=CHISQR+Y(I)**2
   50 CONTINUE
      DO 54 J=1,NTERM
      DO 54 K=1,NTERM
   54 ARRAY(J,K)=SUMX(J+K-1)
      if(.true.)then
      CALL MATINV(ARRAY,NTERM,DET)
      DO 61 L=1,NTERM
      A(L)=0.
      do 61 J=1,NTERM
   61 A(L)=A(L)+ARRAY(L,J)*SUMY(J)
      else
      DELTA=DETERM(ARRAY,NTERM)
      DO 70 L=1,NTERM
      DO 66 J=1,NTERM
      DO 65 K=1,NTERM
   65 ARRAY(J,K)=SUMX(J+K-1)
   66 ARRAY(J,L)=SUMY(J)
   70 A(L)=DETERM(ARRAY,NTERM)/DELTA
      endif
      DO 75 J=1,NTERM
      CHISQR=CHISQR-2.*A(J)*SUMY(J)
      DO 75 K=1,NTERM
   75 CHISQR=CHISQR+A(J)*A(K)*SUMX(J+K-1)
      CHISQR=CHISQR/FLOAT(NPTS-NTERM)
      RETURN
      END
      FUNCTION DETERM(ARRAY,NORDER)
      REAL*4 ARRAY(NORDER,NORDER)
      DETERM=1.
      DO 50 K=1,NORDER
      IF(ARRAY(K,K).NE.0.)GOTO 41
      DO 23 J=K,NORDER
      IF(ARRAY(J,K).NE.0.)GOTO 31
   23 CONTINUE
      DETERM=0.
      RETURN
   31 DO 34 I=K,NORDER
      SAVE=ARRAY(I,J)
      ARRAY(I,J)=ARRAY(I,K)
   34 ARRAY(I,K)=SAVE
      DETERM=-DETERM
   41 DETERM=DETERM*ARRAY(K,K)
      IF(K.GE.NORDER)GOTO 50
      DO 46 I=K+1,NORDER
      DO 46 J=K+1,NORDER
   46 ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE POLTST
      PARAMETER (NPTS=51)
      REAL*4 X(NPTS),Y(NPTS)
      PARAMETER (NORDER=10)
      REAL*4 A(NORDER),ARRAY(NORDER,NORDER),C(NORDER)
      write(*,"(' ENTER NTERM')")
      read(*,*)NTERM
      IF(NTERM.GT.NORDER)STOP 'POLTST NTERM TOO LARGE'
      write(*,"(' ENTER NTERM COEFFICIENTS')")
      read(*,*)(C(N),N=1,NTERM)
      DO 50 I=1,NPTS
      T=-1.+2.*FLOAT(I-1)/FLOAT(NPTS-1)
      X(I)=T
      Y(I)=C(1)
      DO 2 N=2,NTERM
      Y(I)=Y(I)+C(N)*T
    2 T=T*X(I)
   50 CONTINUE
      CALL POLFIT(X,Y,NPTS,A,ARRAY,NTERM,CHISQR)
      DO 10 N=1,NTERM
   10 write(*,101)N,C(N),A(N)
  101 FORMAT(I3,2F10.4)
      write(*,102)CHISQR
  102 FORMAT(' CHISQR',F10.4)
      STOP
      END

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP


text:
! { dg-do compile }
! { dg-options "-std=legacy" }
!
! PR 33269: we used to not simplify format strings before checking if
! they were valid, leading to a missed error.

IMPLICIT CHARACTER*5 (h-z)

CHARACTER*5 f
CHARACTER*5 bad, good
parameter(bad="a", good="(a)")

PRINT ('a'), "hello" ! { dg-error "Missing leading left parenthesis in format string" }
WRITE (*, ("a")) "error"  ! { dg-error "Missing leading left parenthesis in format string" }

PRINT 'a', "hello" ! { dg-error "Missing leading left parenthesis in format string" }
WRITE (*, "a") "error"  ! { dg-error "Missing leading left parenthesis in format string" }
WRITE (*, bad) "error"  ! { dg-error "Missing leading left parenthesis in format string" }

PRINT 'a' // ', a', "err", "or"   ! { dg-error "Missing leading left parenthesis in format string" }

PRINT '(' // 'a'  ! { dg-error "Unexpected end of format string in format string" }

! the following are ok
PRINT "(2f5.3)", bar, foo
PRINT ' (a)', "hello"
WRITE (*, " ((a))") "hello"
print "(a" // ")", "all is fine"
print good, "great"

! verify that we haven't broken non-constant expressions
f = "(f5.3)"
print f, 3.14159
print (f), 2.71813
print implicitly_typed, "something"
write (*, implicitly_typed_as_well) "something else"
END

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
subroutine PGG_AroundUpdateNode_RegularGraph(APlayerPayoff,BPlayerPayoff,nodeChosenUpdate,payoffUpdate, networkMain)
    use global
    implicit none
    real,intent(in)::APlayerPayoff(numPlayerInGame)
    real,intent(in)::BPlayerPayoff(numPlayerInGame)
    type(structuredPopulations),intent(in)::networkMain(Ntotal)
    integer,intent(in)::nodeChosenUpdate
    real,intent(out)::payoffUpdate
    !!
    integer::numAnodeUpdatePGG, numANeighborPGG
    integer::nodeNeighNow, indexNeigh
    integer::strategyUpdateNode, strategyNeighNode

    ! Compute the payoff of the node Chosen For Update
    strategyUpdateNode = strategyPGG(nodeChosenUpdate)
    !
    payoffUpdate = 0.0

    !!! # of A opponents
    numAnodeUpdatePGG = sum(strategyPGG(networkMain(nodeChosenUpdate)%neighborNode))

    !!! calculate the payoff of the focal game
    if(strategyUpdateNode == 1)then
        payoffUpdate = payoffUpdate + APlayerPayoff(numAnodeUpdatePGG + 1)
    else
        payoffUpdate = payoffUpdate + BPlayerPayoff(numAnodeUpdatePGG + 1)
    end if

    !write(*,*) "nodeChosenUpdate", nodeChosenUpdate
    !pause

    if(payoffCollectType == 0)then
        !!! only Single focal game is considered
        return
    end if

    !!!! Calculating the payoffs of games organized by neighbors
    do indexNeigh = 1, degreeMain
        nodeNeighNow = networkMain(nodeChosenUpdate)%neighborNode(indexNeigh)
        strategyNeighNode = strategyPGG(nodeNeighNow)
        !!! # of A opponents of the Update Node
        numANeighborPGG = strategyNeighNode + sum(strategyPGG(networkMain(nodeNeighNow)%neighborNode)) - strategyUpdateNode
        if(strategyUpdateNode == 1)then
            payoffUpdate = payoffUpdate + APlayerPayoff(numANeighborPGG + 1)
        else
            payoffUpdate = payoffUpdate + BPlayerPayoff(numANeighborPGG + 1)
        end if
    end do
    !!!! accumulated payoff (default)
    !!!! average payoff
    payoffUpdate = payoffUpdate / (degreeMain + 1)

    return
end subroutine PGG_AroundUpdateNode_RegularGraph

max_stars_repo_name: leizhougetbetter/IAD

max_issues_repo_name: leizhougetbetter/IAD

max_forks_repo_name: leizhougetbetter/IAD


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!     calculate the maximal admissible pressure ratio pt2/pt1
!
!     1) assuming M2=1 for adiabatic respectively M2=1/dsqrt(kappa)
!        for isotherm pipe choking
!        M1 is calculated iteratively using a dichotomy scheme
!
!     2)the ratio of the critical pressure ratio
!       Qred_1/Qred_2crit=Pt2/Pt1=D(M1)/D(M2_crit)
!       is computed
!       [D(M)=M*(1+0.5*(kappa-1)*M)**(-0.5*(kappa+1)/(kappa-1))]
!       (general gas equation)
!
!     author: Yannick Muller
!
      subroutine pt2zpt1_crit(pt2,pt1,Tt1,lambda,kappa,r,l,d,
     &     inv,pt2zpt1_c,Qred_crit,crit,icase,M1)
!
      implicit none
!
      logical crit
!
      integer inv,icase,i
!
      real*8 pt2,pt1,lambda,kappa,l,d,M1,pt2zpt1,pt2zpt1_c,
     &     km1,kp1,kp1zk,Tt1,r,xflow_crit,Qred_crit,fmin,
     &     f,fmax,M1_min,M1_max,lld
!
      intent(in) pt2,pt1,Tt1,lambda,kappa,r,l,d,
     &     inv,icase
!
      intent(inout) Qred_crit,crit,pt2zpt1_c
!
      crit=.false.
!
!     useful variables and constants
!
      km1=kappa-1.d0
      kp1=kappa+1.d0
      kp1zk=kp1/kappa
      lld=lambda*l/d
!
!     adiabatic case
!
      if(icase.eq.0) then
!
!        computing M1 using dichotomy method (dividing the interval with the function
!        root iteratively by 2)
!
         i=1
!
         M1_min=0.001d0
         M1_max=1
!
         fmin=(1.d0-M1_min**2)*(kappa*M1_min**2)**(-1)
     &        +0.5d0*kp1zk*log((0.5d0*kp1)*M1_min**2
     &        *(1+0.5d0*km1*M1_min**2)**(-1))-lld
!
         fmax=(1.d0-M1_max**2)*(kappa*M1_max**2)**(-1)
     &        +0.5d0*kp1zk*log((0.5d0*kp1)*M1_max**2
     &        *(1+0.5d0*km1*M1_max**2)**(-1))-lld
         do
            i=i+1
            M1=(M1_min+M1_max)*0.5d0
!
            f=(1.d0-M1**2)*(kappa*M1**2)**(-1)
     &           +0.5d0*kp1zk*log((0.5d0*kp1)*M1**2
     &           *(1+0.5d0*km1*M1**2)**(-1))-lld
!
            if(abs(f).le.1E-6) then
               exit
            endif
            if(i.gt.50) then
               exit
            endif
!
            if(fmin*f.le.0.d0) then
               M1_max=M1
               fmax=f
            else
               M1_min=M1
               fmin=f
            endif
         enddo
!
         Pt2zpt1_c=M1*(0.5d0*kp1)**(0.5*kp1/km1)
     &        *(1+0.5d0*km1*M1**2)**(-0.5d0*kp1/km1)
!
!     isothermal case
!
      elseif (icase.eq.1) then
!
!        computing M1 using dichotomy method for choked conditions M2=1/dsqrt(kappa)
!        (1.d0-kappa*M1**2)/(kappa*M1**2)+log(kappa*M1**2)-lambda*l/d=0
!
         i=1
!
         M1_min=0.1d0
         M1_max=1/dsqrt(kappa)
!
         fmin=(1.d0-kappa*M1_min**2)/(kappa*M1_min**2)
     &        +log(kappa*M1_min**2)-lambda*l/d
!
         fmax=(1.d0-kappa*M1_max**2)/(kappa*M1_max**2)
     &        +log(kappa*M1_max**2)-lambda*l/d
!
         do
            i=i+1
            M1=(M1_min+M1_max)*0.5d0
!
            f=(1.d0-kappa*M1**2)/(kappa*M1**2)
     &           +log(kappa*M1**2)-lambda*l/d
!
            if((abs(f).le.1E-5).or.(i.ge.50)) then
               exit
            endif
!
            if(fmin*f.le.0.d0) then
               M1_max=M1
               fmax=f
            else
               M1_min=M1
               fmin=f
            endif
         enddo
!
!        computing the critical pressure ratio in the isothermal case
!        pt=A*dsqrt(kappa)/(xflow*dsqrt(kappa Tt))*
!           M*(1+0.5d0*(kappa-1)M**2)**(-0.5d0*(kappa+1)/(kappa-1))
!        and forming the pressure ratio between inlet and outlet(choked)
!
c         Pt2zPt1_c=dsqrt(Tt2/Tt1)*M1*dsqrt(kappa)*((1+0.5d0*km1/kappa)
c     &        /(1+0.5d0*km1*M1**2))**(0.5d0*(kappa+1)/km1)
         Pt2zPt1_c=M1*dsqrt(kappa)*((1+0.5d0*km1/kappa)
     &        /(1+0.5d0*km1*M1**2))**(0.5d0*(kappa+1)/km1+0.5d0)
!
      endif
!
      pt2zpt1=pt2/pt1
      if(Pt2zPt1.le.Pt2zPt1_c) then
         crit=.true.
      endif
!
      Qred_crit=M1*dsqrt(kappa/r)
     &     *(1+0.5d0*km1*M1**2)**(-0.5d0*kp1/km1)
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!
! $Id$

program prog_dobsdx

!-----------------------------------------------------------
! Numerically differentiate the observation values for small
! change in input vector component. Used to interpolate
! observation.

implicit none
real(8)        :: deltax = 0.1_8, x
real(8),allocatable :: obsx(:), obsxpd(:), dobsdx(:), locx(:)
character(128) :: file_in ='obsval.dat.bck'
character(128) :: file_out ='obsval.dat'
character(128) :: file_x ='model.in.bck'
integer :: nobs,i,j

nobs = 0
open(11,file=trim(file_in),status='old')
do
  read(11,'(1x)',end=123)
  nobs = nobs + 1
enddo
123 rewind(11)
allocate(obsx(nobs),obsxpd(nobs),dobsdx(nobs),locx(nobs))
open(12,file=trim(file_out),status='old')
do i = 1, nobs
   read(11,*)locx(i),j,obsx(i)
   read(12,*)locx(i),j,obsxpd(i)
   dobsdx(i) = (obsxpd(i) - obsx(i))/deltax
end do
close(11)
rewind(12)
open(11,file=trim(file_x),status='old')
read(11,*)
read(11,*)
read(11,*) x
close(11)
do i = 1, nobs
   write(12,*)locx(i),x,obsx(i),dobsdx(i)
enddo
close(12)
end program prog_dobsdx

! <next few lines under version control, do not edit>
! $URL$
! $Id$
! $Revision$
! $Date$

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


module module_math
  use module_globals
  use module_autodiff
  implicit none

  private
  public IntegralTest, load_integral_test, f, Df

  type :: IntegralTest
     real(8) :: st, en, le, cn, seikai
  end type IntegralTest

contains
  !========================================================================================
  subroutine load_integral_test(test)
    type(IntegralTest),intent(out) :: test

    ! int_[a->b](sin(x))dx
    test%st=0.d0
    test%en=30.d0
    test%seikai=-cos(test%en)+cos(test%st)

    !int_[0.01->2](sin(x)/x)dx
    ! test%st=0.01d0
    ! test%en=1.d0
    ! test%seikai=0.93608312592257190411368860853800106052127274d0

    ! !int_[0->30](sin(x)/x)dx
    ! test%st=1.d0
    ! test%en=30.d0
    ! test%seikai=0.620673469663168096042377995183618508711157d0

    ! int_[-1->1] 1/(1+25*x^2) dx
    test%st=-5.d0
    test%en=5.d0
    test%seikai=(atan(5.d0*test%en)-atan(5.d0*test%st))/5.d0

    ! int_[0,1] bessel_j0(x) dx
    ! test%st=0.d0
    ! test%en=1.d0
    ! test%seikai=0.919730410089760239314421194080619970661d0

    ! test%st=0.d0
    ! test%en=2.d0
    ! test%seikai=(test%en**30-test%st**30)/30.d0

    test%le=test%en-test%st
    test%cn=(test%st+test%en)*0.5d0

  end subroutine load_integral_test

  !========================================================================================
  real(8) function f(x)
    real(8),intent(in) :: x
    !f=sin(x)
    !f=sin(x)/x
    !f=sin(x)/x
    f=1.d0/(1.d0+25.d0*x*x)
    !f=bessel_j0(x)
    !f=x**29
  end function f

  !========================================================================================
  type(Dcomplex8) function Df(x)
    type(Dcomplex8),intent(in) :: x
    !Df=sin(x)
    !Df=sin(x)/x
    !Df=sin(x)/x
    Df=1.d0/(1.d0+25.d0*x*x)
    !Df=bessel_j0(x)
    !Df=x**29
  end function Df

end module module_math


text:
module OFP
  type :: uint4
     integer :: s0, s1, s2, s3
  end type uint4

contains

   function local(A)
      implicit none
      real, target, dimension(:,:,:) :: A
      real, pointer, dimension(:,:,:) :: local

      local => A

   end function local

   function local_section(A, halo)
      implicit none
      real, target,  dimension(:,:,:) :: A
      real, pointer, dimension(:,:,:) :: local_section
      integer :: halo(6), lb(3), ub(3)

      lb = lbound(A)
      ub = ubound(A)

      local_section => A( lb(1)+halo(1) : ub(1)-halo(2), &
                          lb(2)+halo(3) : ub(2)-halo(4), &
                          lb(3)+halo(5) : ub(3)-halo(6)  )

   end function local_section

end module OFP

module LIF

   use, intrinsic :: ISO_C_BINDING

   type, bind(C) :: LIF_params
      real(C_FLOAT) :: Vrest
      real(C_FLOAT) :: Vexc
      real(C_FLOAT) :: Vinh
      real(C_FLOAT) :: VinhB

      real(C_FLOAT) :: tau
      real(C_FLOAT) :: tauE
      real(C_FLOAT) :: tauI
      real(C_FLOAT) :: tauIB

      real(C_FLOAT) :: VthRest
      real(C_FLOAT) :: tauVth
      real(C_FLOAT) :: deltaVth

      real(C_FLOAT) :: noiseFreqE
      real(C_FLOAT) :: noiseAmpE
      real(C_FLOAT) :: noiseFreqI
      real(C_FLOAT) :: noiseAmpI
      real(C_FLOAT) :: noiseFreqIB
      real(C_FLOAT) :: noiseAmpIB
   end type LIF_params

   interface
      pure elemental function cl_random_state(state)
         use OFP
         type(uint4), intent(in) :: state
         type(uint4) :: cl_random_state
      end function
      pure elemental real function cl_random_prob(state)
         use OFP
         type(uint4), intent(in) :: state
      end function
   end interface

contains

!$OFP CL_KERNEL :: update_state
subroutine update_state(time, dt, params,         &
                        nf, nx, ny, nb, rnd,      &
                        V, Vth, G_E, G_I, G_IB,   &
                        phiExc, phiInh, phiInhB,  &
                        activity)
   use OFP
   implicit none
   real, intent(in) :: time, dt
   type(LIF_Params), intent(in) :: params
   integer, intent(in) :: nf, nx, ny, nb
   type(uint4), intent(inout) :: rnd(nf,nx,ny)
   real, intent(inout), dimension(nf,nx,ny) :: V, Vth, G_E, G_I, G_IB
   real, intent(inout), dimension(nf,nx,ny) :: phiExc, phiInh, phiInhB
   real, intent(inout), dimension(nf,nx+2*nb,ny+2*nb) :: activity

   ! params
   real :: tau, tauE, tauI, tauIB, Vrest, VthRest, Vexc, Vinh, VinhB, tauVth, deltaVth

   real, parameter :: GMAX = 10.0
   real :: dt_sec

   ! local param variables
   real, pointer, dimension(:,:,:) :: l_activ
   !$OFP LOCAL :: V, Vth, G_E, G_I, G_IB, phiExc, phiInh, phiInhB, l_activ

   ! temporary arrays
   real, dimension(nf,nx,ny) :: tauInf, VmemInf

   !
   ! start of LIF2_update_exact_linear
   !

   ! define local param variables
   !
   tau   = params%tau
   tauE  = params%tauE
   tauI  = params%tauI
   tauIB = params%tauIB

   Vrest = params%Vrest
   Vexc  = params%Vexc
   Vinh  = params%Vinh
   VinhB = params%VinhB

   tauVth   = params%tauVth
   VthRest  = params%VthRest
   deltaVth = params%deltaVth

   ! add noise
   !
   dt_sec = .001 * dt   ! convert to seconds

   where (cl_random_prob(rnd) < dt_sec*params%noiseFreqE)
      rnd = cl_random_state(rnd)
      phiExc = phiExc + params%noiseAmpE*cl_random_prob(rnd)
      rnd = cl_random_state(rnd)
   end where

   where (cl_random_prob(rnd) < dt_sec*params%noiseFreqI)
      rnd = cl_random_state(rnd)
      phiInh = phiInh + params%noiseAmpI*cl_random_prob(rnd)
      rnd = cl_random_state(rnd)
   end where

   where (cl_random_prob(rnd) < dt_sec*params%noiseFreqIB)
      rnd = cl_random_state(rnd)
      phiInhB = phiInhB + params%noiseAmpIB*cl_random_prob(rnd)
      rnd = cl_random_state(rnd)
   end where

   G_E  = phiExc  + G_E *exp(-dt/tauE)
   G_I  = phiInh  + G_I *exp(-dt/tauI )
   G_IB = phiInhB + G_IB*exp(-dt/tauIB)

   tauInf  = (dt/tau) * (1.0 + G_E + G_I + G_IB)
   VmemInf = (Vrest + G_E*Vexc + G_I*Vinh + G_IB*VinhB) / (1.0 + G_E + G_I + G_IB)

   V = VmemInf + (V - VmemInf)*exp(-tauInf)

   !
   ! start of LIF2_update_finish
   !

   phiExc  = 0.0
   phiInh  = 0.0
   phiInhB = 0.0

   Vth = VthRest + (Vth - VthRest)*exp(-dt/tauVth)

   !
   ! start of update_f
   !

   where (G_E  > GMAX) G_E  = GMAX
   where (G_I  > GMAX) G_I  = GMAX
   where (G_IB > GMAX) G_IB = GMAX

   l_activ => local_section(activity, [0,0,nb,nb,nb,nb])

   where (V > Vth)
      l_activ = 1.0
      V = Vrest             ! reset cells that fired
      Vth = Vth + deltaVth  !
      G_IB = G_IB + 1.0     ! add hyperpolarizing current
   elsewhere
      l_activ = 0.0
   end where

   !
   ! These actions must be done outside of kernel
   !    1. set activity to 0 in boundary (if needed)
   !    2. update active indices
   !

end subroutine update_state

end module LIF

max_stars_repo_name: PetaVision/ObsoletePV

max_issues_repo_name: PetaVision/ObsoletePV

max_forks_repo_name: PetaVision/ObsoletePV


text:
program laplace
  integer,parameter :: N1=64, N2=64
  real(8),parameter :: PI = 3.141592653589793238463

  integer :: niter=100
  real(8) :: u(N1,N2),uu(N1,N2)
  real(8) :: value = 0.0

  do j=1,N2
     do i=1,N1
        u(i,j)=0.0
        uu(i,j)=0.0
     end do
  end do

  do j=2,N2-1
     do i=2,N1-1
        u(i,j)=sin(dble(i-1)/N1*PI)+cos(dble(j-1)/N2*PI)
     end do
  end do

  do k=1,niter

     do j=2,N2-1
        do i=2,N1-1
           uu(i,j)=u(i,j)
        end do
     end do

     do j=2,N2-1
        do i=2,N1-1
           u(i,j)=(uu(i-1,j) + uu(i+1,j) + uu(i,j-1) + uu(i,j+1))/4.0
        end do
     end do

  enddo

  do j=2,N2-1
     do i=2,N1-1
        value = value + dabs(uu(i,j) -u(i,j))
     end do
  end do

  print *, 'Verification =', value

end program laplace

max_stars_repo_name: XcalableMP/lecture

max_issues_repo_name: XcalableMP/lecture

max_forks_repo_name: XcalableMP/lecture


text:
c
c ----------------------------------------------------------
c
      subroutine prepc(level,nvar)
c
      implicit double precision (a-h,o-z)
      logical            ovrlap
      logical            graf
      parameter  (maxgr = 192, maxlv=12)
      common /nodal/ rnode(12,maxgr),node(17,maxgr),lstart(maxlv),
     *  newstl(maxlv),
     *  listsp(maxlv),intrat(maxlv),tol,bzonex,bzoney,mstart,ndfree,
     *  lfine,iorder,mxnest,kcheck,lwidth,
     *  maxcl, graf, lhead
      include "calloc.i"
c
c this routine called because regridding just changed the fine grids.
c modify coarse grid boundary lists to store fluxes in appropriate
c fine grids.
c assume new fine grids have node(15) initialized to point to null
c
c  first compute max. possible number of list cells. allocate
c  initially so that one pass through is enough. at end, release
c  unused space
c
      maxsp  = 0
      mkid   = lstart(level+1)
 10   if (mkid .eq. 0) go to 20
         ikeep  = (node(5,mkid)-1)/intrat(level)
         jkeep  = (node(6,mkid)-1)/intrat(level)
         maxsp  = maxsp + 2*(ikeep+jkeep)
      mkid = node(10,mkid)
      go to 10
 20   listsp(level) = maxsp
      if (maxsp .eq. 0) go to 99
c
      mpar = lstart(level)
 30   if (mpar .eq. 0) go to 99
c
       ispot   = 0
       hxpar   = rnode(9, mpar)
       hypar   = rnode(10,mpar)
       mipar   = node(5,mpar)
       mjpar   = node(6,mpar)
       locbc   = igetsp(5*maxsp)
       node(15,mpar) = locbc
c
c  initialize list space to 0 (use 0 terminator to indicate end of bc
c  list processing).
c
       call ziplst(alloc(locbc),maxsp)
c
       mkid = lstart(level+1)
 40    if (mkid .eq. 0) go to 60
       if (.not. ovrlap(mpar,mkid)) go to 50
          mikid  = node(5,mkid)
          mjkid  = node(6,mkid)
          hxkid  = rnode(9,mkid)
          hykid  = rnode(10,mkid)
c
          ist  = (rnode(1,mkid)-rnode(1,mpar))/hxpar + 1.01
          ist  = max0(ist,1)
          xist = rnode(1,mpar)+dfloat(ist-1)*hxpar
          kst  = (xist - rnode(1,mkid))/hxkid + 1.01
          iend = (rnode(7,mkid)-rnode(1,mpar))/hxpar+1.01
          iend = min0(iend,mipar)
          jst  = (rnode(2,mkid)-rnode(2,mpar))/hypar + 1.01
          jst  = max0(jst,1)
          xjst = rnode(2,mpar) + dfloat(jst-1)*hypar
          lst  = (xjst - rnode(2,mkid))/hykid + 1.01
          jend = (rnode(4,mkid)-rnode(2,mpar))/hypar+1.01
          jend = min0(jend,mjpar)
c
          call setuse(alloc(locbc),maxsp,ist,iend,jst,jend,ispot,mkid,
     1          kst,lst,intrat(level),mikid,mjkid,mipar,mjpar)
c
 50     mkid = node(10,mkid)
        go to 40
c
c  done with subgrid cycle. if no cells would need fixing, all done
c  else cycle through again to set up list with info. for bc processing
c
 60     continue
c  for now, leave unused space allocated to the grid. alternative is
c  to return (maxsp-ispot) amount starting at loc. node(15,mpar)+ispot.
c
       mpar = node(10,mpar)
       go to 30
c
 99    return
       end