text:
!***************************************************************************************************
!>Main program (SPKMC algorithm)
!Function: it's used to simulate the process of radiation damage accumulation and subsequent annealing in metals
!          with the Synchronous Sublattice Kinetic Monte Carlo (SPKMC) algorithm.
!Date: 2019/05
!Author: CDD (chendandan@xs.ustb.edu.cn)
!Amendant Record: original code
!***************************************************************************************************
program MISASCD
	use mod_constants
	use mod_structures			!<variable classes for MISASCD
	use mod_globalVariables		!<module containing all global variables
	use mod_randdp			!<module for double precision random number generation
	implicit none
	include 'mpif.h'

	type(defect), pointer :: defectCurrent									!<used to find defects
	type(defectUpdateTracker), pointer :: defectUpdate, defectUpdateCurrent !<used to update reactions
	type(reaction), pointer :: reactionCurrent								!<used to find reactions
	type(reaction), pointer :: reactionChoiceList, reactionChoiceCurrent	!<used to create a list of chosen reactions (for one KMC domain per volume element case)
	type(reaction), pointer :: reactionTemp
	type(cascade), pointer :: CascadeCurrent								!<used to find defects/reactions in fine mesh
	double precision :: tau
	integer :: status(MPI_STATUS_SIZE), sim, outputCounter, nullSteps
	integer :: cascadeCell, i, CascadeCount, TotalCascades
	real :: time1, time2, time3
	double precision :: runTime1, runTime2
	double precision :: commTime1, commTime2
	logical :: releaseToggle, impCascadeToggle
	integer :: timeCounter
	character(14) :: filename1, filename2, filename3, filename4
	integer :: simStatus	!<1: 'irradiation', 0: 'anneal'
	double precision, external :: GenerateTimestep, TotalRateCheck

	!test
	type(defect), pointer :: defectTest
	type(reaction), pointer :: reactionTest
	integer :: impCount

	interface
		subroutine chooseImplantReaction(reactionCurrent, CascadeCurrent)
			use mod_structures
			implicit none
			type(reaction), pointer, intent(inout) :: reactionCurrent
			type(cascade), pointer, intent(inout) :: CascadeCurrent
		end subroutine

		subroutine chooseReaction(reactionCurrent, CascadeCurrent)
			use mod_structures
			implicit none
			type(reaction), pointer, intent(inout) :: reactionCurrent
			type(cascade), pointer, intent(inout) :: CascadeCurrent
		end subroutine

		subroutine addCascadeExplicit(reactionCurrent, CascadeCurrent)
			use mod_structures
			implicit none
			type(reaction), pointer, intent(inout) :: reactionCurrent
			type(cascade), pointer, intent(inout) :: cascadeCurrent
		end subroutine

		subroutine updateDefectList(reactionCurrent, defectUpdateCurrent, CascadeCurrent)
			use mod_structures
			implicit none
			type(reaction), pointer, intent(in) :: reactionCurrent
			type(DefectUpdateTracker), pointer, intent(inout) :: defectUpdateCurrent
			type(cascade), pointer, intent(inout) :: CascadeCurrent
		end subroutine

		subroutine updateReactionList(defectUpdate)
			use mod_structures
			implicit none
			type(DefectUpdateTracker), pointer, intent(inout) :: defectUpdate
		end subroutine

		subroutine releaseFineMeshDefects(CascadeCurrent)
			use mod_structures
			implicit none
			type(cascade), pointer, intent(inout) :: CascadeCurrent
		end subroutine

		double precision function totalRateCascade(CascadeCurrent)
			use mod_structures
			implicit none
			type(Cascade), pointer, intent(in) :: CascadeCurrent
		end function
	end interface

	call cpu_time(time1)
	timeCounter = 1
	impCount = 1
	!***********************************************************************
	!<Initialize MPI interface
	!***********************************************************************
	!periods = .true.
	call MPI_INIT(ierr)
	call MPI_COMM_SIZE(MPI_COMM_WORLD, myProc%numtasks, ierr)		!read number of processors
	call MPI_DIMS_CREATE(myProc%numtasks,3,dims, ierr)
	call MPI_CART_CREATE(MPI_COMM_WORLD,3,dims,periods,.true.,comm,ierr)
	call MPI_CART_GET(comm,3,dims,periods,myProc%coords,ierr)
	call MPI_CART_RANK(comm,myProc%coords,myProc%taskid,ierr)

	call MPI_CART_SHIFT(comm,0,1,myProc%procNeighbor(2),myProc%procNeighbor(1),ierr) !x
	call MPI_CART_SHIFT(comm,1,1,myProc%procNeighbor(4),myProc%procNeighbor(3),ierr) !y
	call MPI_CART_SHIFT(comm,2,1,myProc%procNeighbor(6),myProc%procNeighbor(5),ierr) !z

	if(myProc%taskid==MASTER) then
		write(*,*) 'proc division', dims
	end if

	runTime1=MPI_WTIME()
	commTimeSum=0.0

	!***********************************************************************
	!Initialize input parameters and meshes
	!***********************************************************************
	call ReadInputs()
	call initialMesh()			!open Mesh_xx.txt file and carry out parallel mesh initialization routine
	!<!<Initialize cascadeReactionLimit
	if(implantType=='Cascade' .AND. meshingType=='adaptive') then
		cascadeReactionLimit = 10d0*myMesh(1)%volume*dpaRate/(numDisplacedAtoms*atomSize)	!=10 times the cascade implantation rate
	else
		cascadeReactionLimit = 0d0
	end if

	!<Create fine mesh connectivity
	allocate(cascadeConnectivity(6, numCellsCascade))	!< numCellsCascade=numxFine*numyFine*numzFine
	call  createCascadeConnect()

	!***********************************************************************
	!Multiple simulations begin.
	!***********************************************************************
	do sim=1, numSims

		temperature=tempStore	!< Temperature read in (K)
		!*******************************************************
		!<Initialize outputfile
		!*******************************************************
		if(myProc%taskid==MASTER) then
			if(totdatToggle=='yes') then
				filename1(1:7)='totdat_'
				write(unit=filename1(8:10), fmt='(I3.3)') sim
				filename1(11:14)='.out'
				open(TOTFILE, file=filename1, action='write', status='Unknown')
			end if

			if(defectToggle=='yes') then
				filename2(1:7)='defect_'
				write(unit=filename2(8:10), fmt='(I3.3)') sim
				filename2(11:14)='.out'
				open(DEFFILE , file=filename2, action='write', status='Unknown')
			end if

			if(stadatToggle=='yes') then
				filename3(1:7)='stadat_'
				write(unit=filename3(8:10), fmt='(I3.3)') sim
				filename3(11:14)='.out'
				open(STAFILE , file=filename3, action='write', status='Unknown')
			end if

			if(xyzdatToggle=='yes') then
				filename4(1:7)='xyzdat_'
				write(unit=filename4(8:10), fmt='(I3.3)') sim
				filename4(11:14)='.out'
				open(XYZFILE , file=filename4, action='write', status='Unknown')
			end if
		end if

		call initializeRandomSeeds()		!<set unique random number seeds in each processor
		allocate(defectList(numCells))		!<Create list of defects - array
		allocate(reactionList(numCells))	!<Create list of reactions - array
		allocate(totalRateVol(numCells))	!<Create array of total rates in each volume element
		call cpu_time(time3)
		call initializeVIdefect()			!<initialize point defects in the whole system
		call initializeDefectList()			!<initialize defects within myMesh
		call initializeBoundaryDefectList()	!<initialize defects on boundary of myMesh (in other procs)

		call cpu_time(time2)
		if(myProc%taskid==MASTER) then
			write(*,*) 'initialize VI time', time2-time3
		end if
		call initializeReactionList()		!<initialize reactions within myMesh
		call initializeTotalRate()			!<initialize totalRate and maxRate using reactionList(:)

		if(myProc%taskid==MASTER) then
			write(*,*) 'CuEverMesh', numCuCell
			write(*,*) 'Initial number of V', initialNumV
			write(*,*) 'Initial number of SIA', initialNumI
			write(*,*) 'Equilibrium concentration of V', ceqV
			write(*,*) 'Equilibrium concentration of SIA', ceqI
			write(*,*) 'totalRate', totalRate
		end if

		!***********************************************************************
		!<Restart from a reset file.
		!<If we are restarting from a reset file, then we have to start wtih a nonzero dpa and
		!<at a nonzero time. All of the 'old' implant events are tracked in the master processor.
		!***********************************************************************
		numImpAnn(1)=0     	!<Postprocessing: number of Frenkel pairs / cascades (local)
		numImpAnn(2)=0     	!<Postprocessing: number of annihilation reactions carried out
		elapsedTime=0d0		!<The simulation time that has passed

		numTrapSIA=0		!<Postprocessing: number of SIAs trapped on grain boundary
		numTrapV=0			!<Postprocessing: number of vacancies trapped on grain boundary
		numEmitSIA=0		!<Postprocessing: number of SIAs emitted from grain boundary
		numEmitV=0			!<Postprocessing: number of vacancies emitted from grain boundary

		!<Initialize totalTime
		if(totalDPA > 0d0 .AND. dpaRate > 0d0) then
			totalTime=totalDPA/dpaRate	!simulation time
		else
			totalTime = agingTime
		end if

		!<Initialize Counters
		step=0
		nullSteps=0		!Record the number of steps in which an empty event was selected
		TotalCascades=0
		outputCounter=0
		nullify(ActiveCascades)		! pointer ActiveCascades=NULL
		annealIdentify=.FALSE.		!(.TRUE. if in annealing phase, .FALSE. otherwise) used to determine how to reset reaction rates (should we include implantation or not)
		!<Initizlize maxRate and tau. rateTau(1)=maxRate, rateTau(2)=tau
		rateTau=0d0

		!*************************************************************************************************************
		!*************************************************************************************************************
		!**																		                                    **
		!** 			                               Damage Accumulation Loop								        **
		!** During this process, Frenkerl Pairs or cascades are continuously implanted into system. 				**
		!**																											**
		!**************************************************************************************************************
		!**************************************************************************************************************
		!<KMC LOOP:
		!<1. Calculate total reaction rate and time increment.
		!<2, Choose a reaction
		!<3. Carry out reaction (update defect)->note all relevant defects that have been changed
		!<4. Update reaction lists->use defects deleted and defects added (reactants and products) to update relevant reactions
		!<	(first, delete all reactions associated with the REACTANTS AND PRODUCTS in the cell and all DIFFUSION reactions associated with
		!<	reactants in neighboring cells. Then add all reaction rates associated with REACTANTS AND PRODUCTS in the cell and diffusion
		!<	reactions in neighboring cells.)
		!<5. Repeat 1-4.
		do while(elapsedTime < totalTime)
		!do while(step < 1)

			step=step+1
			!<Logical variable tells us whether cascade communication step needs to be carried out
			!<(0=no cascade, nonzero=number of volume element where cascade event has happened)
			cascadeCell=0
			impCascadeToggle= .FALSE.

			!<find the maximum totalRate in all processors
			commTime1=MPI_WTIME()
			call MPI_REDUCE(totalRate,maxRate,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,comm,ierr)
			commTime2=MPI_WTIME()
			commTimeSum=commTimeSum+commTime2-commTime1

			if(myProc%taskid==MASTER) then
				rateTau(1)=maxRate
			end if

			!<Generate timestep in the master processor and send it to all other processors
			!<If implantScheme=='explicit', cascade implantation needs to be carried out explicitly
			if(implantScheme=='explicit') then
				if(elapsedTime >= numImpAnn(1)*(numDisplacedAtoms*atomSize)/(totalVolume*dpaRate)) then
					!<Do not generate a timestep in this case; this is an explicit (zero-time) reaction
					if(myProc%taskid==MASTER) then
						rateTau(2)=0d0
					end if
				else
					if(myProc%taskid==MASTER) then
						rateTau(2)=GenerateTimestep()
					end if
				end if
			else if(implantScheme=='MonteCarlo') then
				if(myProc%taskid==MASTER) then
					rateTau(2)=GenerateTimestep()
				end if
			end if

			commTime1=MPI_WTIME()
			call MPI_BCAST(rateTau, 2, MPI_DOUBLE_PRECISION, MASTER, comm,ierr)
			commTime2=MPI_WTIME()
			commTimeSum=commTimeSum+commTime2-commTime1

			maxRate=rateTau(1)	!<update maxRate of the processor
			tau=rateTau(2)		!<update time step of the processor

			!***********************************************************************************************
			!Choose from reactions in reactionList(:) (local to this processor). Null events possible.
			!If explicit implant scheme has been chosen, implant cascades when elapsed time has passed
			!a threshold, and do not iterate timestep when doing so.
			!***********************************************************************************************
			!<Initialize defectUpdate. defectUpdate is used to store updated defects
			allocate(defectUpdate)
			allocate(defectUpdate%defectType(SPECIES))
			do i=1,SPECIES
				defectUpdate%defectType(i)=0
			end do
			nullify(defectUpdate%next)
			defectUpdateCurrent=>defectUpdate

			!<choose a reaction in one volume element
			!<If implantScheme=='explicit', cascade implantation needs to be carried out explicitly
			if(implantScheme=='explicit') then
				if(elapsedTime >= numImpAnn(1)*(numDisplacedAtoms*atomSize)/(totalVolume*dpaRate)) then
				!if(elapsedTime >= 2.0d0*impCount) then
				!	impCount=impCount+1

					!<choose a cell in the peocessor to implant cascade
					call addCascadeExplicit(reactionCurrent, CascadeCurrent)
				else
					call chooseReaction(reactionCurrent, CascadeCurrent)
				end if
			else if(implantScheme=='MonteCarlo') then
				if(test3 == 'yes') then
					call initializeOneCascade()
					call chooseImplantReaction(reactionCurrent, CascadeCurrent)
				else
					call chooseReaction(reactionCurrent, CascadeCurrent)
				end if
			else
				write(*,*) 'error choosing reaction main program'
			end if

			!***********************************************************************************************
			!Update defects according to reactions chosen. Communicate defects that have changed on
			!boundaries of other processors and create a list of defects whose reaction rates must be update
			!in this processor due to defects updated in this processor and in neighboring processors.
			!***********************************************************************************************
			if(.NOT. associated(reactionCurrent)) then
				nullSteps=nullSteps+1
			end if

			call updateDefectList(reactionCurrent, defectUpdateCurrent, CascadeCurrent)

			!<If a cascade is chosen, update reaction rates for all defects remaining in coarse mesh element
			if(associated(reactionCurrent)) then
				if(reactionCurrent%numReactants==-10) then
					!Resets reaction list in cell and neighbors (on same processor)
					call resetReactionListSingleCell(reactionCurrent%cellNumber(1))
					!Variable tells us whether cascade communication step needs to be carried out
					!and if so in what volume element in the coarse mesh
					cascadeCell=reactionCurrent%cellNumber(1)
					impCascadeToggle=.TRUE.
				end if
			end if

			!<Update elapsed time based on tau, generated timestep. If cascade implant chosen in explicit scheme, tau=0d0
			elapsedTime=elapsedTime+tau

			call updateReactionList(defectUpdate)

			!If we have chosen an event inside a fine mesh, we check the total reaction rate within that
			!fine mesh. If the total rate is less than a set value, we assume the cascade is annealed and
			!release the defects into the coarse mesh.
			releaseToggle = .FALSE.
			if(associated(CascadeCurrent)) then
				if(totalRateCascade(CascadeCurrent) < cascadeReactionLimit) then
					!Record the coarse mesh cell number of cascade (it is destroyed in releaseFineMeshDefects)
					!Used to reset reaction list and to tell cascadeUpdateStep whether a cascade event has occurred
					cascadeCell=CascadeCurrent%cellNumber
					releaseToggle=.TRUE.
					!Release cascade defects into coarse mesh cell and reset the reaction list within that cell
					call releaseFineMeshDefects(CascadeCurrent)
					call resetReactionListSingleCell(cascadeCell)
				end if
			end if

			!Cascade communication step:
			!Tell neighbors whether a cascade has occurred in a cell that is a boundary of a neighbor.
			!If so, update boundary mesh (send defects to boundary mesh) and update all diffusion reaction rates.
			if(implantType=='Cascade') then
				call cascadeUpdateStep(releaseToggle, cascadeCell)
			end if

			!***test
			!if(cascadeCell /= 0) then
			!	write(*,*) 'myProc%taskid', myProc%taskid, 'cell', cascadeCell
			!	nullify(reactionTest)
			!	reactionTest=>reactionList(cascadeCell)
			!	do while(associated(reactionTest))
			!		write(*,*) 'numReactants', reactionTest%numReactants, 'numProducts', reactionTest%numProducts
			!		write(*,*) 'reactants', reactionTest%reactants
			!		write(*,*) 'products', reactionTest%products
			!		write(*,*) 'cellNumber', reactionTest%cellNumber
			!		write(*,*) 'proc', reactionTest%taskid
			!		write(*,*) 'reactionRate', reactionTest%reactionRate
			!		reactionTest=>reactionTest%next
			!	end do
			!end if

			!****************************************************************************************
			!Every time a cascade is added, we reset the total rate and check how wrong it has become
			!!'Cascade' implant type, check total rate when implantation event is choosen
			!!'FrenkelPair' implant type, check total rate every 1000 steps
			!*****************************************************************************************
			if(associated(reactionCurrent)) then
				if(implantType=='Cascade') then
					if(impCascadeToggle .EQV. .TRUE.) then	!implanted a cascade
						totalRate=TotalRateCheck()
					end if
				else if(implantType=='FrenkelPair') then
					if(mod(step,1000)==0) then
						totalRate=TotalRateCheck()
					end if
				end if
			end if

			!******************************************
			! Optional: count how many cascades are present at step i and compile to find avg. number of cascades present per step
			!******************************************
			TotalCascades=TotalCascades+CascadeCount()

			!********************************************************************************
			! Output according to outputCounter
			!********************************************************************************
			if(elapsedTime >= totalTime/2.0d6*(2.0d0)**(outputCounter)) then
			!if(elapsedTime >= totalTime/50d0*(outputCounter)) then
			!if(mod(step,100000)==0) then
				simStatus=1	!<1: 'irradiation', 0: 'anneal'
				commTime1=MPI_WTIME()
				call MPI_REDUCE(numImpAnn,totalImpAnn, 2, MPI_INTEGER, MPI_SUM, 0,comm, ierr)
				commTime2=MPI_WTIME()
				commTimeSum=commTimeSum+commTime2-commTime1
				if(xyzdatToggle=='yes') call outputDefectsXYZ()	!write xyzdat.out
				if(totdatToggle=='yes' .OR. defectToggle=='yes' .OR. stadatToggle=='yes') then
					call outputDefectsTotal(simStatus)	!write totdat.out, defect.out, stadat.out
				end if
				call cpu_time(time2)
				if(myProc%taskid==MASTER) then
					DPA=dble(totalImpAnn(1))/(systemVol/(numDisplacedAtoms*atomSize))
					write(*,*)
					write(*,*) 'time', elapsedTime, 'DPA', DPA, 'steps', step, 'AverageTimeStep', elapsedTime/dble(step)
					if(implantType=='FrenkelPair') then
						write(*,*) 'FrenkelPairs', totalImpAnn(1), 'computationTime', time2-time1
					else if(implantType=='Cascade')	then
						write(*,*) 'Cascades', totalImpAnn(1), 'computationTime', time2-time1
					else	!Thermal aging
						write(*,*) 'noImplantation', totalImpAnn(1), 'computationTime', time2-time1
					end if
					write(*,*)
				end if

				outputCounter=outputCounter+1
				call MPI_BARRIER(comm,ierr)
			end if

			!******************************************
			!used for test 3
			!******************************************
			if(test3 == 'yes') then
				exit
			end if

			!do while(elapsedTime >= 1.0d-4/dpaRate .AND. timeCounter==1)
			!	if(myProc%taskid==MASTER) then
			!		call cpu_time(time2)
			!		write(*,*) '----1d-4 computationTime: ', time2-time1
			!	end if
			!	timeCounter=timeCounter+1
			!	exit
			!end do
		end do	!end of do while(elapsedTime < totalTime)

		!***********************************************************************
		!Output defects at the end of the implantation loop
		!***********************************************************************
		simStatus=1	!<1: 'irradiation', 0: 'anneal'
		commTime1=MPI_WTIME()
		call MPI_REDUCE(numImpAnn, totalImpAnn, 2, MPI_INTEGER, MPI_SUM,0, comm, ierr)
		commTime2=MPI_WTIME()
		commTimeSum=commTimeSum+commTime2-commTime1
		if(xyzdatToggle=='yes') call outputDefectsXYZ()	!write xyzdat.out
		if(totdatToggle=='yes' .OR. defectToggle=='yes' .OR. stadatToggle=='yes') then
			call outputDefectsTotal(simStatus)	!write totdat.out, defect.out, stadat.out
		end if
		call cpu_time(time2)
		if(myProc%taskid==MASTER) then
			DPA=dble(totalImpAnn(1))/(systemVol/(numDisplacedAtoms*atomSize))
			write(*,*)
			write(*,*) 'Final  step'
			write(*,*) 'time', elapsedTime, 'DPA', DPA, 'steps', step, 'AverageTimeStep', elapsedTime/dble(step)
			if(implantType=='FrenkelPair') then
				write(*,*) 'FrenkelPairs', totalImpAnn(1), 'computationTime', time2-time1
			else if(implantType=='Cascade')	then
				write(*,*) 'Cascades', totalImpAnn(1), 'computationTime', time2-time1
			else
				write(*,*) 'noImplantation', totalImpAnn(1), 'computationTime', time2-time1
			end if
			write(*,*)
		end if

		!*************************************************************************************************************
		!*************************************************************************************************************
		!**																		                                   	**
		!** 							              Annealing Loop								     			**
		!**																		   									**
		!*************************************************************************************************************
		!*************************************************************************************************************

		!***************************************************************************
		!Annealing loop: change the temperature to annealTemp and change all
		!defect production rates to 0 (dpa rate, He implant rate)
		!Then carry out simulation normally until elapsedTime=totalTime+annealTime
		!***************************************************************************
		totalTime=0d0
		if(annealTime > 0d0) then
			call annealInitialization()	!Sets up simulation for annealing (change temp, etc)
			annealIdentify=.TRUE.
			outputCounter=0		!Used to output once for each annealing step
			elapsedTime=0d0
			totalTime=annealTime
			step=0
			nullSteps=0
		end if

		totalRate=TotalRateCheck()

		!if annealTime = 0, the following "do" does not execute
		!!do while(elapsedTime < totalTime+dble(annealIter)*annealTime/dble(annealSteps))
		do while(elapsedTime < totalTime)

			step=step+1
			!Logical variable tells us whether cascade communication step needs to be carried out
			!(0=no cascade, nonzero=number of volume element where cascade event has happened)
			cascadeCell=0

			call MPI_REDUCE(totalRate,maxRate,1,MPI_DOUBLE_PRECISION,MPI_MAX,0,comm,ierr)

			if(myProc%taskid==MASTER) then
				rateTau(1)=maxRate
			end if

			if(myProc%taskid==MASTER) then
				rateTau(2)=GenerateTimestep()
				!if(elapsedTime-totalTime+tau > annealTime/dble(annealSteps)*outputCounter) then
				!we have taken a timestep that moves us past this annealing step
				!	rateTau(2)=annealTime/dble(annealSteps)*outputCounter-(elapsedTime-totalTime)
				!end if
			end if

			call MPI_BCAST(rateTau, 2, MPI_DOUBLE_PRECISION, MASTER, comm,ierr)

			maxRate=rateTau(1)	!update maxRate
			tau=rateTau(2)		!update time step

			!***********************************************************************************************
			!Choose from reactions in in reactionList(:) (local to this processor). Null events possible.
			!***********************************************************************************************
			allocate(defectUpdate)
			allocate(defectUpdate%defectType(SPECIES))
			do i=1,SPECIES
				defectUpdate%defectType(i)=0
			end do
			nullify(defectUpdate%next)
			defectUpdateCurrent=>defectUpdate

			!choose a reaction in one volume element
			call chooseReaction(reactionCurrent, CascadeCurrent)

			!************
			! Optional: count how many steps are null events
			!************
			if(.NOT. associated(reactionCurrent)) then
				nullSteps=nullSteps+1
			end if

			call updateDefectList(reactionCurrent, defectUpdateCurrent, CascadeCurrent)

			if(associated(reactionCurrent)) then
				if(reactionCurrent%numReactants==-10) then
					write(*,*) 'Error chose casacde implantation in annealing phase'
				end if
			end if

			elapsedTime=elapsedTime+tau
			call updateReactionList(defectUpdate)

			releaseToggle=.FALSE.
			if(associated(CascadeCurrent)) then
				if(totalRateCascade(CascadeCurrent) < cascadeReactionLimit) then

					cascadeCell=CascadeCurrent%cellNumber
					releaseToggle=.TRUE.
					call releaseFineMeshDefects(CascadeCurrent)
					call resetReactionListSingleCell(cascadeCell)
					write(*,*) 'taskid', myProc%taskid, 'step', step, 'anneal step', annealIter, &
							'totalRateCascade',totalRateCascade(CascadeCurrent),'release fine meshes'
				end if
			end if

			call cascadeUpdateStep(releaseToggle,cascadeCell)

			!Update the totalRate in order to avoid any truncation error every 1000 steps (this value can be modified)
			if(mod(step,1000)==0) then
				totalRate=TotalRateCheck()
			end if

			!******************************************
			! Optional: count how many cascades are present at step i and compile to find avg. number of cascades present per step
			!******************************************
			TotalCascades=TotalCascades+CascadeCount()

			!******************************************
			! Output
			!******************************************
			!!if((elapsedTime-totalTime) >= annealTime/dble(annealSteps)*outputCounter) then
			if(elapsedTime >= totalTime/2.0d6*(2.0d0)**(outputCounter)) then
				simStatus=0	!<1: 'irradiation', 0: 'anneal'
				!call MPI_REDUCE(numImpAnn, totalImpAnn, 2, MPI_INTEGER, MPI_SUM,0,comm, ierr)
				if(xyzdatToggle=='yes') call outputDefectsXYZ()	!write xyzdat.out
				if(totdatToggle=='yes' .OR. defectToggle=='yes' .OR. stadatToggle=='yes') then
					call outputDefectsTotal(simStatus)	!write totdat.out, defect.out, stadat.out
				end if
				call cpu_time(time2)
				if(myProc%taskid==MASTER) then
					!DPA=dble(totalImpAnn(1))/(systemVol/(numDisplacedAtoms*atomSize))
					write(*,*)
					write(*,*) 'Annealing process'
					write(*,*) 'Anneal time', elapsedTime, 'steps',step, 'AverageTimeStep', elapsedTime/dble(step)
					write(*,*) 'computationTime', time2-time1
					write(*,*)
				end if

				outputCounter=outputCounter+1
				!call MPI_BARRIER(comm，ierr)
			end if
		end do	!Anneal time loop

		!***********************************************************************
		!Output final defect state
		!***********************************************************************
		if(annealTime > 0d0) then
			simStatus=0	!<1: 'irradiation', 0: 'anneal'
			!call MPI_REDUCE(numImpAnn,totalImpAnn, 2, MPI_INTEGER, MPI_SUM, 0, comm, ierr)
			if(xyzdatToggle=='yes') call outputDefectsXYZ()	!write xyzdat.out
			if(totdatToggle=='yes' .OR. defectToggle=='yes' .OR. stadatToggle=='yes') then
				call outputDefectsTotal(simStatus)	!write totdat.out, defect.out, stadat.out
			end if
			call cpu_time(time2)
			!write(*,*) 'Fraction null steps', dble(nullSteps)/dble(step), 'Proc', myProc%taskid
			if(myProc%taskid==MASTER) then
				!DPA=dble(totalImpAnn(1))/(systemVol/(numDisplacedAtoms*atomSize))
				write(*,*)
				write(*,*) 'Annealing process'
				write(*,*) 'Final Defect State'
				write(*,*) 'Anneal time', elapsedTime, 'steps', step, 'AverageTimeStep', elapsedTime/dble(step)
				write(*,*) 'computationTime', time2-time1
				write(*,*)
			end if
		end if

		!***********************************************************************
		!Final step: release all cascades into coarse mesh
		!(deallocate all extra memory here)
		!***********************************************************************
		!<delete defects in fine mesh
		do while(associated(ActiveCascades))
			CascadeCurrent=>ActiveCascades
			cascadeCell=CascadeCurrent%cellNumber
			call releaseFineMeshDefects(CascadeCurrent)
			call resetReactionListSingleCell(cascadeCell)
		end do

		call deallocateDefectList()
		call deallocateReactionList()
		call deallocateBoundarydefectList()
		deallocate(totalRateVol)
		if(allocated(listVI)) then
			deallocate(listVI)
		end if

		if(myProc%taskid==MASTER) then
			if(totdatToggle=='yes') close(TOTFILE)
			if(defectToggle=='yes') close(DEFFILE)
			if(stadatToggle=='yes') close(STAFILE)
			if(xyzdatToggle=='yes') close(XYZFILE)
		end if

	end do	!End of loop for multiple trials

	!***********************************************************************
	!Last step: deallocating memory
	!***********************************************************************
	deallocate(cascadeConnectivity)
	call deallocateCascadeList()
	call deallocateMaterialInput()

	call cpu_time(time2)
	runTime2=MPI_WTIME()
	if(myProc%taskid==MASTER) then
		!write(*,*) 'computation time', time2-time1
		write(*,*) 'run time', runTime2-runTime1
		write(*,*) 'communication time', commTimeSum
		write(*,*) 'Deallocating memory'
	end if

	call MPI_FINALIZE(ierr)		!must be performed at end of simulation to have no error

end program

max_stars_repo_name: MISA-SCD/misa-scd

max_issues_repo_name: MISA-SCD/misa-scd

max_forks_repo_name: MISA-SCD/misa-scd


text:
      subroutine opacit (modeop,waveop)
c******************************************************************************
c     This routine calculates the continuous opacity *kaplam* or
c     *kapref* and optical depth *taulam* as needed. the formulas are
c     ones lifted from the write-up of program *atlas*
c******************************************************************************

      implicit real*8 (a-h,o-z)
      include 'Atmos.com'
      include 'Kappa.com'
      include 'Dummy.com'
      real*8 taucheck(100), tauinteg(100)

c*****initialization
      freq = 2.997925d18/waveop
      freqlg = dlog(freq)
      do i = 1,ntau
         hkt = 6.6256d-27/(1.38065d-16*t(i))
         evhkt(i) = dexp(-freq*hkt)
      enddo
      do i=1,ntau
         aH1(i) =      1.d-99
         aHminus(i) =  1.d-99
         aHeminus(i) = 1.d-99
         aC1(i) =      1.d-99
         aMg1(i) =     1.d-99
         aMg2(i) =     1.d-99
         aAl1(i) =     1.d-99
         aSi1(i) =     1.d-99
         aSi2(i) =     1.d-99
         aFe1(i) =     1.d-99
      enddo

c*****compute the opacities
       call opacH1
       call opacHminus
       call opacHscat
       call opacH2scat
       call opacHeminus
       call opacHescat
       call opacescat
       call opacC1
       call opacMg1
       call opacMg2
       call opacAl1
       call opacSi1
       call opacSi2
       call opacFe1

c*****sum up all the opacities
      do i=1,ntau
         kaplamabs(i) = aH1(i) + aHminus(i) + aHeminus(i) + aC1(i) +
     .                  aMg1(i) + aMg2(i) + aAl1(i) + aSi1(i) +
     .                  aSi2(i) + aFe1(i)
         kaplamsca(i) = sigH(i) + sigH2(i) + sigHe(i) + sigel(i)
         kaplam(i)    = kaplamabs(i) + kaplamsca(i)
      enddo

c*****write out the opacities
      if (modprintopt .gt. 1) then
         write (nf1out,1001) waveop
         do i=1,ntau
            write (nf1out,1002) i, nint(t(i)), dlog10(kaplam(i)),
     .      dlog10(aH1(i)), dlog10(aHminus(i)), dlog10(sigH(i)),
     .      dlog10(aHeminus(i)), dlog10(sigHe(i)), dlog10(sigel(i)),
     .      dlog10(sigH2(i))
         enddo
         write (nf1out,1003) waveop
         do i=1,ntau
            write (nf1out,1002) i, nint(t(i)), dlog10(kaplam(i)),
     .      dlog10(aC1(i)),  dlog10(aMg1(i)), dlog10(aMg2(i)),
     .      dlog10(aAl1(i)), dlog10(aSi1(i)), dlog10(aSi2(i)),
     .      dlog10(aFe1(i))
         enddo
      endif

c*****add in an arbitrary amount of "extra" continuous opacity;
c     this is a pure fudge factor, so experiment with care
      if (fudge .gt. 0.0) then
         do i=1,ntau
            kaplam(i) = kaplam(i)*((fudge*10000)/t(i))
            if (i .eq. 1) then
               write(nf1out,1005)
               write(nf1out,1006) fudge
            endif
         enddo
      endif

c*****compute an optical depth array at this wavelength, and exit
      if (modeop .ne. 1) then
         do  i=1,ntau
            dummy1(i) = tauref(i)*kaplam(i)/(0.4343*kapref(i))
         enddo
         first = tauref(1)*kaplam(1)/kapref(1)
         dummy1(1) = rinteg(xref,dummy1,taulam,ntau,first)
         taulam(1) = first
         do i=2,ntau
            taulam(i) = taulam(i-1) + taulam(i)
         enddo
         if (modprintopt .gt. 1) write(nf1out,1004) (taulam(i),i=1,ntau)
         return
      endif

c*****here is the assignment of opacities kaplam to kapref; exit normally
      do i=1,ntau
         kapref(i) = kaplam(i)
      enddo
      if(modprintopt .lt. 2) return

c*****here is an internal tauref check on the externally read in tau scale;
c     then exit normally
      do i=1,ntau
         dummy1(i) = dlog10(tauref(i))
         tauinteg(i) = tauref(i)*kaplam(i)/(0.4343*kapref(i))
      enddo
      first = tauref(1)
      tauinteg(1) = rinteg(dummy1,tauinteg,taucheck,ntau,first)
      taucheck(1) = first
      do i=2,ntau
         taucheck(i) = taucheck(i-1) + taucheck(i)
      enddo
      write (nf1out,1004) (taucheck(i),i=1,ntau)
      return

c*****format statements
1001  format (/'log opacities due to H, He, e- at a wavelength of ',
     .       f10.2,' A'/' i  T(i)  kaplam     aH1     aH-    sigH',
     .       '    aHe-   sigHe   sigel   sigH2')
1002  format (i2,i6,10f8.2)
1003  format (/'log opacities due to "metal" edges at a wavelength',
     .       ' of ',f10.2,' A'/' i  T(i)  kaplam     aC1    aMg1',
     .       '    aMg2    aAl1    aSi1    aSi2    aFe1')
1004  format(/'COMPUTED TAULAM ARRAY'/(1p7e11.3))
1005  format ('Opacities artificially are now ',
     .        'altered by scaling law of: ')
1006  format ('      kaplam(i) = kaplam(i) * ',
     .        '((', f6.2, ' * 10000)/t(i))')

      end

ext: f

max_stars_repo_name: ipelupessy/FASMA-synthesis

max_issues_repo_name: ipelupessy/FASMA-synthesis

max_forks_repo_name: ipelupessy/FASMA-synthesis


text:
c take code from zchoose4 to make nice plots of postage stamps
c from spectra

      program plotpstamp

      parameter (maxcol=10000,maxrow=400)
      parameter (maxpstampx=1000, maxpstampy=500)
      parameter (maxaxis=5)
c finagle plot order (works only for columns varying faster than
c rows)
      parameter (IFCOLTRICK=0)
c skip plotting in bottom row?
      parameter (ISKIPBOTTOM=1)
c put filenames above imaging pstamp?
      parameter (IPRINTNAME=0)
c change stretch for each spectrum?
      parameter (ISPECSTRETCH=1)
c change stretch for each image?
      parameter (IIMSTRETCH=0)
c fatter, shorter spectra?
      parameter (IFATSPEC=1)

      real specdata(maxcol*maxrow)
      real spec2d(maxcol,maxrow)
      real data1d(maxcol)
      real dwave1(maxcol,maxrow)
      real wave1(maxcol)
      integer laxes1(maxaxis),laxes2(maxaxis)
      logical flag

      real pstamp(maxpstampx,maxpstampy)
c      real pstasmooth(maxpstampx,maxpstampy)
      real trans(6)

      character fname1*160,fname2*160,xlabel*30,toplabel*50

c      include 'pcredshift.h'

c change this if need to invert background/out of bounds color
      badflag = -10000.

c coord transf for image plotting
      trans(1) = 0.
      trans(2) = 1.
      trans(3) = 0.
      trans(4) = 0.
      trans(5) = 0.
      trans(6) = 1.

c image scales, arcsec/pixel
      deimospix = 0.119
      cfhtpix = 0.206

      nplot=0
      ierr=0

      idwin=pgopen('?')
      write(*,'("pgplot nx by ny: ",$)')
      read(*,*) nxplot,nyplot
      call pgsubp(nxplot,nyplot)
      nxplotold = nxplot
      nxplot = abs(nxplot)
      call pgpage()
      call pgscf(2)
      call pgsch(2.7)
      call pgscir(16,254)

      ixpl = 1
      iypl = 1
      call pgpanl(ixpl,iypl)

      write(*,'("y radius (Deimos pixels): ",$)')
      read (*,*) yrad
      write(*,'("km/sec radius: " ,$)')
      read(*,*) velrad
      write(*,'("Deimos flux min, max: ",$)')
      read(*,*) fmin,fmax
      write(*,'("CFHT flux min, max: ",$)')
      read(*,*) fmin_im,fmax_im

 100  continue
      write(*,'("2d spectrum file [quit]: ",$)')
      read(*,'(a)') fname1
      if (fname1(1:3) .eq. '   ') go to 666

      write(*,'("ypos (pixels): ",$)')
      read(*,*) ypos
      iypos = int(ypos)
      iyrad = int(yrad)
c      write(*,'("wavelength radius (A): ",$)')
c      read(*,*) wrad
      write(*,'("rest wavelength and z: ",$)')
      read(*,*) wrest,z
      wcen = wrest*(1.+z)

      if (ISPECSTRETCH .eq. 1) then
         write(*,'("spectrum fmin, fmax: ",$)')
         read(*,*) fmin,fmax
      end if

      call ftgiou(im1,ierr)
      call doerr(ierr)
      call ftopen(im1,fname1,0,block,ierr)
      if (ierr .ne. 0) then
         write(*,'(a,a)') "Couldnt open image ",fname1
         go to 100
      end if
      call doerr(ierr)

c retrieve 2-d data from hdu 2
      ihdu = 2

      call readfield(im1,ihdu,'FLUX',itype1,isize1,
     $     naxis1,laxes1,specdata)
      call reformat(isize1,specdata,naxis1,laxes1,maxcol,maxrow,spec2d)
      nx1=laxes1(1)
      ny1=laxes1(2)
      call readfield(im1,ihdu,'DLAMBDA',itype1,isize1,
     $     naxis1,laxes1,specdata)
      call reformat(isize1,specdata,naxis1,laxes1,maxcol,maxrow,dwave1)
      call readfield(im1,ihdu,'LAMBDA0',itype1,iw0size1,
     $     nw0axis1,lwaxis1,wave1)

      call ftclos(im1,ierr)
      call ftfiou(im1,ierr)
      call doerr(ierr)
      do i=1,nx1
         wave1(i) = wave1(i) + dwave1(i,iypos)
      end do
      dw1 = (wave1(nx1) - wave1(1)) / (nx1-1.)
      w01 = wave1(1) - dw1
c      write(*,'(a80,2x,f7.2,2x,f5.3)') fname1,w01,dw1

      call findwave(nx1,wave1,wcen,ixcen1)
c      xrad = wrad /dw1
      xrad = velrad / 3.0e5 * wcen / dw1
      ixrad = int(xrad)
      write(*,*) ixcen1,ixrad,iypos,iyrad

      ioffset = ixcen1-ixrad-1
      joffset = iypos-iyrad-1
      do j=iypos-iyrad,iypos+iyrad
         jpsta = j-joffset
         do i=ixcen1-ixrad,ixcen1+ixrad
            ipsta = i-ioffset
            if (i .ge. 1 .and. i .le. nx1 .and.
     $           j .ge. 1 .and. j .le. ny1) then
               pstamp(ipsta,jpsta) = spec2d(i,j)
            else
               pstamp(ipsta,jpsta) = badflag
            end if
         end do
      end do

      nplotx = 2*ixrad+1
      nploty = 2*iyrad+1

c      call pgqvp(0,x1,x2,y1,y2)
c      write(*,*) x1,x2,y1,y2
c      call pgsvp(0.1,0.5,0.18,0.82)
      if (IFATSPEC .eq. 0) then
         call pgsvp(0.1,0.62,0.3,0.7)
      else
         call pgsvp(0.1,0.53,0.2,0.8)
      end if
      nplot=nplot+1

c      call pgenv(0.5,nplotx+0.5,0.5,nploty+0.5,1,0)
c square pixels
c      call pgenv(0.5,nplotx+0.5,0.5,nploty+0.5,1,-2)
c non-square pixels
c      call pgenv(0.5,nplotx+0.5,0.5,nploty+0.5,0,-2)
c      call pgenv(0.5,nplotx+0.5,0.5,nploty+0.5,0,-2)
      call pgswin(0.5,nplotx+0.5,0.5,nploty+0.5)
c arcsec per y pixel
      arcppix = deimospix
      yarcsecrad = yrad*arcppix
c      trans(1) = wcen-wrad
c      trans(2) = dw1
c      trans(3) = 0.0
c      trans(4) = -yarcsecrad
c      trans(5) = 0.0
c      trans(6) = arcppix
c      call pgenv(wcen-wrad,wcen+wrad,-yarcsecrad,yarcsecrad,0,0)
c 0 for linear
      call pgsitf(2)
      call pggray(pstamp,maxpstampx,maxpstampy,
c     $     ixpmin,ixpmax,iypmin,iypmax,
     $     1,nplotx,1,nploty,
     $     fmin,fmax,trans)
c      call pgswin(wcen-wrad,wcen+wrad,-yarcsecrad,yarcsecrad)
      call pgswin(-velrad,velrad,-yarcsecrad,yarcsecrad)
      write(xlabel,'("km/sec")')
      write(toplabel,'(i4," at z = ",f6.4)') nint(wrest),z
c      call pglabel(xlabel, 'arcsec', toplabel)
c      call pglabel(' ','arcsec',' ')
c if plot is on left hand column, write y label
      if (ixpl .eq. 1)
     $     call pgmtxt('L',2.0,0.5,0.5,"arcsec")
c if plot is on bottom, write x label
      if ((ISKIPBOTTOM .eq. 0 .and. iypl .eq. nyplot) .or.
     $     (ISKIPBOTTOM .eq. 1 .and. iypl .eq. nyplot-1))
     $     call pgmtxt('B',2.6,0.5,0.5,'km/sec')
      call pgmtxt('T',1.1,0.0,0.0,toplabel)

c ticks out
c      call pgbox('BCINST',0.0,0,'BCINST',0.0,0)
      call pgbox('BCNST',0.0,0,'BCNST',0.0,0)

c temp
c      go to 100

c now read in an image
 200  continue
      write(*,'("image fits file: ",$)')
      read(*,'(a)') fname1

      call ftgiou(im1,ierr)
      call doerr(ierr)
      call ftopen(im1,fname1,0,block,ierr)
      if (ierr .ne. 0) then
         write(*,'(a,a)') "Couldnt open image ",fname1
         go to 200
      end if
      call doerr(ierr)

      call ftgknj(im1,'NAXIS',1,7,laxes1,idim,ierr)
      call doerr(ierr)
      nx = laxes1(1)
      ny = laxes1(2)
c      write(*,'("xmin, xmax, ymin, ymax: ",$)')
c      read (*,*) xmin,xmax,ymin,ymax

c or compute it automatically, centering the image and making
c its size match the deimos slit
      yrad2 = yrad * deimospix/cfhtpix
      yarcsecrad2 = yrad2 * cfhtpix
      xrad2 = yrad2
      xmin = (nx+1)/2. -xrad2
      xmax = (nx+1)/2. +xrad2
      ymin = (ny+1)/2. -yrad2
      ymax = (ny+1)/2. +yrad2

      ixmin = int(xmin)
      ixmax = int(xmax)+1
      iymin = int(ymin)
      iymax = int(ymax)+1
      nplotx = ixmax-ixmin+1
      nploty = iymax-iymin+1

      if (IIMSTRETCH .eq. 1) then
         write(*,'("image fmin, fmax: ",$)')
         read(*,*) fmin_im,fmax_im
      end if

c      write(*,*) nx,ny,nplotx,nploty,yrad2,yarcsecrad2,
c     $     ixmin,ixmax,iymin,iymax

      do j=iymin,iymax
         jpsta = j-iymin+1
         if (j .lt. 1 .or. j. gt. ny) then
            do i=ixmin,ixmax
               ipsta = i-ixmin+1
               pstamp(ipsta,jpsta) = badflag
            end do
         else
            ipixel = (j-1)*nx + 1
            call ftgpve(im1,0,ipixel,nx,0.0,data1d,flag,ierr)
            call doerr(ierr)
            do i=ixmin,ixmax
               ipsta = i-ixmin+1
               if (i .lt. 1 .or. i .gt. nx) then
                  pstamp(ipsta,jpsta) = badflag
               else
                  pstamp(ipsta,jpsta) = data1d(i)
               end if
            end do
         end if
      end do

c      call pgsvp(0.6,0.9,0.18,0.82)
      if (IFATSPEC .eq. 0) then
         call pgsvp(0.68,0.92,0.3,0.7)
      else
         call pgsvp(0.58,0.92,0.2,0.8)
      end if
c      call pgswin(0.5,nplotx+0.5,0.5,nploty+0.5)
c want square pixels, but don't want it to skip to next panel ...
c      call pgenv(0.5,nplotx+0.5,0.5,nploty+0.5,1,-2)
      call pgwnad(0.5,nplotx+0.5,0.5,nploty+0.5)
c      call pgswin(0.5,nplotx+0.5,0.5,nploty+0.5)

c pgsitf(0) = linear, 1= log, 2= square root
      call pgsitf(2)
      call pggray(pstamp,maxpstampx,maxpstampy,
c     $     ixpmin,ixpmax,iypmin,iypmax,
     $     1,nplotx,1,nploty,
     $     fmin_im,fmax_im,trans)
      call pgswin(-yarcsecrad2,+yarcsecrad2,-yarcsecrad2,+yarcsecrad2)
c ticks out
c      call pgbox('BCINST',0.0,0,'BCIMST',0.0,0)
      call pgbox('BCNST',0.0,0,'BCMST',0.0,0)
c      call pglabel('arcsec' , '', '')

c      irow = mod(nplot,nyplot)
c      if (irow .eq. 0) irow = nyplot
c      icol = int(nplot / nyplot) +1
c      write(*,*) "icol, irow = ",icol,irow
      write(*,*) "icol, irow = ",ixpl,iypl

c if plot is on bottom, write x label
      if ((ISKIPBOTTOM .eq. 0 .and. iypl .eq. nyplot) .or.
     $     (ISKIPBOTTOM .eq. 1 .and. iypl .eq. nyplot-1))
     $     call pgmtxt('B',2.6,0.5,0.5,'arcsec')
c temporary for id'ing object
      if (IPRINTNAME .eq. 1) call pgmtxt('T',1.2,0.0,0.0,fname1)

      if (IFCOLTRICK .eq. 0) then
c don't finagle plot order
         call pgpage()
         if (nxplotold .ge. 0) then
c fill rows first
            if (ixpl .lt. nxplot) then
               ixpl=ixpl+1
            else
               ixpl=1
               if(iypl .lt. nyplot) then
                  iypl = iypl+1
               else
                  iypl=1
               end if
            end if
         else
c fill columns first
            if((ISKIPBOTTOM .eq. 0 .and. iypl .lt. nyplot) .or.
     $           (ISKIPBOTTOM .eq. 1 .and. iypl .lt. nyplot-1)) then
               iypl = iypl+1
            else
               iypl=1
               if(ixpl .lt. nxplot) then
                  ixpl=ixpl+1
               else
                  ixpl=1
               end if
            end if
         end if
      else
c finagling, column order
c if we are at the next-to-bottom row
         if (ISKIPBOTTOM .eq. 1 .and. iypl .eq. nyplot-1) then
            iypl=1
            if (ixpl .lt. nxplot) then
               ixpl=ixpl+1
            else
c advance to next page by going to last panel and calling pgpage
               call pgpanl(nxplot,nyplot)
               call pgpage()
               ixpl=1
            end if
            call pgpanl(ixpl,iypl)
         else
            call pgpage()
            if(iypl .lt. nyplot) then
               iypl = iypl+1
            else
               iypl = 1
               if (ixpl .lt. nxplot) then
                  ixpl = ixpl+1
               else
                  ixpl = 1
               end if
            end if
         end if
      end if

      go to 100

 666  continue

      call pgclos

      end

cccccccccccccccccccccccccccccc
c  find the index closest to a given wave in the wavelength array,
c  assumed increasing
      subroutine findwave(nw,warr,w,iw)

      real warr(nw)

c  define a buffer in Angstroms, so we can plot even if slightly off end
c  of spectrum
      wbuff = 5.0
      if (warr(1)-wbuff .gt. w) then
c iw=-1 means w not in array
         iw=-1
         return
      end if
      j=2
 100  continue
      if(warr(j) .lt. w) then
         j=j+1
c test for reaching the end
         if (j .gt. nw) then
            if (warr(nw)+wbuff .gt. w) then
c just barely off
               iw=nw
            else
               iw=-1
            end if
            return
         end if
c if not at end, try again
         go to 100
      else
c warr(j) > w, we've gone past w.  Pick the closest one
         if (abs(warr(j)-w) .lt. abs(warr(j-1)-w)) then
            iw = j
         else
            iw=j-1
         end if
         return
      end if

      end

cccccccccccccccccccccccccccccc
c  reformat a 1-d data array into proper 2-d
      subroutine reformat(n1,data1d,naxis,laxes,nxphys,nyphys,data2d)

      real data1d(n1),data2d(nxphys,nyphys)
      integer laxes(naxis)

      if (naxis .le. 0) then
         write(*,*) "Error in reformat, naxis <= 0"
         return
      else if (naxis .eq. 1) then
         nx = laxes(1)
         ny = 1
      else
         nx = laxes(1)
         ny = laxes(2)
      end if
c bounds checking
      if (nx .gt. nxphys) then
         write(*,'("Warning! 2d image x is ",i6," max ",i6)') nx,nxphys
         nx=nxphys
      end if
      if (ny .gt. nyphys) then
         write(*,'("Warning! 2d image y is ",i6," max ",i6)') ny,nyphys
         ny=nyphys
      end if

      do j=1,ny
         ioffset = (j-1)*nx
         do i=1,nx
            data2d(i,j) = data1d(ioffset+i)
         end do
      end do
      return
      end


! Copyright 2005-2016 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
!
! In applying this licence, ECMWF does not waive the privileges and immunities granted to it by
! virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
!
!
!
!  Description: how to copy a message in memory
!
!
!
program copy
  use grib_api
  implicit none
  integer                            :: err, centre
  integer(kind=kindOfSize)           :: byte_size
  integer                            :: infile,outfile
  integer                            :: igrib_in
  integer                            :: igrib_out
  character(len=1), dimension(:), allocatable :: message

  call grib_open_file(infile,'../../data/constant_field.grib1','r')
  call grib_open_file(outfile,'out.copy.grib1','w')

  !     a new grib message is loaded from file
  !     igrib is the grib id to be used in subsequent calls
  call grib_new_from_file(infile,igrib_in)

  call grib_get_message_size(igrib_in, byte_size)
  allocate(message(byte_size), stat=err)

  call grib_copy_message(igrib_in,message)

  call grib_new_from_message(igrib_out, message)

  centre=80
  call grib_set(igrib_out,"centre",centre)

!  write messages to a file
  call grib_write(igrib_out,outfile)

  call grib_release(igrib_out)

  call grib_release(igrib_in)

  call grib_close_file(infile)
  call grib_close_file(outfile)
  deallocate(message)

end program copy


MODULE cbmz_bb_Jacobian
  USE cbmz_bb_Parameters
  USE cbmz_bb_JacobianSP
  IMPLICIT NONE
CONTAINS
SUBROUTINE Jac_SP_Vec ( JVS, UV, JUV )
  REAL(kind=dp) :: JVS(LU_NONZERO)
  REAL(kind=dp) :: UV(NVAR)
  REAL(kind=dp) :: JUV(NVAR)
  JUV(1) = JVS(1)*UV(1)+JVS(2)*UV(7)+JVS(3)*UV(50)
  JUV(2) = JVS(4)*UV(2)
  JUV(3) = JVS(5)*UV(3)
  JUV(4) = JVS(6)*UV(4)+JVS(7)*UV(23)+JVS(8)*UV(34)+JVS(9)*UV(51)
  JUV(5) = JVS(10)*UV(5)+JVS(11)*UV(34)+JVS(12)*UV(35)+JVS(13)*UV(39)+JVS(14)*UV(45)+JVS(15)*UV(51)+JVS(16)*UV(53)&
             &+JVS(17)*UV(55)
  JUV(6) = JVS(18)*UV(6)+JVS(19)*UV(51)
  JUV(7) = JVS(20)*UV(7)+JVS(21)*UV(50)
  JUV(8) = JVS(22)*UV(8)+JVS(23)*UV(50)
  JUV(9) = JVS(24)*UV(9)+JVS(25)*UV(50)+JVS(26)*UV(55)
  JUV(10) = JVS(27)*UV(10)+JVS(28)*UV(52)+JVS(29)*UV(53)
  JUV(11) = JVS(30)*UV(11)+JVS(31)*UV(50)
  JUV(12) = JVS(32)*UV(12)+JVS(33)*UV(52)+JVS(34)*UV(54)
  JUV(13) = JVS(35)*UV(13)+JVS(36)*UV(50)
  JUV(14) = JVS(37)*UV(14)+JVS(38)*UV(34)+JVS(39)*UV(39)+JVS(40)*UV(50)+JVS(41)*UV(51)
  JUV(15) = JVS(42)*UV(15)+JVS(43)*UV(26)+JVS(44)*UV(50)+JVS(45)*UV(52)+JVS(46)*UV(54)
  JUV(16) = JVS(47)*UV(16)+JVS(48)*UV(50)+JVS(49)*UV(52)+JVS(50)*UV(55)
  JUV(17) = JVS(51)*UV(17)+JVS(52)*UV(34)+JVS(53)*UV(39)+JVS(54)*UV(40)+JVS(55)*UV(50)+JVS(56)*UV(51)
  JUV(18) = JVS(57)*UV(11)+JVS(58)*UV(13)+JVS(59)*UV(18)+JVS(60)*UV(50)+JVS(61)*UV(56)
  JUV(19) = JVS(62)*UV(19)+JVS(63)*UV(29)+JVS(64)*UV(34)+JVS(65)*UV(39)+JVS(66)*UV(41)+JVS(67)*UV(47)+JVS(68)*UV(48)&
              &+JVS(69)*UV(49)+JVS(70)*UV(50)+JVS(71)*UV(51)+JVS(72)*UV(54)+JVS(73)*UV(56)
  JUV(20) = JVS(74)*UV(20)+JVS(75)*UV(44)+JVS(76)*UV(50)+JVS(77)*UV(55)
  JUV(21) = JVS(78)*UV(21)+JVS(79)*UV(40)+JVS(80)*UV(50)+JVS(81)*UV(55)
  JUV(22) = JVS(82)*UV(22)+JVS(83)*UV(50)+JVS(84)*UV(52)+JVS(85)*UV(55)+JVS(86)*UV(56)
  JUV(23) = JVS(87)*UV(23)+JVS(88)*UV(50)+JVS(89)*UV(51)
  JUV(24) = JVS(90)*UV(24)+JVS(91)*UV(34)+JVS(92)*UV(39)+JVS(93)*UV(44)+JVS(94)*UV(50)+JVS(95)*UV(51)
  JUV(25) = JVS(96)*UV(6)+JVS(97)*UV(25)+JVS(98)*UV(51)+JVS(99)*UV(52)+JVS(100)*UV(54)+JVS(101)*UV(56)
  JUV(26) = JVS(102)*UV(11)+JVS(103)*UV(13)+JVS(104)*UV(26)+JVS(105)*UV(50)+JVS(106)*UV(54)
  JUV(27) = JVS(107)*UV(12)+JVS(108)*UV(26)+JVS(109)*UV(27)+JVS(110)*UV(36)+JVS(111)*UV(42)+JVS(112)*UV(43)+JVS(113)&
              &*UV(45)+JVS(114)*UV(50)+JVS(115)*UV(52)+JVS(116)*UV(54)+JVS(117)*UV(55)
  JUV(28) = JVS(118)*UV(23)+JVS(119)*UV(28)+JVS(120)*UV(30)+JVS(121)*UV(33)+JVS(122)*UV(34)+JVS(123)*UV(35)+JVS(124)&
              &*UV(36)+JVS(125)*UV(39)+JVS(126)*UV(42)+JVS(127)*UV(43)+JVS(128)*UV(45)+JVS(129)*UV(50)+JVS(130)*UV(51)&
              &+JVS(131)*UV(54)+JVS(132)*UV(56)
  JUV(29) = JVS(133)*UV(13)+JVS(134)*UV(19)+JVS(135)*UV(29)+JVS(136)*UV(31)+JVS(137)*UV(32)+JVS(141)*UV(45)+JVS(145)&
              &*UV(50)+JVS(147)*UV(54)+JVS(148)*UV(55)+JVS(149)*UV(56)
  JUV(30) = JVS(150)*UV(18)+JVS(151)*UV(26)+JVS(152)*UV(30)+JVS(153)*UV(50)+JVS(154)*UV(51)+JVS(156)*UV(56)
  JUV(31) = JVS(157)*UV(31)+JVS(158)*UV(35)+JVS(159)*UV(54)+JVS(160)*UV(55)+JVS(161)*UV(56)
  JUV(32) = JVS(162)*UV(32)+JVS(163)*UV(35)+JVS(164)*UV(50)+JVS(165)*UV(55)+JVS(166)*UV(56)
  JUV(33) = JVS(167)*UV(33)+JVS(168)*UV(45)+JVS(169)*UV(50)+JVS(170)*UV(55)+JVS(171)*UV(56)
  JUV(34) = JVS(172)*UV(34)+JVS(173)*UV(50)+JVS(174)*UV(51)+JVS(175)*UV(54)
  JUV(35) = JVS(176)*UV(35)+JVS(177)*UV(50)+JVS(178)*UV(51)+JVS(179)*UV(54)
  JUV(36) = JVS(180)*UV(20)+JVS(181)*UV(23)+JVS(182)*UV(24)+JVS(183)*UV(30)+JVS(184)*UV(32)+JVS(185)*UV(33)+JVS(186)&
              &*UV(34)+JVS(187)*UV(35)+JVS(188)*UV(36)+JVS(190)*UV(41)+JVS(191)*UV(44)+JVS(192)*UV(45)+JVS(193)*UV(46)&
              &+JVS(194)*UV(50)+JVS(195)*UV(51)+JVS(196)*UV(54)+JVS(198)*UV(56)
  JUV(37) = JVS(199)*UV(11)+JVS(200)*UV(13)+JVS(201)*UV(23)+JVS(202)*UV(26)+JVS(203)*UV(30)+JVS(204)*UV(34)+JVS(205)&
              &*UV(35)+JVS(206)*UV(37)+JVS(207)*UV(39)+JVS(208)*UV(42)+JVS(209)*UV(45)+JVS(210)*UV(47)+JVS(211)*UV(48)&
              &+JVS(212)*UV(49)+JVS(213)*UV(50)+JVS(214)*UV(51)+JVS(215)*UV(54)+JVS(216)*UV(55)+JVS(217)*UV(56)
  JUV(38) = JVS(218)*UV(33)+JVS(219)*UV(38)+JVS(220)*UV(39)+JVS(221)*UV(45)+JVS(222)*UV(46)+JVS(223)*UV(47)+JVS(224)&
              &*UV(48)+JVS(225)*UV(49)+JVS(226)*UV(50)+JVS(227)*UV(51)+JVS(228)*UV(54)+JVS(230)*UV(56)
  JUV(39) = JVS(231)*UV(39)+JVS(232)*UV(50)+JVS(233)*UV(51)+JVS(234)*UV(54)
  JUV(40) = JVS(235)*UV(17)+JVS(236)*UV(21)+JVS(237)*UV(34)+JVS(238)*UV(39)+JVS(239)*UV(40)+JVS(240)*UV(47)+JVS(241)&
              &*UV(48)+JVS(242)*UV(49)+JVS(243)*UV(50)+JVS(244)*UV(51)+JVS(245)*UV(54)+JVS(246)*UV(55)+JVS(247)*UV(56)
  JUV(41) = JVS(248)*UV(34)+JVS(249)*UV(39)+JVS(250)*UV(41)+JVS(251)*UV(49)+JVS(252)*UV(50)+JVS(254)*UV(54)+JVS(255)&
              &*UV(55)+JVS(256)*UV(56)
  JUV(42) = JVS(257)*UV(13)+JVS(258)*UV(30)+JVS(259)*UV(33)+JVS(260)*UV(34)+JVS(261)*UV(39)+JVS(262)*UV(42)+JVS(263)&
              &*UV(45)+JVS(264)*UV(46)+JVS(265)*UV(47)+JVS(266)*UV(50)+JVS(267)*UV(51)+JVS(268)*UV(54)+JVS(270)*UV(56)
  JUV(43) = JVS(271)*UV(8)+JVS(272)*UV(21)+JVS(273)*UV(23)+JVS(274)*UV(30)+JVS(275)*UV(31)+JVS(276)*UV(33)+JVS(277)&
              &*UV(34)+JVS(278)*UV(35)+JVS(279)*UV(39)+JVS(280)*UV(40)+JVS(281)*UV(41)+JVS(282)*UV(43)+JVS(283)*UV(45)&
              &+JVS(284)*UV(47)+JVS(285)*UV(48)+JVS(286)*UV(49)+JVS(287)*UV(50)+JVS(288)*UV(51)+JVS(289)*UV(54)+JVS(291)&
              &*UV(56)
  JUV(44) = JVS(292)*UV(14)+JVS(293)*UV(20)+JVS(294)*UV(34)+JVS(295)*UV(38)+JVS(296)*UV(39)+JVS(297)*UV(43)+JVS(298)&
              &*UV(44)+JVS(299)*UV(45)+JVS(304)*UV(50)+JVS(305)*UV(51)+JVS(306)*UV(53)+JVS(307)*UV(54)+JVS(308)*UV(55)&
              &+JVS(309)*UV(56)
  JUV(45) = JVS(310)*UV(31)+JVS(311)*UV(32)+JVS(312)*UV(35)+JVS(313)*UV(45)+JVS(314)*UV(50)+JVS(315)*UV(51)+JVS(316)&
              &*UV(54)+JVS(318)*UV(56)
  JUV(46) = JVS(319)*UV(38)+JVS(320)*UV(39)+JVS(322)*UV(46)+JVS(326)*UV(50)+JVS(327)*UV(51)+JVS(328)*UV(54)+JVS(329)&
              &*UV(55)+JVS(330)*UV(56)
  JUV(47) = JVS(331)*UV(32)+JVS(332)*UV(33)+JVS(335)*UV(46)+JVS(336)*UV(47)+JVS(337)*UV(48)+JVS(339)*UV(50)+JVS(342)&
              &*UV(55)
  JUV(48) = JVS(344)*UV(29)+JVS(347)*UV(34)+JVS(349)*UV(39)+JVS(352)*UV(47)+JVS(353)*UV(48)+JVS(355)*UV(50)+JVS(356)&
              &*UV(51)+JVS(357)*UV(54)+JVS(358)*UV(55)+JVS(359)*UV(56)
  JUV(49) = JVS(360)*UV(15)+JVS(361)*UV(18)+JVS(363)*UV(31)+JVS(364)*UV(32)+JVS(366)*UV(41)+JVS(367)*UV(45)+JVS(368)&
              &*UV(48)+JVS(369)*UV(49)+JVS(370)*UV(50)+JVS(372)*UV(52)+JVS(373)*UV(54)+JVS(374)*UV(55)+JVS(375)*UV(56)
  JUV(50) = JVS(376)*UV(6)+JVS(377)*UV(7)+JVS(378)*UV(8)+JVS(379)*UV(9)+JVS(380)*UV(11)+JVS(381)*UV(13)+JVS(382)*UV(14)&
              &+JVS(383)*UV(16)+JVS(384)*UV(17)+JVS(385)*UV(20)+JVS(386)*UV(21)+JVS(387)*UV(22)+JVS(388)*UV(23)+JVS(389)&
              &*UV(24)+JVS(390)*UV(26)+JVS(391)*UV(27)+JVS(392)*UV(28)+JVS(393)*UV(29)+JVS(394)*UV(30)+JVS(398)*UV(34)&
              &+JVS(399)*UV(35)+JVS(400)*UV(36)+JVS(401)*UV(38)+JVS(402)*UV(39)+JVS(405)*UV(42)+JVS(406)*UV(43)+JVS(408)&
              &*UV(45)+JVS(410)*UV(47)+JVS(412)*UV(49)+JVS(413)*UV(50)+JVS(414)*UV(51)+JVS(415)*UV(52)+JVS(417)*UV(54)&
              &+JVS(418)*UV(55)+JVS(419)*UV(56)
  JUV(51) = JVS(420)*UV(23)+JVS(421)*UV(25)+JVS(422)*UV(30)+JVS(423)*UV(34)+JVS(424)*UV(35)+JVS(425)*UV(39)+JVS(426)&
              &*UV(45)+JVS(427)*UV(50)+JVS(428)*UV(51)+JVS(429)*UV(52)+JVS(430)*UV(53)+JVS(432)*UV(55)+JVS(433)*UV(56)
  JUV(52) = JVS(434)*UV(10)+JVS(435)*UV(12)+JVS(436)*UV(15)+JVS(437)*UV(16)+JVS(438)*UV(18)+JVS(439)*UV(22)+JVS(440)&
              &*UV(25)+JVS(442)*UV(27)+JVS(443)*UV(31)+JVS(444)*UV(32)+JVS(445)*UV(33)+JVS(448)*UV(37)+JVS(450)*UV(40)&
              &+JVS(451)*UV(41)+JVS(454)*UV(44)+JVS(456)*UV(46)+JVS(458)*UV(48)+JVS(459)*UV(49)+JVS(460)*UV(50)+JVS(461)&
              &*UV(51)+JVS(462)*UV(52)+JVS(463)*UV(53)+JVS(464)*UV(54)+JVS(465)*UV(55)+JVS(466)*UV(56)
  JUV(53) = JVS(467)*UV(10)+JVS(468)*UV(30)+JVS(469)*UV(34)+JVS(470)*UV(35)+JVS(471)*UV(38)+JVS(472)*UV(39)+JVS(473)&
              &*UV(42)+JVS(474)*UV(43)+JVS(475)*UV(45)+JVS(476)*UV(46)+JVS(480)*UV(50)+JVS(481)*UV(51)+JVS(482)*UV(52)&
              &+JVS(483)*UV(53)+JVS(484)*UV(54)+JVS(485)*UV(55)+JVS(486)*UV(56)
  JUV(54) = JVS(487)*UV(12)+JVS(488)*UV(25)+JVS(489)*UV(26)+JVS(490)*UV(27)+JVS(491)*UV(34)+JVS(492)*UV(35)+JVS(493)&
              &*UV(36)+JVS(494)*UV(37)+JVS(495)*UV(39)+JVS(496)*UV(40)+JVS(497)*UV(41)+JVS(498)*UV(42)+JVS(499)*UV(43)&
              &+JVS(500)*UV(44)+JVS(501)*UV(45)+JVS(502)*UV(46)+JVS(504)*UV(48)+JVS(506)*UV(50)+JVS(507)*UV(51)+JVS(508)&
              &*UV(52)+JVS(509)*UV(53)+JVS(510)*UV(54)+JVS(511)*UV(55)+JVS(512)*UV(56)
  JUV(55) = JVS(513)*UV(7)+JVS(514)*UV(8)+JVS(515)*UV(9)+JVS(516)*UV(11)+JVS(517)*UV(13)+JVS(518)*UV(16)+JVS(519)*UV(18)&
              &+JVS(520)*UV(20)+JVS(521)*UV(21)+JVS(522)*UV(23)+JVS(523)*UV(24)+JVS(524)*UV(26)+JVS(525)*UV(28)+JVS(526)&
              &*UV(30)+JVS(527)*UV(31)+JVS(528)*UV(32)+JVS(529)*UV(33)+JVS(530)*UV(34)+JVS(531)*UV(35)+JVS(532)*UV(36)&
              &+JVS(533)*UV(37)+JVS(534)*UV(39)+JVS(535)*UV(40)+JVS(536)*UV(41)+JVS(537)*UV(42)+JVS(538)*UV(43)+JVS(539)&
              &*UV(44)+JVS(540)*UV(45)+JVS(541)*UV(46)+JVS(542)*UV(47)+JVS(543)*UV(48)+JVS(544)*UV(49)+JVS(545)*UV(50)&
              &+JVS(546)*UV(51)+JVS(547)*UV(52)+JVS(548)*UV(53)+JVS(549)*UV(54)+JVS(550)*UV(55)+JVS(551)*UV(56)
  JUV(56) = JVS(552)*UV(18)+JVS(553)*UV(22)+JVS(554)*UV(25)+JVS(555)*UV(31)+JVS(556)*UV(32)+JVS(557)*UV(33)+JVS(559)&
              &*UV(37)+JVS(561)*UV(40)+JVS(562)*UV(41)+JVS(564)*UV(44)+JVS(566)*UV(46)+JVS(568)*UV(48)+JVS(570)*UV(50)&
              &+JVS(571)*UV(51)+JVS(572)*UV(52)+JVS(573)*UV(53)+JVS(574)*UV(54)+JVS(575)*UV(55)+JVS(576)*UV(56)
END SUBROUTINE Jac_SP_Vec
SUBROUTINE JacTR_SP_Vec ( JVS, UV, JTUV )
  REAL(kind=dp) :: JVS(LU_NONZERO)
  REAL(kind=dp) :: UV(NVAR)
  REAL(kind=dp) :: JTUV(NVAR)
  JTUV(1) = JVS(1)*UV(1)
  JTUV(2) = JVS(4)*UV(2)
  JTUV(3) = JVS(5)*UV(3)
  JTUV(4) = JVS(6)*UV(4)
  JTUV(5) = JVS(10)*UV(5)
  JTUV(6) = JVS(18)*UV(6)+JVS(96)*UV(25)+JVS(376)*UV(50)
  JTUV(7) = JVS(2)*UV(1)+JVS(20)*UV(7)+JVS(377)*UV(50)+JVS(513)*UV(55)
  JTUV(8) = JVS(22)*UV(8)+JVS(271)*UV(43)+JVS(378)*UV(50)+JVS(514)*UV(55)
  JTUV(9) = JVS(24)*UV(9)+JVS(379)*UV(50)+JVS(515)*UV(55)
  JTUV(10) = JVS(27)*UV(10)+JVS(434)*UV(52)+JVS(467)*UV(53)
  JTUV(11) = JVS(30)*UV(11)+JVS(57)*UV(18)+JVS(102)*UV(26)+JVS(199)*UV(37)+JVS(380)*UV(50)+JVS(516)*UV(55)
  JTUV(12) = JVS(32)*UV(12)+JVS(107)*UV(27)+JVS(435)*UV(52)+JVS(487)*UV(54)
  JTUV(13) = JVS(35)*UV(13)+JVS(58)*UV(18)+JVS(103)*UV(26)+JVS(133)*UV(29)+JVS(200)*UV(37)+JVS(257)*UV(42)+JVS(381)&
               &*UV(50)+JVS(517)*UV(55)
  JTUV(14) = JVS(37)*UV(14)+JVS(292)*UV(44)+JVS(382)*UV(50)
  JTUV(15) = JVS(42)*UV(15)+JVS(360)*UV(49)+JVS(436)*UV(52)
  JTUV(16) = JVS(47)*UV(16)+JVS(383)*UV(50)+JVS(437)*UV(52)+JVS(518)*UV(55)
  JTUV(17) = JVS(51)*UV(17)+JVS(235)*UV(40)+JVS(384)*UV(50)
  JTUV(18) = JVS(59)*UV(18)+JVS(150)*UV(30)+JVS(361)*UV(49)+JVS(438)*UV(52)+JVS(519)*UV(55)+JVS(552)*UV(56)
  JTUV(19) = JVS(62)*UV(19)+JVS(134)*UV(29)
  JTUV(20) = JVS(74)*UV(20)+JVS(180)*UV(36)+JVS(293)*UV(44)+JVS(385)*UV(50)+JVS(520)*UV(55)
  JTUV(21) = JVS(78)*UV(21)+JVS(236)*UV(40)+JVS(272)*UV(43)+JVS(386)*UV(50)+JVS(521)*UV(55)
  JTUV(22) = JVS(82)*UV(22)+JVS(387)*UV(50)+JVS(439)*UV(52)+JVS(553)*UV(56)
  JTUV(23) = JVS(7)*UV(4)+JVS(87)*UV(23)+JVS(118)*UV(28)+JVS(181)*UV(36)+JVS(201)*UV(37)+JVS(273)*UV(43)+JVS(388)*UV(50)&
               &+JVS(420)*UV(51)+JVS(522)*UV(55)
  JTUV(24) = JVS(90)*UV(24)+JVS(182)*UV(36)+JVS(389)*UV(50)+JVS(523)*UV(55)
  JTUV(25) = JVS(97)*UV(25)+JVS(421)*UV(51)+JVS(440)*UV(52)+JVS(488)*UV(54)+JVS(554)*UV(56)
  JTUV(26) = JVS(43)*UV(15)+JVS(104)*UV(26)+JVS(108)*UV(27)+JVS(151)*UV(30)+JVS(202)*UV(37)+JVS(390)*UV(50)+JVS(489)&
               &*UV(54)+JVS(524)*UV(55)
  JTUV(27) = JVS(109)*UV(27)+JVS(391)*UV(50)+JVS(442)*UV(52)+JVS(490)*UV(54)
  JTUV(28) = JVS(119)*UV(28)+JVS(392)*UV(50)+JVS(525)*UV(55)
  JTUV(29) = JVS(63)*UV(19)+JVS(135)*UV(29)+JVS(344)*UV(48)+JVS(393)*UV(50)
  JTUV(30) = JVS(120)*UV(28)+JVS(152)*UV(30)+JVS(183)*UV(36)+JVS(203)*UV(37)+JVS(258)*UV(42)+JVS(274)*UV(43)+JVS(394)&
               &*UV(50)+JVS(422)*UV(51)+JVS(468)*UV(53)+JVS(526)*UV(55)
  JTUV(31) = JVS(136)*UV(29)+JVS(157)*UV(31)+JVS(275)*UV(43)+JVS(310)*UV(45)+JVS(363)*UV(49)+JVS(443)*UV(52)+JVS(527)&
               &*UV(55)+JVS(555)*UV(56)
  JTUV(32) = JVS(137)*UV(29)+JVS(162)*UV(32)+JVS(184)*UV(36)+JVS(311)*UV(45)+JVS(331)*UV(47)+JVS(364)*UV(49)+JVS(444)&
               &*UV(52)+JVS(528)*UV(55)+JVS(556)*UV(56)
  JTUV(33) = JVS(121)*UV(28)+JVS(167)*UV(33)+JVS(185)*UV(36)+JVS(218)*UV(38)+JVS(259)*UV(42)+JVS(276)*UV(43)+JVS(332)&
               &*UV(47)+JVS(445)*UV(52)+JVS(529)*UV(55)+JVS(557)*UV(56)
  JTUV(34) = JVS(8)*UV(4)+JVS(11)*UV(5)+JVS(38)*UV(14)+JVS(52)*UV(17)+JVS(64)*UV(19)+JVS(91)*UV(24)+JVS(122)*UV(28)&
               &+JVS(172)*UV(34)+JVS(186)*UV(36)+JVS(204)*UV(37)+JVS(237)*UV(40)+JVS(248)*UV(41)+JVS(260)*UV(42)+JVS(277)&
               &*UV(43)+JVS(294)*UV(44)+JVS(347)*UV(48)+JVS(398)*UV(50)+JVS(423)*UV(51)+JVS(469)*UV(53)+JVS(491)*UV(54)&
               &+JVS(530)*UV(55)
  JTUV(35) = JVS(12)*UV(5)+JVS(123)*UV(28)+JVS(158)*UV(31)+JVS(163)*UV(32)+JVS(176)*UV(35)+JVS(187)*UV(36)+JVS(205)&
               &*UV(37)+JVS(278)*UV(43)+JVS(312)*UV(45)+JVS(399)*UV(50)+JVS(424)*UV(51)+JVS(470)*UV(53)+JVS(492)*UV(54)&
               &+JVS(531)*UV(55)
  JTUV(36) = JVS(110)*UV(27)+JVS(124)*UV(28)+JVS(188)*UV(36)+JVS(400)*UV(50)+JVS(493)*UV(54)+JVS(532)*UV(55)
  JTUV(37) = JVS(206)*UV(37)+JVS(448)*UV(52)+JVS(494)*UV(54)+JVS(533)*UV(55)+JVS(559)*UV(56)
  JTUV(38) = JVS(219)*UV(38)+JVS(295)*UV(44)+JVS(319)*UV(46)+JVS(401)*UV(50)+JVS(471)*UV(53)
  JTUV(39) = JVS(13)*UV(5)+JVS(39)*UV(14)+JVS(53)*UV(17)+JVS(65)*UV(19)+JVS(92)*UV(24)+JVS(125)*UV(28)+JVS(207)*UV(37)&
               &+JVS(220)*UV(38)+JVS(231)*UV(39)+JVS(238)*UV(40)+JVS(249)*UV(41)+JVS(261)*UV(42)+JVS(279)*UV(43)+JVS(296)&
               &*UV(44)+JVS(320)*UV(46)+JVS(349)*UV(48)+JVS(402)*UV(50)+JVS(425)*UV(51)+JVS(472)*UV(53)+JVS(495)*UV(54)&
               &+JVS(534)*UV(55)
  JTUV(40) = JVS(54)*UV(17)+JVS(79)*UV(21)+JVS(239)*UV(40)+JVS(280)*UV(43)+JVS(450)*UV(52)+JVS(496)*UV(54)+JVS(535)&
               &*UV(55)+JVS(561)*UV(56)
  JTUV(41) = JVS(66)*UV(19)+JVS(190)*UV(36)+JVS(250)*UV(41)+JVS(281)*UV(43)+JVS(366)*UV(49)+JVS(451)*UV(52)+JVS(497)&
               &*UV(54)+JVS(536)*UV(55)+JVS(562)*UV(56)
  JTUV(42) = JVS(111)*UV(27)+JVS(126)*UV(28)+JVS(208)*UV(37)+JVS(262)*UV(42)+JVS(405)*UV(50)+JVS(473)*UV(53)+JVS(498)&
               &*UV(54)+JVS(537)*UV(55)
  JTUV(43) = JVS(112)*UV(27)+JVS(127)*UV(28)+JVS(282)*UV(43)+JVS(297)*UV(44)+JVS(406)*UV(50)+JVS(474)*UV(53)+JVS(499)&
               &*UV(54)+JVS(538)*UV(55)
  JTUV(44) = JVS(75)*UV(20)+JVS(93)*UV(24)+JVS(191)*UV(36)+JVS(298)*UV(44)+JVS(454)*UV(52)+JVS(500)*UV(54)+JVS(539)&
               &*UV(55)+JVS(564)*UV(56)
  JTUV(45) = JVS(14)*UV(5)+JVS(113)*UV(27)+JVS(128)*UV(28)+JVS(141)*UV(29)+JVS(168)*UV(33)+JVS(192)*UV(36)+JVS(209)&
               &*UV(37)+JVS(221)*UV(38)+JVS(263)*UV(42)+JVS(283)*UV(43)+JVS(299)*UV(44)+JVS(313)*UV(45)+JVS(367)*UV(49)&
               &+JVS(408)*UV(50)+JVS(426)*UV(51)+JVS(475)*UV(53)+JVS(501)*UV(54)+JVS(540)*UV(55)
  JTUV(46) = JVS(193)*UV(36)+JVS(222)*UV(38)+JVS(264)*UV(42)+JVS(322)*UV(46)+JVS(335)*UV(47)+JVS(456)*UV(52)+JVS(476)&
               &*UV(53)+JVS(502)*UV(54)+JVS(541)*UV(55)+JVS(566)*UV(56)
  JTUV(47) = JVS(67)*UV(19)+JVS(210)*UV(37)+JVS(223)*UV(38)+JVS(240)*UV(40)+JVS(265)*UV(42)+JVS(284)*UV(43)+JVS(336)&
               &*UV(47)+JVS(352)*UV(48)+JVS(410)*UV(50)+JVS(542)*UV(55)
  JTUV(48) = JVS(68)*UV(19)+JVS(211)*UV(37)+JVS(224)*UV(38)+JVS(241)*UV(40)+JVS(285)*UV(43)+JVS(337)*UV(47)+JVS(353)&
               &*UV(48)+JVS(368)*UV(49)+JVS(458)*UV(52)+JVS(504)*UV(54)+JVS(543)*UV(55)+JVS(568)*UV(56)
  JTUV(49) = JVS(69)*UV(19)+JVS(212)*UV(37)+JVS(225)*UV(38)+JVS(242)*UV(40)+JVS(251)*UV(41)+JVS(286)*UV(43)+JVS(369)&
               &*UV(49)+JVS(412)*UV(50)+JVS(459)*UV(52)+JVS(544)*UV(55)
  JTUV(50) = JVS(3)*UV(1)+JVS(21)*UV(7)+JVS(23)*UV(8)+JVS(25)*UV(9)+JVS(31)*UV(11)+JVS(36)*UV(13)+JVS(40)*UV(14)+JVS(44)&
               &*UV(15)+JVS(48)*UV(16)+JVS(55)*UV(17)+JVS(60)*UV(18)+JVS(70)*UV(19)+JVS(76)*UV(20)+JVS(80)*UV(21)+JVS(83)&
               &*UV(22)+JVS(88)*UV(23)+JVS(94)*UV(24)+JVS(105)*UV(26)+JVS(114)*UV(27)+JVS(129)*UV(28)+JVS(145)*UV(29)&
               &+JVS(153)*UV(30)+JVS(164)*UV(32)+JVS(169)*UV(33)+JVS(173)*UV(34)+JVS(177)*UV(35)+JVS(194)*UV(36)+JVS(213)&
               &*UV(37)+JVS(226)*UV(38)+JVS(232)*UV(39)+JVS(243)*UV(40)+JVS(252)*UV(41)+JVS(266)*UV(42)+JVS(287)*UV(43)&
               &+JVS(304)*UV(44)+JVS(314)*UV(45)+JVS(326)*UV(46)+JVS(339)*UV(47)+JVS(355)*UV(48)+JVS(370)*UV(49)+JVS(413)&
               &*UV(50)+JVS(427)*UV(51)+JVS(460)*UV(52)+JVS(480)*UV(53)+JVS(506)*UV(54)+JVS(545)*UV(55)+JVS(570)*UV(56)
  JTUV(51) = JVS(9)*UV(4)+JVS(15)*UV(5)+JVS(19)*UV(6)+JVS(41)*UV(14)+JVS(56)*UV(17)+JVS(71)*UV(19)+JVS(89)*UV(23)&
               &+JVS(95)*UV(24)+JVS(98)*UV(25)+JVS(130)*UV(28)+JVS(154)*UV(30)+JVS(174)*UV(34)+JVS(178)*UV(35)+JVS(195)&
               &*UV(36)+JVS(214)*UV(37)+JVS(227)*UV(38)+JVS(233)*UV(39)+JVS(244)*UV(40)+JVS(267)*UV(42)+JVS(288)*UV(43)&
               &+JVS(305)*UV(44)+JVS(315)*UV(45)+JVS(327)*UV(46)+JVS(356)*UV(48)+JVS(414)*UV(50)+JVS(428)*UV(51)+JVS(461)&
               &*UV(52)+JVS(481)*UV(53)+JVS(507)*UV(54)+JVS(546)*UV(55)+JVS(571)*UV(56)
  JTUV(52) = JVS(28)*UV(10)+JVS(33)*UV(12)+JVS(45)*UV(15)+JVS(49)*UV(16)+JVS(84)*UV(22)+JVS(99)*UV(25)+JVS(115)*UV(27)&
               &+JVS(372)*UV(49)+JVS(415)*UV(50)+JVS(429)*UV(51)+JVS(462)*UV(52)+JVS(482)*UV(53)+JVS(508)*UV(54)+JVS(547)&
               &*UV(55)+JVS(572)*UV(56)
  JTUV(53) = JVS(16)*UV(5)+JVS(29)*UV(10)+JVS(306)*UV(44)+JVS(430)*UV(51)+JVS(463)*UV(52)+JVS(483)*UV(53)+JVS(509)&
               &*UV(54)+JVS(548)*UV(55)+JVS(573)*UV(56)
  JTUV(54) = JVS(34)*UV(12)+JVS(46)*UV(15)+JVS(72)*UV(19)+JVS(100)*UV(25)+JVS(106)*UV(26)+JVS(116)*UV(27)+JVS(131)&
               &*UV(28)+JVS(147)*UV(29)+JVS(159)*UV(31)+JVS(175)*UV(34)+JVS(179)*UV(35)+JVS(196)*UV(36)+JVS(215)*UV(37)&
               &+JVS(228)*UV(38)+JVS(234)*UV(39)+JVS(245)*UV(40)+JVS(254)*UV(41)+JVS(268)*UV(42)+JVS(289)*UV(43)+JVS(307)&
               &*UV(44)+JVS(316)*UV(45)+JVS(328)*UV(46)+JVS(357)*UV(48)+JVS(373)*UV(49)+JVS(417)*UV(50)+JVS(464)*UV(52)&
               &+JVS(484)*UV(53)+JVS(510)*UV(54)+JVS(549)*UV(55)+JVS(574)*UV(56)
  JTUV(55) = JVS(17)*UV(5)+JVS(26)*UV(9)+JVS(50)*UV(16)+JVS(77)*UV(20)+JVS(81)*UV(21)+JVS(85)*UV(22)+JVS(117)*UV(27)&
               &+JVS(148)*UV(29)+JVS(160)*UV(31)+JVS(165)*UV(32)+JVS(170)*UV(33)+JVS(216)*UV(37)+JVS(246)*UV(40)+JVS(255)&
               &*UV(41)+JVS(308)*UV(44)+JVS(329)*UV(46)+JVS(342)*UV(47)+JVS(358)*UV(48)+JVS(374)*UV(49)+JVS(418)*UV(50)&
               &+JVS(432)*UV(51)+JVS(465)*UV(52)+JVS(485)*UV(53)+JVS(511)*UV(54)+JVS(550)*UV(55)+JVS(575)*UV(56)
  JTUV(56) = JVS(61)*UV(18)+JVS(73)*UV(19)+JVS(86)*UV(22)+JVS(101)*UV(25)+JVS(132)*UV(28)+JVS(149)*UV(29)+JVS(156)&
               &*UV(30)+JVS(161)*UV(31)+JVS(166)*UV(32)+JVS(171)*UV(33)+JVS(198)*UV(36)+JVS(217)*UV(37)+JVS(230)*UV(38)&
               &+JVS(247)*UV(40)+JVS(256)*UV(41)+JVS(270)*UV(42)+JVS(291)*UV(43)+JVS(309)*UV(44)+JVS(318)*UV(45)+JVS(330)&
               &*UV(46)+JVS(359)*UV(48)+JVS(375)*UV(49)+JVS(419)*UV(50)+JVS(433)*UV(51)+JVS(466)*UV(52)+JVS(486)*UV(53)&
               &+JVS(512)*UV(54)+JVS(551)*UV(55)+JVS(576)*UV(56)
END SUBROUTINE JacTR_SP_Vec
END MODULE cbmz_bb_Jacobian


text:
! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
!
! --- End License

!
!////////////////////////////////////////////////////////////////////////
!
!      Main.f90
!      Created: February 22, 2013 6:18 PM
!      By: NocturnalAviationSoftware
!!
!////////////////////////////////////////////////////////////////////////
!
      PROGRAM TestObjectsMain
      USE TestSuiteManagerClass
      IMPLICIT NONE

      TYPE(TestSuiteManager) :: testSuite
      INTEGER                :: numberOfFailedTests

      EXTERNAL :: FTDictionaryClassTests
      EXTERNAL :: FTExceptionClassTests
      EXTERNAL :: FTValueClassTests
      EXTERNAL :: FTValueDictionaryClassTests
      EXTERNAL :: FTLinkedListClassTests
      EXTERNAL :: StackClassTests
      EXTERNAL :: MutableArrayClassTests
      EXTERNAL :: SparseMatrixTests
      EXTERNAL :: MultiIndexTableTests
      EXTERNAL :: FTStringSetTests
      EXTERNAL :: OptionalDataTest

      CHARACTER(LEN=1), POINTER :: optData(:)

!      CALL setvbuf3f(6,2,0) !PGIFortran only
!
!     -----
!     Setup
!     -----
!
      CALL testSuite % init()

      CALL testSuite % addTestSubroutineWithName(FTValueClassTests,"FTValueClass Tests")
      CALL testSuite % addTestSubroutineWithName(FTDictionaryClassTests,"FTDictionaryClass Tests")
      CALL testSuite % addTestSubroutineWithName(FTValueDictionaryClassTests,"FTValueDictionaryClass Tests")
      CALL testSuite % addTestSubroutineWithName(FTLinkedListClassTests,"FTLinkedListClass Tests")
      CALL testSuite % addTestSubroutineWithName(StackClassTests,"StackClass Tests")
      CALL testSuite % addTestSubroutineWithName(MutableArrayClassTests,"Mutable Array Tests")
      CALL testSuite % addTestSubroutineWithName(FTExceptionClassTests,"FTExceptionClass Tests")
      CALL testSuite % addTestSubroutineWithName(SparseMatrixTests,"SparseMatrixClass Tests")
      CALL testSuite % addTestSubroutineWithName(MultiIndexTableTests,"MultiIndexTable Tests" )
      CALL testSuite % addTestSubroutineWithName(FTStringSetTests,"String set Tests" )

      ALLOCATE(optData(2))
      optData(1) = "a"
      optdata(2) = "b"
      CALL testSuite % addTestSubroutineWithName(OptionalDataTest,"Optional Data Tests", optData)

!
!     -------------
!     Run the tests
!     -------------
!
      CALL testSuite % performTests(numberOfFailedTests)

!
!     ---------------------------------------
!     Exit with error in case of failed tests
!     ---------------------------------------
!
      IF (numberOfFailedTests .gt. 0) THEN
        ERROR STOP 'At least one test has failed'
      END IF

      END PROGRAM TestObjectsMain

max_stars_repo_name: HigherOrderMethods/FTObjectLibrary

max_issues_repo_name: HigherOrderMethods/FTObjectLibrary

max_forks_repo_name: HigherOrderMethods/FTObjectLibrary


text:
      SUBROUTINE jacprod(c, h, nots, nb)
      USE vspline
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER nots, nb, jmax
      REAL(rprec), DIMENSION(*) :: c, h
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      REAL(rprec), DIMENSION(nots) ::
     1   aspline, bspline, dspline, dum1
C-----------------------------------------------

!
!       THIS ROUTINE COMPUTES THE INNER PRODUCT COUT(I) = CIN(J)*JACOBIAN(J,I)
!       WHERE JACOBIAN(J,I) = D[G(J)]/D[F(I)]
!       HERE, G(J) ARE SECOND-DERIVATIVE KNOTS, F(I) FUNCTION KNOTS
!
!       COMPUTE COEFFICIENT ARRAY ELEMENTS A*X(I+1) + D*X(I) + B*X(I-1)
!       (TO BE SAFE, RECOMPUTE EACH TIME, SINCE IOTA, P SPLINES MAY
!        DIFFER FROM CALL TO CALL)
!
      aspline(1) = h(1)
      dspline(1) = 2.0*h(1)
      aspline(2:nots-1) = h(2:nots-1)
      bspline(2:nots-1) = h(:nots-2)
      dspline(2:nots-1) = 2.0*(h(2:nots-1)+h(:nots-2))

      jspmin(1) = 2
      IF (nb .eq. ideriv) jspmin(1) = 1
      jmax = nots - 1
      CALL tridslv(aspline,dspline,bspline,c,jspmin,jmax,0,nots,1)
      dum1(1) = 6.0*(c(2)-c(1))/h(1)
      dum1(2:nots) = 6.0*(c(:nots-1)-c(2:nots))/h(:nots-1)
      c(2:nots-1) = dum1(2:nots-1) - dum1(3:nots)
      c(1) = dum1(1)
      c(nots) = dum1(nots)

      END SUBROUTINE jacprod


! { dg-do run }
! PR24224 Test formatted input/output to/from character arrays with strides
! other than 1.  Contributed by Jerry DeLisle <jvdelisle@verizon.net>.
      program arrayio_6
      implicit none
      integer         :: i(3),j,k(3)
      character(12)  :: r(4,4,4) = '0123456789AB'
      character(12)  :: s(64)
      equivalence(r,s)

      i = (/(j,j=1,3)/)
      write(r(1:4:2,2:4:1,3:4:2),'(3(2x,i4/)/3(3x,i6/))') i

      if (s(36).ne.'0123456789AB') STOP 1
      if (s(37).ne.'     1      ') STOP 2
      if (s(38).ne.'0123456789AB') STOP 3
      if (s(39).ne.'     2      ') STOP 4
      if (s(40).ne.'0123456789AB') STOP 5
      if (s(41).ne.'     3      ') STOP 6
      if (s(42).ne.'0123456789AB') STOP 7
      if (s(43).ne.'            ') STOP 8
      if (s(44).ne.'0123456789AB') STOP 9
      if (s(45).ne.'            ') STOP 10
      if (s(46).ne.'0123456789AB') STOP 11

      k = i
      i = 0
      read(r(1:4:2,2:4:1,3:4:2),'(3(2x,i4/)/3(3x,i6/))') i
      if (any(i.ne.k)) STOP 12

      end program arrayio_6


program TestBPWriteAttributes
    use small_test_data
    use mpi
    use adios2
    implicit none

    type(adios2_adios) :: adios
    type(adios2_io) :: ioWrite, ioRead
    type(adios2_engine) :: bpWriter, bpReader
    type(adios2_attribute), dimension(14) :: attributes, attributes_in

    integer :: ierr, i
    character(len=23):: iString_value
    integer(kind=1):: i8_value
    integer(kind=2):: i16_value
    integer(kind=4):: i32_value
    integer(kind=8):: i64_value
    real :: r32_value
    real(kind=8):: r64_value

    character(len=15), dimension(3):: iString_array
    integer(kind=1), dimension(3):: i8_array
    integer(kind=2), dimension(3):: i16_array
    integer(kind=4), dimension(3):: i32_array
    integer(kind=8), dimension(3):: i64_array
    real, dimension(3) :: r32_array
    real(kind=8), dimension(3):: r64_array

    ! Launch MPI
    call MPI_Init(ierr)

    ! Create adios handler passing the communicator, debug mode and error flag
    call adios2_init(adios, MPI_COMM_WORLD, adios2_debug_mode_on, ierr)

    !!!!!!!!!!!!!!!!!!!!!!!!! WRITER !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Declare an IO process configuration inside adios
    call adios2_declare_io(ioWrite, adios, "ioWrite", ierr)

    do i=1,14
        if( attributes(i)%valid .eqv. .true. ) stop 'Invalid attribute default'
    end do

    ! single value
    call adios2_define_attribute(attributes(1), ioWrite, 'att_String', &
                                 'ADIOS2 String attribute', ierr)

    call adios2_define_attribute(attributes(2), ioWrite, 'att_i8', &
                                 data_I8(1), ierr)

    call adios2_define_attribute(attributes(3), ioWrite, 'att_i16', &
                                 data_I16(1), ierr)

    call adios2_define_attribute(attributes(4), ioWrite, 'att_i32', &
                                 data_I32(1), ierr)

    call adios2_define_attribute(attributes(5), ioWrite, 'att_i64', &
                                 data_I64(1), ierr)

    call adios2_define_attribute(attributes(6), ioWrite, 'att_r32', &
                                 data_R32(1), ierr)

    call adios2_define_attribute(attributes(7), ioWrite, 'att_r64', &
                                 data_R64(1), ierr)

    ! arrays
    call adios2_define_attribute(attributes(8), ioWrite, 'att_Strings_array', &
                                 data_Strings, 3, ierr)

    call adios2_define_attribute(attributes(9), ioWrite, 'att_i8_array', &
                                 data_I8, 3, ierr)

    call adios2_define_attribute(attributes(10), ioWrite, 'att_i16_array', &
                                 data_I16, 3, ierr)

    call adios2_define_attribute(attributes(11), ioWrite, 'att_i32_array', &
                                 data_I32, 3, ierr)

    call adios2_define_attribute(attributes(12), ioWrite, 'att_i64_array', &
                                 data_I64, 3, ierr)

    call adios2_define_attribute(attributes(13), ioWrite, 'att_r32_array', &
                                 data_R32, 3, ierr)

    call adios2_define_attribute(attributes(14), ioWrite, 'att_r64_array', &
                                 data_R64, 3, ierr)

    do i=1,14
        if( attributes(i)%valid .eqv. .false. ) stop 'Invalid adios2_define_attribute'
    end do

    call adios2_open(bpWriter, ioWrite, "fattr_types.bp", adios2_mode_write, &
                     ierr)

    call adios2_close(bpWriter, ierr)

    call MPI_Barrier(MPI_COMM_WORLD, ierr)

    !!!!!!!!!!!!!!!!!!!!!!!!! READER !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    call adios2_declare_io(ioRead, adios, "ioRead", ierr)

    call adios2_open(bpReader, ioRead, 'fattr_types.bp', adios2_mode_read, ierr)

    call adios2_inquire_attribute(attributes_in(1), ioRead, 'att_String', ierr)
    call adios2_inquire_attribute(attributes_in(2), ioRead, 'att_i8', ierr)
    call adios2_inquire_attribute(attributes_in(3), ioRead, 'att_i16', ierr)
    call adios2_inquire_attribute(attributes_in(4), ioRead, 'att_i32', ierr)
    call adios2_inquire_attribute(attributes_in(5), ioRead, 'att_i64', ierr)
    call adios2_inquire_attribute(attributes_in(6), ioRead, 'att_r32', ierr)
    call adios2_inquire_attribute(attributes_in(7), ioRead, 'att_r64', ierr)

    if(attributes_in(1)%valid .eqv. .false.) stop 'attribute iString not found'
    if(attributes_in(1)%type /= adios2_type_string) stop 'attribute iString wrong type'
    if(attributes_in(1)%length /= 1) stop 'attribute iString length is not 1'
    if(attributes_in(1)%is_value .eqv. .false.) stop 'attribute iString must be value'
    call adios2_attribute_data( iString_value, attributes_in(1), ierr)
    if( iString_value /=  'ADIOS2 String attribute' ) stop 'attribute iString data error'

    if(attributes_in(2)%valid .eqv. .false.) stop 'attribute i8 not found'
    if(attributes_in(2)%type /= adios2_type_integer1) stop 'attribute i8 wrong type'
    if(attributes_in(2)%length /= 1) stop 'attribute i8 length is not 1'
    if(attributes_in(2)%is_value .eqv. .false.) stop 'attribute i8 must be value'
    call adios2_attribute_data( i8_value, attributes_in(2), ierr)
    if( i8_value /=  data_I8(1) ) stop 'attribute i8 data error'

    if(attributes_in(3)%valid .eqv. .false.) stop 'attribute i16 not found'
    if(attributes_in(3)%type /= adios2_type_integer2) stop 'attribute i16 wrong type'
    if(attributes_in(3)%length /= 1) stop 'attribute i16 length is not 1'
    if(attributes_in(3)%is_value .eqv. .false.) stop 'attribute i16 must be value'
    call adios2_attribute_data( i16_value, attributes_in(3), ierr)
    if( i16_value /=  data_I16(1) ) stop 'attribute i16 data error'

    if(attributes_in(4)%valid .eqv. .false.) stop 'attribute i32 not found'
    if(attributes_in(4)%type /= adios2_type_integer4) stop 'attribute i32 wrong type'
    if(attributes_in(4)%length /= 1) stop 'attribute i32 length is not 1'
    if(attributes_in(4)%is_value .eqv. .false.) stop 'attribute i32 must be value'
    call adios2_attribute_data( i32_value, attributes_in(4), ierr)
    if( i32_value /=  data_I32(1) ) stop 'attribute i32 data error'

    if(attributes_in(5)%valid .eqv. .false.) stop 'attribute i64 not found'
    if(attributes_in(5)%type /= adios2_type_integer8) stop 'attribute i64 wrong type'
    if(attributes_in(5)%length /= 1) stop 'attribute i64 length is not 1'
    if(attributes_in(5)%is_value .eqv. .false.) stop 'attribute i64 must be value'
    call adios2_attribute_data( i64_value, attributes_in(5), ierr)
    if( i64_value /=  data_I64(1) ) stop 'attribute i64 data error'

    if(attributes_in(6)%valid .eqv. .false.) stop 'attribute r32 not found'
    if(attributes_in(6)%type /= adios2_type_real) stop 'attribute r32 wrong type'
    if(attributes_in(6)%length /= 1) stop 'attribute r32 length is not 1'
    if(attributes_in(6)%is_value .eqv. .false.) stop 'attribute r32 must be value'
    call adios2_attribute_data( r32_value, attributes_in(6), ierr)
    if( r32_value /=  data_R32(1) ) stop 'attribute r32 data error'

    if(attributes_in(7)%valid .eqv. .false.) stop 'attribute r64 not found'
    if(attributes_in(7)%type /= adios2_type_dp) stop 'attribute r64 wrong type'
    if(attributes_in(7)%length /= 1) stop 'attribute r64 length is not 1'
    if(attributes_in(7)%is_value .eqv. .false.) stop 'attribute r64 must be value'
    call adios2_attribute_data( r64_value, attributes_in(7), ierr)
    if( r64_value /=  data_R64(1) ) stop 'attribute r64 data error'

    ! Array
    call adios2_inquire_attribute(attributes_in(8), ioRead, 'att_Strings_array', ierr)
    call adios2_inquire_attribute(attributes_in(9), ioRead, 'att_i8_array', ierr)
    call adios2_inquire_attribute(attributes_in(10), ioRead, 'att_i16_array', ierr)
    call adios2_inquire_attribute(attributes_in(11), ioRead, 'att_i32_array', ierr)
    call adios2_inquire_attribute(attributes_in(12), ioRead, 'att_i64_array', ierr)
    call adios2_inquire_attribute(attributes_in(13), ioRead, 'att_r32_array', ierr)
    call adios2_inquire_attribute(attributes_in(14), ioRead, 'att_r64_array', ierr)

    if(attributes_in(8)%valid .eqv. .false.) stop 'attribute string array not found'
    if(attributes_in(8)%type /= adios2_type_string) stop 'attribute string array wrong type'
    if(attributes_in(8)%length /= 3) stop 'attribute string array length is not 3'
    if(attributes_in(8)%is_value .eqv. .true.) stop 'attribute string array must be array'
    call adios2_attribute_data( iString_array, attributes_in(8), ierr)
    do i=1,3
       if( iString_array(i) /= data_Strings(i) ) stop 'attribute string array data error'
    end do

    if(attributes_in(9)%valid .eqv. .false.) stop 'attribute i8 array not found'
    if(attributes_in(9)%type /= adios2_type_integer1) stop 'attribute i8 array wrong type'
    if(attributes_in(9)%length /= 3) stop 'attribute i8 array length is not 3'
    if(attributes_in(9)%is_value .eqv. .true.) stop 'attribute i8 array must be array'
    call adios2_attribute_data( i8_array, attributes_in(9), ierr)
    do i=1,3
       if( i8_array(i) /=  data_I8(i) ) stop 'attribute i8 array data error'
    end do

    if(attributes_in(10)%valid .eqv. .false.) stop 'attribute i16 array not found'
    if(attributes_in(10)%type /= adios2_type_integer2) stop 'attribute i16 array wrong type'
    if(attributes_in(10)%length /= 3) stop 'attribute i16 array length is not 3'
    if(attributes_in(10)%is_value .eqv. .true.) stop 'attribute i16 array must be array'
    call adios2_attribute_data( i16_array, attributes_in(10), ierr)
    do i=1,3
       if( i16_array(i) /=  data_I16(i) ) stop 'attribute i16 array data error'
    end do

    if(attributes_in(11)%valid .eqv. .false.) stop 'attribute i32 array not found'
    if(attributes_in(11)%type /= adios2_type_integer4) stop 'attribute i32 array wrong type'
    if(attributes_in(11)%length /= 3) stop 'attribute i32 array length is not 3'
    if(attributes_in(11)%is_value .eqv. .true.) stop 'attribute i32 array must be array'
    call adios2_attribute_data( i32_array, attributes_in(11), ierr)
    do i=1,3
       if( i32_array(i) /=  data_I32(i) ) stop 'attribute i32 array data error'
    end do

    if(attributes_in(12)%valid .eqv. .false.) stop 'attribute i64 array not found'
    if(attributes_in(12)%type /= adios2_type_integer8) stop 'attribute i64 array wrong type'
    if(attributes_in(12)%length /= 3) stop 'attribute i64 array length is not 3'
    if(attributes_in(12)%is_value .eqv. .true.) stop 'attribute i64 array must be array'
    call adios2_attribute_data( i64_array, attributes_in(12), ierr)
    do i=1,3
       if( i64_array(i) /=  data_I64(i) ) stop 'attribute i64 array data error'
    end do

    if(attributes_in(13)%valid .eqv. .false.) stop 'attribute r32 array not found'
    if(attributes_in(13)%type /= adios2_type_real) stop 'attribute r32 array wrong type'
    if(attributes_in(13)%length /= 3) stop 'attribute r32 array length is not 3'
    if(attributes_in(13)%is_value .eqv. .true.) stop 'attribute r32 array must be array'
    call adios2_attribute_data( r32_array, attributes_in(13), ierr)
    do i=1,3
       if( r32_array(i) /=  data_R32(i) ) stop 'attribute r32 array data error'
    end do

    if(attributes_in(14)%valid .eqv. .false.) stop 'attribute r64 array not found'
    if(attributes_in(14)%type /= adios2_type_dp) stop 'attribute r64 array wrong type'
    if(attributes_in(14)%length /= 3) stop 'attribute r64 array length is not 3'
    if(attributes_in(14)%is_value .eqv. .true.) stop 'attribute r64 array must be array'
    call adios2_attribute_data( r64_array, attributes_in(14), ierr)
    do i=1,3
       if( r64_array(i) /=  data_R64(i) ) stop 'attribute r64 array data error'
    end do

    call adios2_close(bpReader, ierr)

    call adios2_finalize(adios, ierr)

    call MPI_Finalize(ierr)

end program TestBPWriteAttributes


text:
! (C) Copyright 1996- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation nor
! does it submit to any jurisdiction.

!> @brief Variables needed for the mark5 calculation
!> @author Di Giuseppe, F., ECMWF
!> @author Maciel, P., ECMWF
MODULE mo_mark5

    USE mo_constants

    IMPLICIT NONE

    ! intermediate (fuel) variables
    TYPE mark5_fuel_type
        REAL :: moist             = rfillvalue  !< fuel moisture content
        REAL :: weight            = rfillvalue  !< fuel weight
        REAL :: curing            = rfillvalue  !< degree of curing
        REAL :: kb_drought_index  = rfillvalue  !< Keetch-Byram drought index (0 no drought - 203 severe drought)
        REAL :: drought_factor    = rfillvalue  !< drought factor
        REAL :: timesincerain     = rfillvalue
    END TYPE

    ! Mark 5 Fuel properties and behaviours
    TYPE mark5_prop_type
        REAL :: ros_theta0        = rfillvalue  !< spread component
        REAL :: ros_theta         = rfillvalue  !< spread rate
        REAL :: flame_height      = rfillvalue
        REAL :: flame_distance    = rfillvalue  !< distance at which flame can be seen
    END TYPE
    ! Mark 5 Fire danger
    TYPE mark5_prob_type
        REAL :: fire_danger_index = rfillvalue
    END TYPE

    TYPE(mark5_fuel_type), ALLOCATABLE :: mark5_fuel(:)
    TYPE(mark5_prob_type), ALLOCATABLE :: mark5_prob(:)
    TYPE(mark5_prop_type), ALLOCATABLE :: mark5_prop(:)

END MODULE mo_mark5

max_stars_repo_name: ecmwf-projects/geff

max_issues_repo_name: ecmwf-projects/geff

max_forks_repo_name: ecmwf-projects/geff


text: +CHECKER=fwrisc_exec_formal_smoke_checker


text:
      subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
      integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
      double precision t, y, pd, alph1, alph2
      dimension y(*), pd(nrowpd,*)
      dimension neq(*)
      data alph1/1.0d0/, alph2/1.0d0/, ng/5/
      mband = ml + mu + 1
      mu1 = mu + 1
      mu2 = mu + 2
      do 10 j = 1,neq(1)
        pd(mu1,j) = -2.0d0
        pd(mu2,j) = alph1
 10     pd(mband,j) = alph2
      do 20 j = ng,neq(1),ng
 20     pd(mu2,j) = 0.0d0
      return
      end

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


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

subroutine calcAext
!$ use omp_lib
use salmon_parallel, only: nproc_id_global
use salmon_communication, only: comm_is_root, comm_summation
use scf_data
implicit none
integer :: j2
integer :: ipulse
real(8) :: env_trigon_1
real(8) :: env_trigon_2

select case(ikind_eext)
  case(0)
    A_ext(1:3,:)=0.d0
    A_ext(1,:)=Fst*epdir_re1(1)
    A_ext(2,:)=Fst*epdir_re1(2)
    A_ext(3,:)=Fst*epdir_re1(3)
  case(1,15)
    A_ext(1:3,:)=0.d0
    do itt=Miter_rt+1,itotNtime+1
      if(dt*dble(itt) <= pulse_tw1)then
        ipulse=1
        call calc_env_trigon(ipulse,env_trigon_1)
        do j2=1,3
          A_ext(j2,itt)=A_ext(j2,itt)+amplitude1*epdir_re1(j2)*env_trigon_1
        end do
      end if
      if(dt*dble(itt)-t1_t2 >= 1.d-12 .and. dt*dble(itt)-t1_t2 <= pulse_tw2)then
        ipulse=2
        call calc_env_trigon(ipulse,env_trigon_2)
        do j2=1,3
          A_ext(j2,itt)=A_ext(j2,itt)+amplitude1*epdir_re1(j2)*env_trigon_2
        enddo
      end if
    end do
end select

return

end subroutine calcAext

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON