!********************************
!***Quality Assurance of Method**
!********************************
!	INCLUDE "c:\Program Files (x86)\Intel\oneAPI\mkl\2023.2.0\include\MKL_VSL.F90"

	SUBROUTINE QASIMULATIONS()
!!	USE MKL_VSL_TYPE
!	USE MKL_VSL
	USE WINTERACTER
	USE RESOURCE
	USE IONTOXMODEL_MODELVARS
	USE IONTOXMODEL_INTERFACEVARS
 	INTEGER NSIM/1000/,ISIM,JSIM,ICOUNT(4),IDAT,IPARM
	REAL*8 PARMTRUE(4),PARMSIM(1000,4),PARMMEAN(4),P,Z,RAND,MODELCALC,TEMP,NORMPZ
	REAL*8 ZSUM,ZSUM2,ZN
	EXTERNAL NORMPZ,MODELCALC,COMPLEXOPTIMIZE,CONFLIMITS,PARMESTIMATE,RAND
	ZN=0
	ZSUM=0
	ZSUM2=0
!***
!***Set "true" parameter values and parameter guesses to parameter solutions from basic analysis
!***
	DO IPARM=1,4
		PARMTRUE(IPARM)=PARMSOLUT(IPARM)
		PARMGUESS(IPARM)=PARMSOLUT(IPARM)
		PARMMEAN(IPARM)=0
		ICOUNT(IPARM)=0
		TYPE "('***QASIM*** PARMTRUE',I5,3F10.5)",IPARM,PARMTRUE(IPARM)
	ENDDO
!***
!***Start Simulation Loop
!***
	P=RAND(538595)
	DO ISIM=1,NSIM
		TYPE "('***QASIM*** ISIM',I5)",ISIM
!***
!***For same concentrations (XDAT), generate random numbers of YDAT
!***
		DO IPARM=1,4
			PARMTEMP(IPARM)=PARMTRUE(IPARM)
		ENDDO
        	DO IDAT=1,NDAT
			P=RAND(0)
			Z=NORMPZ(P)
			ZN=ZN+1
			ZSUM=ZSUM+Z
			ZSUM2=ZSUM2+Z**2
			YDAT(IDAT)=MODELCALC(XDAT(IDAT))+Z*PARMTRUE(4)
			TYPE "('***QASIM*** SIMULATED DATA',I5,4F10.5)",IDAT,P,Z,XDAT(IDAT),YDAT(IDAT)
		ENDDO
!***
!***Call analysis routines
!***
		IPHASE=3
		NPOINT=10
		PRINTFLAG=.FALSE.
		CALL COMPLEXOPTIMIZE(NPOINT,NPARM,PARMGUESS,PARMDELTA,PARMMIN,PARMMAX,PARMSOLUT,FUNCSOLUT,MODELERROR,NITER,NHIST,HISTITER,HISTFUNC,HISTPARM,PRINTFLAG)
		IPHASE=4
		CALL CONFLIMITS()
		DO IPARM=1,4
			PARMSIM(ISIM,IPARM)=PARMSOLUT(IPARM)
			TYPE "('***QASIM*** PARMSOLUT AND CLs',I5,3F10.5)",IPARM,PARMSOLUT(IPARM),PARMLCL(IPARM),PARMUCL(IPARM)
		ENDDO
!***
!***Evaluate whether parameter solution is within confidence limits and increment counters
!***
        	DO IPARM=1,4
        		IF(PARMTRUE(IPARM).GE.PARMLCL(IPARM).AND.PARMTRUE(IPARM).LE.PARMUCL(IPARM)) ICOUNT(IPARM)=ICOUNT(IPARM)+1
        	ENDDO
		TYPE "('***QASIM*** CONFLIM COUNTS',6I5)",ICOUNT
!***
!***End simulation loop and output QA results
!***
	TYPE "('STATS FOR Z',f6.0,2f10.3,2f10.5)",ZN,ZSUM,ZSUM2,ZSUM/ZN,(ZSUM2-(ZSUM/ZN)**2)/ZN
	ENDDO
	DO ISIM=1,NSIM-1
		DO JSIM=ISIM+1,NSIM
			DO IPARM=1,4
				IF(PARMSIM(JSIM,IPARM).LT.PARMSIM(ISIM,IPARM)) THEN
					TEMP=PARMSIM(JSIM,IPARM)
					PARMSIM(JSIM,IPARM)=PARMSIM(ISIM,IPARM)
					PARMSIM(ISIM,IPARM)=TEMP
				ENDIF
			ENDDO
		ENDDO
	ENDDO
	DO ISIM=1,NSIM
		TYPE "('***QASIM*** ORDERED PARMSET',I5,6F10.5)",ISIM,(PARMSIM(ISIM,IPARM),IPARM=1,4)
		DO IPARM=1,4
			PARMMEAN(IPARM)=PARMMEAN(IPARM)+PARMSIM(ISIM,IPARM)/NSIM
		ENDDO
	ENDDO
	DO IPARM=1,4
		CALL WGridPutCellDouble(IDF_MAIN_PARAMETERS_GRID,19,IPARM,1D2*ICOUNT(IPARM)/NSIM)
		CALL WGridPutCellDouble(IDF_MAIN_PARAMETERS_GRID,20,IPARM,PARMMEAN(IPARM))
	ENDDO
!***
!***Because this QA simulation uses the same analysis routines as the original analysis, variables that would
!***be needed in the output of the original analyses would be altered. Rather than create separate analysis   
!***routines for the simulations or save dozens of variables from the original analyses for the output routines, 
!***the original analysis is simply redone here, using the the initial guesses/limits still in the parameter grid.  
!***  
	CALL PARMESTIMATE()
	CALL CONFLIMITS()
	RETURN
	END SUBROUTINE QASIMULATIONS 

	FUNCTION NORMPZ(P)
	REAL*8 T,B1,B2,B3,B4,B5,C,PI
      REAL*8 NORMPZ,P,PTARG,PP(3),ZZ(3)
	INTEGER IPT,IT
	DATA B1,B2,B3,B4,B5,C,PI/0.319381530D0,-0.356563782D0,1.781477937D0,&
	-1.821255978D0,1.330274429D0,0.2316419D0,0.3989422804D0/
	IF (P.LT.0.0000003) THEN
		NORMPZ=-5.0
		RETURN
	ELSEIF (P.GT.0.9999997) THEN
		NORMPZ=5.0
		RETURN
      ELSEIF(P.GT.0.5D0) THEN
		PTARG=P
	ELSE
		PTARG=1-P
	ENDIF
	ZZ(1)=0.0D0
	ZZ(2)=1.0D0
	ZZ(3)=5.0D0
	DO IPT=1,3      
		T=1.0D0/(1.0D0+C*ZZ(IPT))
		PP(IPT)=1.0D0-PI*EXP(-ZZ(IPT)**2/2.0D0)*T*(B1+T*(B2+T*(B3+T*(B4+T*B5))))
	ENDDO
	DO IT=1,20
            IF(PP(2).GT.PTARG) THEN
			ZZ(3)=ZZ(2)
			PP(3)=PP(2)
		ELSE
			ZZ(1)=ZZ(2)
			PP(1)=PP(2)
		ENDIF
		ZZ(2)=(ZZ(1)+ZZ(3))/2
		T=1.0D0/(1.0D0+C*ZZ(2))
		PP(2)=1.0D0-PI*EXP(-ZZ(2)**2/2.0D0)*T*(B1+T*(B2+T*(B3+T*(B4+T*B5))))
	ENDDO			
	IF(P.GT.0.5D0) THEN
		NORMPZ=ZZ(2)
	ELSE
		NORMPZ=-ZZ(2)
	ENDIF
   	RETURN
	END

