module write_mod
      implicit none
!     unit conversion
      double precision ,parameter :: h2icm = 219474.69d0
      double precision, parameter :: au2Debye = 2.541746d0
      character(len=250), parameter :: sep_line = '(250("-"))'
      character(len=250), parameter :: block_line = '(250("="))'

      contains

!     <Subroutine for writing the Output
      subroutine write_output
     >     (q,x1,x2,y,wt,par,p_act,p_spread,nset,npar,
     >     flag,lauf)
      use adia_mod, only: adia
      use dim_parameter,only: qn,ntot,numdatpt,ndiab
      use ctrans_mod,only: ctrans
      implicit none
!     IN: variables
      integer lauf
      integer flag              !< 0= initial output 1=fit not converged 2= Fit Converged, 3= max iteration reached
      integer npar,nset
      double precision par(npar,nset),p_spread(npar)
      integer p_act(npar)
      double precision q(qn,numdatpt),x1(qn,numdatpt),x2(qn,numdatpt)
      double precision y(ntot,numdatpt),wt(ntot,numdatpt)

!     INTERNAL: Variables
      integer,parameter  :: id_out = 20 , std_out = 6
      integer pt
      integer i, id_print
      double precision, allocatable ::  ymod(:,:)
      double precision, allocatable ::  ew(:,:)
      double precision, allocatable ::  ev(:,:,:)

      logical skip

      allocate(ymod(ntot,numdatpt))
      allocate(ew(ndiab,numdatpt))
      allocate(ev(ndiab,ndiab,numdatpt))

      skip=.false.

!     get Model Outputs for all geometries for current best parameter set par(:,1)
      do pt=1,numdatpt
         call adia(pt,par(1:npar,1),npar,ymod(1:ntot,pt),
     >        ew(1:ndiab,pt),ev(1:ndiab,1:ndiab,pt),skip)
         call ctrans(q(:,pt),x1(:,pt),x2(:,pt))
      enddo

!     Initial write print everything you want to see before the fit and return
      if(flag.eq.0) then
         call print_parameterstate(std_out,par(:,1),p_act,npar)
         call print_ErrorSummary(std_out,y,ymod,wt)
!        print Data into the plotfiles
         return
      endif
!     open output files for individual makro iterations
      call open_outfile(id_out,lauf)
!     print Data into the plotfiles
      call print_plotfiles(x1,y,wt,ymod)

!     print Genetic output into files
      do i=1, 2
         if (i.eq.1) then
            id_print= std_out
         else
            id_print= id_out
         endif
         write(id_print,'("Writing Iteration: ",i4)') lauf
         write(id_print,block_line)
!        write data information only in outfile
         if(i.eq.2) then
            call print_data(id_print,x1,y,ymod,wt)
            call print_Set_Errors(id_print,y,ymod,wt)
         endif
         call print_parameterblock
     >        (id_print,par(:,1),p_act,p_spread,npar)
         call print_ErrorSummary(id_print,y,ymod,wt)

      enddo

      call print_fortranfile(par(:,1),npar)
      
      ! write the type of calc at the end of the output


      close (id_out)
      deallocate(ymod,ev,ew)
      end subroutine
!----------------------------------------------------------------------------------------------------
!     <subroutine for scan seperated Error analysis>
      subroutine print_Set_Errors(id_out,y, ymod, wt)
      use io_parameters,only: llen
      use dim_parameter,only: ndata,nstat,ntot,numdatpt,sets
      integer , intent(in) :: id_out
      double precision, intent(in) ::  y(ntot,numdatpt),
     >     ymod(ntot,numdatpt), wt(ntot,numdatpt)
      integer :: set, setpoint, pt
      double precision :: Set_rms(sets,ntot), Set_num(sets,ntot)
      double precision :: Total_rms, Total_Energy_rms,Energy_rms(nstat)
      character(len=llen) fmt
      write(id_out,'(A)') 'Errors in icm for individual Sets' //
     >     '(specified by sets: and npoints:)'
      write(id_out,'(A5,3A16)')'Set','Total',
     >     'Total_Energy', 'Energy[nstat]'
      write(id_out,sep_line)
      write(fmt,'("(I5,2f16.1,",I2,"f16.1)")') nstat
      Set_rms = 0.d0
      pt = 0
      do set=1, sets
         do setpoint=1, ndata(set)
            pt = pt + 1
            where(wt(:,pt) > 0.d0)
               Set_rms(set,:) = Set_rms(set,:)+(ymod(:,pt)-y(:,pt))**2
               Set_num(set,:) = Set_num(set,:) + 1
            end where
         enddo
         Total_rms
     >        = dsqrt(sum(Set_rms(set,:))
     >        / (sum(Set_num(set,:))))
         Total_Energy_rms
     >        = dsqrt(sum(Set_rms(set,1:nstat))
     >        / (sum(Set_num(set,1:nstat))))
         Energy_rms(1:nstat)
     >        = dsqrt(Set_rms(set,1:nstat)
     >        / (Set_num(set,1:nstat)))
         write(id_out,fmt) set, Total_rms*h2icm, Total_Energy_rms*h2icm,
     >        Energy_rms(1:nstat)*h2icm
      enddo
      write(id_out,block_line)
      write(id_out,*) ''
      end subroutine print_Set_Errors

!----------------------------------------------------------------------------------------------------
!     <subroutine for printing the parameter and the pst vector in fortran readable style for including the fitted parameters in other programs
      subroutine print_fortranfile(p,npar)
      use io_parameters,only: maxpar_keys
      use dim_parameter,only: pst
      implicit none
!     IN: variables
      integer npar
      double precision p(npar)
!     INTERNAL: variables
      integer i
      integer, parameter :: id_out = 49
      character(len=32), parameter :: fname ='fit_genric_bend_no3.f90'

      open(id_out,file=fname)

 30   format(6x,A2,i3,A2,d18.9)
 31   format(6x,A6,i3,A2,i3)
      
      write(id_out,'(2X,A)') "Module dip_param"
      write(id_out,'(5X,A)') "IMPLICIT NONE"
      write(id_out,'(5X,A,I0)') "Integer,parameter :: np=",npar 
      write(id_out,'(5X,A,I0,A)') "Double precision :: p(",npar,")"
      write(id_out,'(5X,A,I0,A)') "integer :: pst(2,",maxpar_keys,")"
      write(id_out,'(5X,A)') "contains"
      write(id_out,*)''
  
      write (id_out,'(5x,a)') "SUBROUTINE init_dip_planar_data()"
      write (id_out,'(8X,A)') "implicit none"
      do i=1,npar
         write(id_out,30) 'p(',i,')=',p(i)
      enddo
      do i=1,maxpar_keys
         write(id_out,31) 'pst(1,',i,')=',pst(1,i)
         write(id_out,31) 'pst(2,',i,')=',pst(2,i)
      enddo
      

      write(id_out,"(A)") "End SUBROUTINE init_dip_planar_data"
      write(id_out,"(A)") "End Module dip_param"

      close(id_out)
      end subroutine

!----------------------------------------------------------------------------------------------------
!     <subroutine print_ErrorSummary: calculates the rms errros and prints them in the corresponding file
      subroutine print_ErrorSummary(id_out,y,ymod,wt)
      use dim_parameter,only: nstat,rms_thr,ntot,numdatpt
      use io_parameters,only: llen
      implicit none
!     IN: variables
      integer id_out
      double precision y(ntot,numdatpt),ymod(ntot,numdatpt)
      double precision wt(ntot,numdatpt)
!     INTERNAL: variables
!     Counter and RMS variables
      double precision  Cut_thr(nstat)
      double precision  Output_rms(ntot),Cut_rms(nstat),Weighted_rms
      integer Output_num(ntot),Cut_num(nstat)
      double precision  Weighted_wt
      double precision  Total_rms,Total_Weighted_rms
      double precision  Total_Energie_rms,Total_State_rms(nstat)
      double precision  Cut_Energie_rms, Cut_State_rms(nstat)

      ! Variables for computing the NRMSE 
      !double precision:: ymean(ntot),ysum(ntot),NRMSE 

!     loop control
      integer j,pt

!     Fabian
      character(len=llen) fmt
!     initialize RMS variables
      Output_rms(1:ntot) = 0.d0
      Output_num(1:ntot) = 0
      Weighted_rms  = 0.d0
      Weighted_wt = 0.d0
      Cut_rms(1:nstat)= 0.d0
      Cut_num(1:nstat)= 0

!     Define Threshold for Cut_* RMS Values
      Cut_thr(1:nstat) = rms_thr(1:nstat)
!     SUMM!
!     Loop over all Datapoints
      do pt=1,numdatpt
!        get unweighted rms for each output value and count their number
         do j=1,ntot
            if(wt(j,pt).gt.0.d0) then
               Output_rms(j) = Output_rms(j) +
     >              (ymod(j,pt)-y(j,pt))**2
               Output_num(j)=Output_num(j) + 1
            endif
         enddo
!        get the unweighted rms under the given threshold and count their number
         do j=1,nstat
            if(wt(j,pt).gt.0.d0) then
               if(y(j,pt).le.Cut_thr(j)) then
                  Cut_rms(j) = Cut_rms(j) +
     >                 (ymod(j,pt)-y(j,pt))**2
                  Cut_num(j) = Cut_num(j) + 1
               endif
            endif
         enddo
!        get the weighted rms over all output values
         Weighted_rms = Weighted_rms +
     >        sum(((ymod(1:ntot,pt)-y(1:ntot,pt))**2)
     >        *(wt(1:ntot,pt)**2))
         Weighted_wt = Weighted_wt + sum(wt(1:ntot,pt)**2)
      enddo

!     NORM!
!     TOTAL RMS:
!     unweighted
      Total_rms =
     >     dsqrt(sum(Output_rms(1:ntot)) /(sum(Output_num(1:ntot))))

!     Weighted
      Total_Weighted_rms = dsqrt(Weighted_rms/Weighted_wt)

!     unweighted, considering only first nstat values
      Total_Energie_rms =
     >     dsqrt(sum(Output_rms(1:nstat)) /(sum(Output_num(1:nstat))))

!     unweighted,for each of the first nstat values separatly
      Total_State_rms(1:nstat) =
     >     dsqrt(Output_rms(1:nstat) / Output_num(1:nstat))

!     unweighted,first nstat values only counting points under given threshold
      Cut_Energie_rms =
     >     dsqrt(sum(Cut_rms(1:nstat)) /(sum(Cut_num(1:nstat))))

!     unweighted,each nstat values seperatly only counting points under threshold
      Cut_State_rms(1:nstat) =
     >     dsqrt(Cut_rms(1:nstat)/Cut_num(1:nstat))

!     WRITE!
!     make the actual writing into the file
      write(id_out,39)
      write(id_out,40)
      write(id_out,41) Total_rms, Total_rms*au2Debye!Total_rms*h2icm
      write(id_out,42) sum(Output_num(1:ntot))
      write(id_out,43) Total_Weighted_rms, Total_Weighted_rms*h2icm
      write(id_out,44) Weighted_wt
      write(id_out,45) Total_Energie_rms, Total_Energie_rms*h2icm
      write(id_out,42) sum(Output_num(1:nstat))
      write(fmt,'("(A,10x,A,",I2,"f8.1)")') nstat
      write(id_out,fmt) '#','State resolved RMS(icm): ',
     $     Total_State_rms(1:nstat)*h2icm
      write(fmt,'("(A,10x,A,",I2,"i8)")') nstat
      write(id_out,fmt) '#','No. of Points per State: ',
     $     Output_num(1:nstat)
      write(id_out,51)

!     write the errors under a given threshold if there were any points
      if(any(Cut_num(1:nstat).gt.0)) then
         write(id_out,48) Cut_Energie_rms, Cut_Energie_rms*h2icm
         write(id_out,42) sum(Cut_num(1:nstat))

         write(fmt,'("(A,10x,A,",I2,"f8.1,A)")') nstat
         write(id_out,fmt) '#','Red. State resolved RMS: ',
     $        Cut_State_rms(1:nstat)*h2icm,' icm'
         write(fmt,'("(A,10x,A,",I2,"i8)")') nstat
         write(id_out,fmt) '#','No. of Points per State: ',
     $        Cut_num(1:nstat)
         write(fmt,'("(A,10x,A,",I2,"f8.1,A)")') nstat
         write(id_out,fmt) '#','Threshold per State:     ',
     $        Cut_thr(1:nstat)*h2icm,' icm above Reference Point.'         

      endif
      write(id_out,39)

!     FORMAT! specifications for the writing
 39   format(250('#'))
 40   format('#',10x,'ERROR SUMMARY: ')
 41   format('#',10x,'Total RMS:            ',g16.8, '(',g16.8,
     >     ' Debye)')
 42   format('#',10x,'No. of Points:        ',i10)
 43   format('#',10x,'Total weighted RMS:   ',g16.8, '(',f8.1,' icm)')
 44   format('#',10x,'Sum of point weights: ',f16.8)
 45   format('#',10x,'Total Energie RMS:    ',g16.8, '(',f8.1,' icm)')

 48   format('#',10x,'Red. Energie RMS:        ',g16.8,'(',f8.1,' icm)')
 51   format('#')

      end subroutine

!----------------------------------------------------------------------------------------------------
      subroutine print_plotfiles(x,y,wt,ymod)
      use dim_parameter,only: ndata,sets,qn,ntot,numdatpt,plot_coord
      implicit none
!     IN: variables
      double precision x(qn,numdatpt),y(ntot,numdatpt)
      double precision wt(ntot,numdatpt), ymod(ntot,numdatpt)
!     INTERNAL: variables
      integer sstart,ssend,set,id_plot

!     Initialize position pointer
      ssend=0
!     loop over datasets and print the plotfiles
      do set=1 ,sets
         if(ndata(set).eq.0) cycle
         id_plot=50+set
         call open_plotfile(id_plot,set)
         write(id_plot,'(A)') '# -*- truncate-lines: t -*-'
!        get start and end point of each set
         sstart=ssend+1
         ssend=ssend+ndata(set)
         if (plot_coord(set).eq.0) then
            call print_plotwalk(x(:,sstart:ssend),y(:,sstart:ssend),
     >        wt(:,sstart:ssend),ymod(:,sstart:ssend),
     >           ndata(set),id_plot,set)
         else
            call print_plotcoord(plot_coord(set),
     >           x(:,sstart:ssend),y(:,sstart:ssend),
     >           wt(:,sstart:ssend),ymod(:,sstart:ssend),
     >           ndata(set),id_plot,set)
         endif
         close(id_plot)
      enddo

      end subroutine

!----------------------------------------------------------------------------------------------------
      subroutine print_plotwalk(x,y,wt,ymod,npt,id_plot,set)
      use dim_parameter,only: qn,ntot
      use io_parameters,only: llen
      implicit none
!     IN: variables
      integer id_plot,npt,set
      double precision x(qn,npt),y(ntot,npt),ymod(ntot,npt),wt(ntot,npt)
!     INTERNAL: variables
      double precision xdiff(qn),walktime
      double precision walknorm
!     loop control
      integer i,j

      character(len=llen) fmt
      j=ntot-1

      call print_plotheader(id_plot,0,npt,set)

      call getwalknorm(x,walknorm,npt)
      walktime = 0.d0
      do i=1,npt
         if(i.gt.1) then
            xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1)
            walktime = walktime + dsqrt(sum(xdiff(1:qn)**2))/walknorm
         endif
         write(id_plot,"(ES16.8,*(3(ES16.8),:))")
     >        walktime ,ymod(:,i),y(:,i),(wt(:,i))
      enddo

      end subroutine

!----------------------------------------------------------------------------------------------------
      subroutine print_plotcoord(coord,x,y,wt,ymod,npt,id_plot,set)
      use dim_parameter,only: qn,ntot
      use io_parameters,only: llen
      implicit none
!     IN: variables
      integer, intent(in) :: id_plot,npt,set,coord
      double precision, intent(in) :: x(qn,npt),y(ntot,npt)
      double precision, intent(in) :: ymod(ntot,npt),wt(ntot,npt)
!     loop control
      integer i

      call print_plotheader(id_plot,coord,npt,set)
      do i=1,npt
!         write(id_plot,"(ES16.8,*(3(ES16.8),:))")
!     >        x(coord,i), ymod(:,i),y(:,i),(wt(:,i))
         write(id_plot,"(2ES16.8,*(3(ES16.8),:))")
     >        x(coord,i), x(coord+1,i),y(:,i)
      enddo

      end subroutine

!----------------------------------------------------------------------------------------------------
      subroutine print_plotheader(id_plot,coord,npt,set)
      use dim_parameter,only: qn,ntot
      use io_parameters,only: llen
      implicit none
      integer, intent(in) :: id_plot,npt,set,coord

      character(len=llen) fmt

      write(id_plot,'("#SET: ",i5)') set
      write(id_plot,'("#OUTPUT VALUES",i4)') ntot
      write(id_plot,'("#DATA POINTS: ",i4)') npt
      if (coord.le.0) then
      write(id_plot,'("#t(x) = WALK")')
      else
         write(id_plot,'("#t(x) = x(",I0,")")') coord
      endif
      write(id_plot,'("#UNIT: hartree")')
      write(id_plot,'()')
      write(id_plot,'("#",A15)',advance='no') "t(x)"
      write(fmt,'("(3(7X,A9,",I3,"(16x)))")') ntot-1
      write(id_plot,fmt) 'ymod(p,x)','y(x)  ','wt(x) '


      end subroutine

!----------------------------------------------------------------------------------------------------
!     <subroutine walknorm calulates the distance in coordinate space for each set
      subroutine getwalknorm(x,walknorm,npt)
      use dim_parameter,only: qn
      implicit none
!     IN: variables
      integer npt
      double precision x(qn,npt)
      double precision walknorm
!     INTERNAL: variables
      double precision xdiff(qn)
      integer i

      walknorm =0.d0
      do i=2,npt
         xdiff(1:qn) = x(1:qn,i) - x(1:qn,i-1)
         walknorm = walknorm + dsqrt(sum(xdiff(1:qn)**2))
      enddo

      end subroutine

!----------------------------------------------------------------------------------------------------
!     <Subroutine for generating output filenames and openeing the correspondign files
      subroutine open_plotfile(id_plot,set)
      implicit none
!     IN: Variables
      integer id_plot,set
!     INTERNAL: Variables
      character(len=30) name    !name of output file

!     define name sheme for plot files
      if (set .lt. 10 ) then 
            write(name,203) set 
      else
            write(name,202) set 
      endif
      
 202  format('scan',I2,'.dat')
 203  format('scan0',I1,'.dat')
      !write (name,202) set

c     open plotfile
      open(id_plot,file=name)

      end subroutine

!----------------------------------------------------------------------------------------------------
!     <Subroutine for generating output filenames and openeing the correspondign files
      subroutine open_outfile(id_out,it_makro)
      implicit none
      integer id_out,it_makro
      character(len=30) outname !name of output file

 543  format('mnlfit-',i1,'.out')
 544  format('mnlfit-',i2,'.out')
 545  format('mnlfit-',i3,'.out')

      if(it_makro.lt.10) then
         write(outname,543) it_makro
      else if (it_makro.lt.100) then
         write(outname,544) it_makro
      else if (it_makro.lt.1000) then
         write(outname,545) it_makro
      else
         write(6,*)
     >        'ERROR: No rule for Outputfile naming for MAXIT >= 1000'
         stop
      endif

      open (id_out,file=outname)

      end subroutine

!----------------------------------------------------------------------------------------------------
!     <Subroutine for printing the Parameterkeys for use in Input File
!     < prints the keystring given in keys.incl and the corresponding parameters when there was atleast one parameter given in the input for the spcific key
!     < how many parameters and spreads per line are printed can be specified with the hardcoded parameters np and nsp but they must be atleast >=2
!     <@param id_out specifies the file in which the Parameters are Printed
!     <@param p vector containing one set of parameter values
!     <@param p_act vector containing the active state 0 (inactive) or 1 (active) for each parameter
!     <@param p_spread vector containing the spreads for each parameter
!     <@param npar lenght of the parmeter vectors (p,p_act,p_spread)
!     <@TODO extract subroutine for printing the multiline values, would make this more readable
      subroutine print_parameterblock(id_out,p,p_act,p_spread,npar)
      use dim_parameter,only: pst, facspread
      use io_parameters,only: key, parkeynum,parkeylen,llen
      implicit none
!     IN: Variables
      integer  id_out,npar,p_act(npar)
      double precision  p(npar),p_spread(npar)

!     INTERNAL: variables
!     loop index
      integer i,k,l,t,n         !< internal variables for loops and positions in parameter vectors

!     number of values per line, values must be atleast 2 set this to personal preference
      integer, parameter :: np=5,nsp=5

      character(len=llen) fmt


!     Write header for Parameter block
 1    format('!',200('='))
      write(id_out,1)
      write(id_out,'(A2,5x,A11,i3)') '! ','PARAMETER: ',npar
      write(id_out,1)

!     loop over all Parameter Keys
      do i = 1, parkeynum
!        save start and end of parameter block for specific key
         k = pst(1,i)
         l = pst(1,i)+pst(2,i)-1
!        print only used keys with atleast one parameter
         if(pst(2,i).gt.0) then
            write(fmt,'("(a",I3,"'' ''i3)")') parkeylen
            write(id_out,fmt) adjustl(key(1,i)), pst(2,i)

!           write the actual parameters -> subroutine print_parameterlines()?
            if(l-k.le.(np-1)) then
               write(fmt,'("(a",I3,"'' ''",I3,"g24.15)")') parkeylen,np
               write(id_out,fmt) key(2,i),(p(n), n=k,l)

            else
!              start of multi line parameter print, number of values per line specified by np
               write(fmt,'("(a",I3,"'' ''",I3,"g24.15'' &'')")')
     $              parkeylen,np
               write(id_out,fmt) key(2,i),(p(n), n=k,k+(np-1))

               t=k+np
!              write continuation lines till left parameters fit on last line
               do while(t.le.l)
                  if(l-t.le.(np-1)) then
                     write(fmt,'("(",I3,"x'' ''",I3,"g24.15)")')
     $                    parkeylen,np
                     write(id_out,fmt) (p(n), n=t, l)

                  else
                     write(fmt,'("(",I3,"x'' ''",I3,"g24.15'' &'')")')
     $                    parkeylen,np
                     write(id_out,fmt) (p(n), n=t, t+(np-1))

                  endif
                  t=t+np
               enddo

            endif               !-> end subroutine print_parameterlines

!           write parameter active state in one line
            write(fmt,'("(a",I3,"'' ''","50i3)")') parkeylen
            write(id_out,fmt) key(3,i),(p_act(n),n=k,l)

!           write the spreads for each parameter
            if(l-k.le.(np-1)) then
               write(fmt,'("(a",I3,"'' ''",I3,"g24.8)")') parkeylen,nsp
               write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,l)

            else
!              start of multiline spread values
               write(fmt,'("(a",I3,"'' ''",I3,"g24.8'' &'')")')
     $              parkeylen,nsp
               write(id_out,fmt) key(4,i),(p_spread(n)/facspread, n=k,k
     >              +(np-1))

               t=k+nsp
!              write continuation lines till left spreads fit on last line
               do while(t.le.l)
                  if(l-t.le.(np-1)) then
                     write(fmt,'("(",I3,"x'' ''",I3,"g24.8)")')
     $                    parkeylen,nsp
                     write(id_out,fmt) (p_spread(n)/facspread, n=t, l)
                  else
                     write(fmt,'("(",I3,"x'' ''",I3,"g24.8'' &'')")')
     $                    parkeylen,nsp
                     write(id_out,fmt) (p_spread(n)/facspread, n=t, t
     >                    +(np-1))

                  endif
                  t=t+np
               enddo

            endif
!           print empty line between diffrent parameter blocks for better readability
            write(id_out,'(" ")')
         endif

      enddo

      end subroutine

!----------------------------------------------------------------------------------------------------
!     <Subroutine for printing the current Parameters and their active state
!     < prints only the numeric values of the parameters and does not specify the corresponding key
!     <@param npar number of parameter
!     <@param id_out specifies the output file
!     <@param p,p_act parameter vectors containing the values and the activity state of parameters
      subroutine print_parameterstate(id_out,p,p_act,npar)
      implicit none

!     IN: Variables
      integer npar,id_out
      double precision p(npar)
      integer p_act(npar)

!     INTERNAL: Variables
      integer i                 !< loop control
      integer nopt              !< number of counted active parameters
      character(len=16) opt(npar) !< string for optimisation state

!     initialize number of opt parameters and the string vector opt
      nopt=0
      opt = '   not opt.   '
!     loop over all parameters and check their active state count if active and set string to opt
      do i=1,npar
!     Nicole: change due to value 2 of p_act
!         if(p_act(i).eq.1) then
         if(p_act(i).ge.1) then
            opt(i) = '     opt.     '
            nopt=nopt+1
         endif
      enddo
!     print the Parameters and their active state within separating lines
      write(id_out,*)''
      write(id_out,block_line)
      write(id_out,*) 'Parameters:'
      write(id_out,sep_line)
      write(id_out,'(5g14.6)') (p(i),i=1,npar)
      write(id_out,'(5a14)')   (opt(i),i=1,npar)
      write(id_out,sep_line)
      write(id_out,'("No. of optimized parameters: ",i6)') nopt
      write(id_out,block_line)
      write(id_out,*)''
      end subroutine

!----------------------------------------------------------------------------------------------------
!     <Subroutine for printing coordinates,refdata,modeldata,diffrence between them and the weights
!     <@param id_out identiefies the output file
!     <@param x vector of input pattern for each datapoint
!     <@param y vector of expected output patterns for each datapoint
!     <@param ymod vector of output patterns generated by the model depending on paramerters
!     <@param wt vector of weights for each datapoint
!     <@param qn number of input patterns
!     <@param ntot total number of output patterns for each datapoint
!     <@param numdatpt number of totatl datapoints
!     <@param sets number of sets the datapoints are divided into
!     <@param ndata vector containing the number of included datapoints for each set
!     <@param i,j,point internal variables for loop controll and datapoint counting
      subroutine print_data(id_out,x,y,ymod,wt)
      use dim_parameter,only: sets,ndata,qn,ntot,numdatpt,qn_read
      implicit none
!     IN: Variables
      integer id_out
      double precision  x(qn,numdatpt)
      double precision  y(ntot,numdatpt),ymod(ntot,numdatpt)
      double precision  wt(ntot,numdatpt)

!     INTERNAL: Variables
      integer i,j,point

 18   format(A8,i6)
 19   format (3(A15,3x), 2x, A18  , 4x, A12)

!     print seperating line and header for Data output
      write(id_out,*) 'Printing Data Sets:'

      write(id_out,19) adjustl('y(x)'),adjustl('ymod(x)'),
     >     adjustl('y(x)-ymod(x)'),adjustl('weight'),
     >     adjustl('x(1:qn_read) ')
      write(id_out,sep_line)
!     loop over all datapoints for each set and count the actual datapointnumber with point
      point=0
      do i=1,sets
         write(id_out,18) 'Set: ', i
         do j=1,ndata(i)
            write(id_out,18) 'Point: ', j
            point=point+1
!           print all data for one datapoint
            call print_datapoint(id_out,x(:,point),y(:,point),
     >           ymod(:,point),wt(:,point))
            write(id_out,sep_line)

         enddo
      enddo
!     write end of data statement and two seperating lines
      write(id_out,block_line)
      write(id_out,*) ''
      end subroutine
!----------------------------------------------------------------------------------------------------
!     <Subroutine prints a single Datapoint splits Data in nstat nci(ndiab) blocks for readability
!     <@param id_out identiefies the output file
!     <@param x vector of input pattern for each datapoint
!     <@param y vector of expected output patterns for each datapoint
!     <@param ymod vector of output patterns generated by the model depending on paramerters
!     <@param wt vector of weights for each datapoint
!     <@param qn number of input patterns
!     <@param ntot total number of output patterns for each datapoint
!     <@param i,j,k internal variables for loop controll and counting
      subroutine print_datapoint(id_out,x,y,ymod,wt)
      use dim_parameter,only: nstat,ndiab,nci,qn,ntot,qn_read
      use io_parameters,only: llen
      implicit none
      integer id_out
      double precision  x(qn),y(ntot),ymod(ntot),wt(ntot)

      integer i,j,k

 18   format(A10,i3)
 19   format(3F18.8, 2X, F18.6, 4X,*(F12.6))

!     print the nstat output patterns
      do i=1,nstat
         write(id_out,19)y(i),ymod(i),ymod(i)-y(i), wt(i), x(1:qn)
      enddo
!     loop over number (nci) of  metadata with lenght (ndiab)
      do i=1,nci
         write(id_out,18) 'nci: ',i
         do j=1,ndiab
            k=nstat + (i-1)*ndiab + j
            write(id_out,19) y(k),ymod(k),(ymod(k)-y(k)),
     >           wt(k), x(1:qn_read)
         enddo
      enddo

      end subroutine

      

      end module write_mod