module ThermalBuilderAndSolverM

  use UtilitiesM
  use DebuggerM

  use SparseKit

  use NodePtrM
  use ElementPtrM
  use ConditionPtrM

  use LeftHandSideM

  use ThermalModelM

  use BuilderAndSolverM

  use LinearSolverM
  use mklPardisoM

  implicit none

  private
  public :: ThermalBuilderAndSolverDT, applyDirichlet, applyNewmann

  type, extends(NewBuilderAndSolverDT) :: ThermalBuilderAndSolverDT
   contains
     procedure :: buildAndSolve
  end type ThermalBuilderAndSolverDT

contains

  subroutine buildAndSolve(this, model)
    implicit none
    class(ThermalBuilderAndSolverDT), intent(inout) :: this
    class(ThermalModelDT)          , intent(inout) :: model
    integer(ikind) :: i
    write(*,*) '*** Direct Builder And Solver ***'
    call debugLog('  Assembling stiffness matrix and right hand side vector')
    print'(A)', 'Assembling stiffness matrix and right hand side vector'
    call assembleSystem(model)
    call applyBC(model)
  end subroutine buildAndSolve

  subroutine assembleSystem(model)
    implicit none
    class(ThermalModelDT)    , intent(inout) :: model
    integer(ikind)                           :: i, j, iElem, nElem, nNode, row, col
    type(LeftHandSideDT)                     :: localLHS
    real(rkind), dimension(:)  , allocatable :: localRHS
    real(rkind)                              :: adder
    type(ElementPtrDT)                       :: element
    nElem = model%getnElement()
    model%rhs = 0._rkind
    do iElem = 1, nElem
       element = model%getElement(iElem)
       nNode = element%getnNode()
       call element%calculateLocalSystem(model%processInfo, localLHS, localRHS)
       do i = 1, nNode
          row = element%getNodeID(i)
          adder = 0._rkind
          do j = 1, nNode
             col = element%getNodeID(j)
             call model%LHS%append(val = localLHS%stiffness(i,j)  &
                  , row = row                                    &
                  , col = col                                    )
             adder = adder + localLHS%mass(i,j)
          end do
          model%rhs(row) = model%rhs(row) + localRHS(i)
          call model%mass%append(val = adder                     &
               , row = row                                      &
               , col = row                                      )
       end do
       call localLHS%free()
       deallocate(localRHS)
    end do
    call model%lhs%makeCRS()
    call model%mass%makeCRS()
  end subroutine assembleSystem

  subroutine applyBC(model)
    implicit none
    class(ThermalModelDT), intent(inout) :: model
    call applyNewmann(model)
    call applyDirichlet(model)
  end subroutine applyBC

  subroutine applyNewmann(model)
    implicit none
    class(ThermalModelDT)      , intent(inout) :: model
    integer(ikind)                             :: i, j, iCond, nCond, nNode, row, col
    type(LeftHandSideDT)                       :: localLHS
    real(rkind), dimension(:)  , allocatable   :: localRHS
    type(ConditionPtrDT)                       :: condition
    nCond = model%getnCondition()
    do iCond = 1, nCond
       condition = model%getCondition(iCond)
       nNode = condition%getnNode()
       call condition%calculateLocalSystem(model%processInfo, localLHS, localRHS)
       if(condition%getAffectsLHS() == .true.) then
          do i = 1, nNode
             row = condition%getNodeID(i)
             do j = 1, nNode
                col = condition%getNodeID(j)
                call model%LHS%appendPostCRS(val = localLHS%stiffness(i,j)  &
                     , row = row                                           &
                     , col = col                                           )
             end do
          end do
       call localLHS%free()
       end if
       if(condition%getAffectsRHS() == .true.) then
          do i = 1, nNode
             row = condition%getNodeID(i)
             model%rhs(row) = model%rhs(row) + localRHS(i)
          end do
          deallocate(localRHS)
       end if
    end do
  end subroutine applyNewmann

  subroutine applyDirichlet(model)
    implicit none
    class(ThermalModelDT), intent(inout) :: model
    integer(ikind)                       :: i, nNode, nodeID
    type(NodePtrDT)                      :: node
    nNode = model%getnNode()
    do i = 1, nNode
       node = model%getNode(i)
       if(node%ptr%dof(1)%isFixed) then
          nodeID = node%ptr%getID()
          call model%lhs%setDirichlet(nodeID)
          model%rhs(nodeID) = node%ptr%dof(1)%fixedVal
       end if
    end do
  end subroutine applyDirichlet

end module ThermalBuilderAndSolverM


text:
      subroutine epimulti(n,epi, gamaEpi, sh,school,wght,nparam,
     $     par,parmin,parmax, step,ilog,imask,iseed,nsamps ,
     $     ithin,ndata,tps,rtn,curMin,coef,
     $     parCPL,pminCPL,pmaxCPL,stepCPL,ilogCPL,imaskCPL,fRij,tab,
     $     tabCPL,profiles,nRnd, ndays,epi_model, Temp)

      implicit none

      integer iday_per_week, nblock
      parameter(iday_per_week=7,nblock=3)

      integer nCPL
      parameter(nCPL=2)
      integer nstep
      parameter (nstep = 6)
! Y and gamay hold the synthetic profiles
      integer n, nparam, iseed
      integer ndays, ithin, ndata, nsamps
      integer epi_model, nRnd
      integer ilog(nparam),imask(nparam)
      real*8 epi(ndata,n), gamaEpi(ndata,n)
      real*8 sh(ndata,n), school(ndata,n)
      real*8 wght(ndata,n)
      real*8 par(nparam,n)
      real*8 parmin(nparam,n), parmax(nparam,n),step(nparam,n)
      real*8 Temp
!!    real*8 ymu(nparam,n), sigma(nparam,n)
      real*8 fN(n)
      real*8 pC(n)
      real*8 tps(ndata)
      real*8 e_bckgrnd(n)
      real*8 rtn(ndata,n)
      real  tab(nsamps/ithin,nparam+1,n)
      real  tabCPL(nsamps/ithin,nCPL)
      real  profiles(nRnd, ndata, n)
      real*8 rtnBest(ndata,n),rtnNew(ndata,n)
      real*8 dsdt(((ndays)*nstep),n)
      real*8 curpars(nparam,n),copypar(nparam,n),savepar(nparam,n)
      real*8 parupdt(nparam,n),parBest(nparam,n)
      real*8 curLLK(n), curMin(n)
      real*8 fnewLLK(n)
      real*8 ran1
      real*8 accept_rate(nblock)
      logical plog(nparam)
      integer i, j, k, m, icount
      integer noptv, ioptv(nparam)
      integer nopt(nblock),iopt(nparam,nblock)
      integer iaccept(nblock)
      integer iadapt(nblock)
      real*8 AUX
      real*8  range_min, range_max, scale, myaccept
      real*8 step_max, step_min
      real*8 tmp((nblock+1))
      integer ionep
      integer imid((ndata+1))
! Unique to the coupling between regions
      real*8 fRij(n,n), fMij(n,n), coef(n)
      real*8 parCPL(nCPL),pminCPL(nCPL)
      real*8 pmaxCPL(nCPL),stepCPL(nCPL)
      real*8 savefMij(n,n),parBestCPL(nCPL)
      real*8 saveparCPL(nCPL), copyparCPL(nCPL)
      real*8 curparsCPL(nCPL),parupdtCPL(nCPL)
      real*8 fnewLLK1D, curLLK1D,curMin1D
!      real*8 rtn1D(ndata),rtnNew1D(ndata)

      integer imaskCPL(nCPL), ilogCPL(nCPL)
      integer ioptCPL(nCPL), noptCPL
      integer iacceptCPL, iadaptCPL
      logical plogCPL(nCPL)

      real*8 tempCPL(sum(imaskCPL))
      integer iburn, iline, nlines
      external calcFit,ran1

      nlines = nsamps / ithin
! For an adaptive size MCMC - decide if we need to update step size every 1%
! or 1,000 steps
      scale = 2.0d0
      ionep = int(nsamps * 0.01)
      ionep = min(1000, ionep)
      range_min = 0.20d0
      range_max = 0.30d0
      step_max = 1.d0
      step_min = 1e-5

! ran1 works with -iseed

      iseed = -abs(iseed)

! ilog can have values of
! 0 - uniform
! 1 - log uniform
! 2 - Gaussian

! set an auxilary array with the indices of the parameters we are optimizing

      noptv = 0
      ioptv = 0
      nopt  = 0
      iopt = 0
      do i=1,nparam
         if (imask(i) > 0) Then
            noptv = noptv + 1
            ioptv(noptv) = i
         endif
      enddo

!
! Separate R0 and pC into their own block.  It helps..
!
      ilog(3) = 2
      ilog(5) = 2

      do i = 1,nparam
         if (imask(i) > 0) Then
            k = ilog(i) + 1
            nopt(k) = nopt(k) + 1
            iopt(nopt(k),k) = i
            do j = 1, n
            if (par(i,j) .le. parmin(i,j) .or. par(i,j)
     $              .ge. parmax(i,j)) Then
               par(i,j) = parmin(i,j) +
     $              0.50d0 * (parmax(i,j) - parmin(i,j))

            endif
         enddo
         endif
      enddo

      icount = 0

! convert 0 and 1 to true and false

      do i=1,nparam
         if (ilog(i) .eq. 0) Then
            plog(i) = .false.
         else
            plog(i) = .true.
         endif
      enddo

! For the coupling parameters
      noptCPL = 0
      do i =1,nCPL
         if(imaskCPL(i) > 0) Then
            noptCPL=noptCPL+1
            ioptCPL(noptCPL) = i
           if (parCPL(i) .le. pminCPL(i) .or. parCPL(i) .ge. pmaxCPL(i))
     $           Then
               parCPL(i) = pminCPL(i) +
     $         0.50d0 * (pmaxCPL(i) - pminCPL(i))
            endif
         endif
         if(ilogCPL(i) .eq. 0) Then
            plogCPL(i) = .false.
         else
            plogCPL(i) = .true.
         endif
      enddo

!
! Need to find the mid-point
! This will need to be further cleaned - the CDC data is absolute day numbers..
!

      if (tps(1) > (iday_per_week * 10)) Then
         call BuildIMIDforCDC(ndata, nstep, tps, imid)
      else
         call BuildIMID(ndata, nstep, tps, imid)
      endif

! the paramters are in the following order
! t0, pC, R0, seed,

      curpars = par

! integrate the ODEs to get the profiles

      fN = curpars(1, 1:n)
      pC = curpars(5, 1:n)
      e_bckgrnd = curpars(8, 1:n)

      fMij = 0.0d0
      tab = 0.0d0
      tabCPL = 0.0d0

      call buildMij(n,nCPL,fN,fRij,parCPL,fMij)

       select case (epi_model)
         case (1)
             call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $            curpars,fMij,dsdt)
         case (2)
             call RK4SEIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $            curpars,fMij,dsdt)
          case default
             call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $            curpars,fMij,dsdt)
      end select

      call weekly(n,ndata, ndays, nstep, imid, dsdt, pC,
     $     e_bckgrnd, rtn)

! calculate Likelihood of the solution

      call calcFit(n,epi,gamaEpi,rtn,wght,ndata,curLLK)

      curLLK1D = dot_product(curLLK, coef) /dble(n *Temp)

      curMIN = curLLK
      curMIN1D = curLLK1D
      rtnBest = rtn

      iaccept = 0
      iadapt = 0
      icount = 0

      iacceptCPL = 0
      iadaptCPL = 0

      parupdt = curpars
      parupdtCPL = parCPL
      curparsCPL = parCPL

      tab = 0.0d0
      tabCPL = 0.0d0

! MCMC loop starts here
      do i =1, nsamps

         do k = 1, nblock

            if (nopt(k) .eq. 0) go to 102
            savepar = curpars
            copypar = curpars

            call fnProposeParamUpdatesBlock(n,nparam,copypar,
     $           parmin,parmax,step,plog,parupdt,
     $           iseed,nopt(k),iopt(1:nopt(k),k))

!     update all the parameters and calculate a new solution and LLK

            curpars(iopt(1:nopt(k),k),:)=parupdt(iopt(1:nopt(k),k),:)

            fN = curpars(1, 1:n)
            pC = curpars(5, 1:n)
            e_bckgrnd = curpars(8, 1:n)

            select case (epi_model)
            case (1)
              call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $            curpars,fMij,dsdt)
            case (2)
               call RK4SEIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $            curpars,fMij,dsdt)
            case default
               call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $              curpars,fMij,dsdt)
            end select

            rtnNew = 0.0d0
            call weekly(n, ndata, ndays, nstep, imid, dsdt, pC,
     $           e_bckgrnd,rtnNew)

            call calcFit(n, epi,gamaEpi,rtnNew,wght,ndata,
     $           fnewLLK)

            fnewLLK1D = dot_product(fnewLLK, coef)/dble(n * Temp)

            AUX = 1.0d0
c$$$            if (k .eq. 3) Then
c$$$               AUX = Ratios(nparam,iopt(1:nopt(k),k),
c$$$     $              curpars,savepar,ymu,sigma,Temp(1))
c$$$            endif

            call MH(n,fnewLLK1D,curLLK1D,AUX,curMin1D,nparam,ndata,
     $           iseed,iaccept(k),iadapt(k),curpars,
     $           savepar,parBest,rtnNew,rtnBest,fnewLLK,curLLK,curMin)

 102        continue

         enddo ! End of loop over blocks

!! Now updated the coupling parameters
         if (noptCPL .eq. 0) go to 103
         savefMij = fMij
         saveparCPL = curparsCPL
         copyparCPL = curparsCPL

         call fnProposeParamUpdates(nCPL,copyparCPL,
     $        pminCPL,pmaxCPL,stepCPL,plogCPL,
     $        parupdtCPL,iseed,noptCPL,ioptCPL(1:noptCPL))

         curparsCPL(ioptCPL(1:noptCPL))=
     $        parupdtCPL(ioptCPL(1:noptCPL))

         fN = curpars(1, 1:n)
         pC = curpars(5, 1:n)
         e_bckgrnd = curpars(8, 1:n)

         call buildMij(n,nCPL,fN,fRij,curparsCPL,fMij)
         select case (epi_model)
         case (1)
             call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $            curpars,fMij,dsdt)
         case (2)
            call RK4SEIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $           curpars,fMij,dsdt)
         case default
            call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $           curpars,fMij,dsdt)
         end select

         rtnNew = 0.0d0
         call weekly(n, ndata, ndays, nstep, imid, dsdt, pC,e_bckgrnd,
     $        rtnNew)

         call calcFit(n, epi,gamaEpi,rtnNew,wght,ndata,
     $        fnewLLK)

         fnewLLK1D = dot_product(fnewLLK , coef) /dble(n * Temp)

         call MH_CPL(n,fnewLLK1D,curLLK1D,AUX,curMin1D,nCPL,ndata,
     $        iseed,iacceptCPL,iadaptCPL,curparsCPL,
     $        saveparCPL,parBestCPL,fMij,savefMij,rtnNew,rtnBest,
     $        fnewLLK,curLLK,curMin)

 103     continue

!! End of adjusting coupling parameters

         if (mod(i,ithin) .eq. 0) Then
            icount = icount + 1
            tab(icount,1:nparam,:)   = real(curpars)
            tab(icount,(nparam+1),:) = real(curLLK)
            tabCPL(icount,:) = real(curparsCPL)
         endif

         if (mod(i,ionep) .eq. 0) Then
            m = 0
            do k = 1,nblock
               if (nopt(k) .eq. 0) go to 202
               myaccept = dble(iadapt(k))/dble(ionep)

               if (myaccept > range_max .and.
     $         all(step(iopt(1:nopt(k),k),:)*scale < step_max)) Then
                 step(iopt(1:nopt(k),k),:) =
     $                 step(iopt(1:nopt(k),k),:) *scale
               endif
               if (myaccept < range_min .and.
     $         all(step(iopt(1:nopt(k),k),:)/scale > step_min)) Then
                 step(iopt(1:nopt(k),k),:) =
     $                 step(iopt(1:nopt(k),k),:) /scale
               endif
               iadapt(k) = 0    !reset acceptance number
               m = m + 1
               tmp(m) = dble(iaccept(k))/dble(i)*100.0
 202           continue

            enddo

!! Repeat for coupling
            if (noptCPL. eq. 0) go to 302
            myaccept = dble(iadaptCPL)/dble(ionep)

            if (myaccept > range_max) Then
              tempCPL = stepCPL(ioptCPL)
              where (tempCPL*scale < step_max)
               tempCPL = tempCPL*scale
              end where
              stepCPL(ioptCPL) = tempCPL
            endif
            if (myaccept < range_min ) Then
              tempCPL = stepCPL(ioptCPL)
              where (tempCPL/scale > step_min)
                tempCPL = tempCPL/scale
              end where
              stepCPL(ioptCPL) = tempCPL
            endif

            iadaptCPL = 0       !reset acceptance number

            m = m + 1
            tmp(m) = dble(iacceptCPL)/dble(i)*100.0

 302        continue
            call intpr('MCMC Step Number:',-1,i,1)
            call dblepr('curLLK1D',-1,curLLK1D,1)
            call dblepr('accept%',-1,tmp,m)
         endif

      enddo

!update change information only every ithin iterations

      accept_rate = dble(iaccept)/dble(nsamps) *100.0d0

      call dblepr('acceptence %',-1,accept_rate,nblock)

! generate the profiles _ could really be done on the fly in the MCMC procedure
! see aove

      iburn = nlines/5
      do i = 1, nRnd
         iline = iburn + floor(ran1(iseed)*(nlines - iburn))
         curpars = tab(iline, 1:nparam, 1:n)
         parCPL  = tabCPL(iline,:)

         fN = curpars(1, 1:n)
         pC = curpars(5, 1:n)
         e_bckgrnd = curpars(8, 1:n)

         fMij = 0.0d0

         call buildMij(n,nCPL,fN,fRij,parCPL,fMij)

         select case (epi_model)
         case (1)
            call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $           curpars,fMij,dsdt)
         case (2)
            call RK4SEIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $           curpars,fMij,dsdt)
         case default
            call RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $           curpars,fMij,dsdt)
         end select

         call weekly(n,ndata, ndays, nstep, imid, dsdt, pC,
     $        e_bckgrnd, rtn)

         profiles(i, 1:ndata, 1:n) = real(rtn)
      enddo

      rtn = rtnBest
      par = parBest

      return

      end subroutine epimulti

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

      subroutine RK4SIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $     curpars,fMij,dsdt)

      implicit none
      integer iday_per_week
      parameter(iday_per_week=7)
      integer n, ndata, ndays, nstep, nparam
      real*8 sh(ndata,n), school(ndata,n)
      real*8 curpars(nparam,n)
      real*8 deltaR(n), Aparam(n), alpha(n)
      real*8 delta(n), ts(n), dur(n)
      real*8 dt
      real*8 t0(n), R0(n), Tg, sigma(n)
      real*8 fN(n), fMij(n,n)
      real*8 dsdt((ndays)*nstep,n)
      real*8 seed(n) !initial number of infectious
      real*8 tps(ndata) !time-series
      real*8 tps2(0:(ndata+1))
      real*8 y(4,n), tmpY(4,n)
      real*8 dy1(4,n), dy2(4,n), dy3(4,n), dy4(4,n)
      real*8 lambda(n)
      real*8 t_cur, Rt(n)
      real*8 p5, p3, p6
      integer np(n), index
      integer k, iweek, istep, iday
      integer icount, istart, iend

      dt = 1.0d0 / dble(nstep)

      p5 = 1.0d0 / 2.0d0
      p3 = 1.0d0 / 3.0d0
      p6 = 1.0d0 / 6.0d0

      fN = curpars(1,:)
      Tg  = curpars(2,1)
      R0 = curpars(3,:)
      sigma = curpars(4,:)
!pC = curpars(5,:)
      t0 = curpars(6,:)
      seed = curpars(7,:)
!e_bckgrnd = curpars(8,:)
      deltaR = curpars(9,:)
      Aparam = curpars(10,:)
      alpha  = curpars(11,:)
      delta = curpars(12,:)
      ts = curpars(13,:)
      dur = curpars(14,:)

! initialize things - this is the Y vector

      y(1,1:n) = fN - seed !fN * fraction-1.0
      y(2,1:n) = seed
      y(3,1:n) = fN   !* (1.0d0 - fraction)
      y(4,1:n) = 0.0d0

      dsdt = 0.0d0

      tps2 = 0.0d0
      tps2(1:ndata) = tps
      tps2(0) = tps2(1) - (tps(2) - tps(1))
      tps2((ndata+1)) = tps2(ndata) +  (tps(ndata) - tps((ndata-1)))

c$$$      t_cur = tps(1)
      t_cur = tps2(0)

      icount = 0

!     Loop over ndays and in each day loop over nsteps
!     need to asign each day to a week because the sh and school data are weekly
!     tps has an arbitrary start (could be week 27 for example) but sh/school
!     go from week 1 to week ndata hence we substract the start

       istart = int(tps2(0))
       iend   = int(tps2((ndata+1)))

       do iday = istart, iend
          iweek = (iday-int(tps2(0)))/iday_per_week + 1
          iweek = min(iweek,ndata)
          iweek = max(1,iweek)

          do istep=1,nstep
             t_cur = t_cur+dt
             index = 0
             dY1 = 0.0d0
             dY2 = 0.0d0
             dY3 = 0.0d0
             dY4 = 0.0d0
             tmpY = 0.0d0
! calculate the time dependent R(t)
             Rt = R0
     $            * (1.0d0+deltaR*exp(-Aparam*sh(iweek,1:n)))
     $            * (1.0d0-alpha*school(iweek,1:n))

             do k = 1,n
                if (t_cur > t0(k)) then
                   index = index+1
                   np(index) = k
                end if

                if (t_cur >= ts(k) .and. t_cur <= (ts(k)+dur(k))) Then
                   Rt(k) = Rt(k) * (1.0d0 + delta(k))
                endif
             enddo

             Rt = Rt / Tg

             call calcForce(n,y(2,:),fN,fMij,Rt,lambda)

             call diffSIR(n, np,index, lambda, Tg, Y, dY1)

             tmpY = Y + dt * dY1 * p5

             call calcForce(n,tmpY(2,:),fN,fMij,Rt,lambda)

             call diffSIR(n, np,index, lambda, Tg, tmpY, dY2)

             tmpY = Y + dt * dY2 * p5

             call calcForce(n,tmpY(2,:),fN,fMij,Rt,lambda)

             call diffSIR(n, np,index, lambda, Tg, tmpY, dY3)

             tmpY = Y + dt * dY3

             call calcForce(n,tmpY(2,:),fN,fMij,Rt,lambda)

             call diffSIR(n, np,index, lambda, Tg, tmpY, dY4)

             tmpY = Y + dt * ( dY1 * p6 + dY2 * p3 +
     $              dY3 * p3 + dY4 * p6 )
             Y = tmpY

!     Update -dS/dt
             icount= icount + 1

             dsdt(icount,1:n) = y(4,1:n)
          enddo                 ! End of a single day
        enddo                    ! End of all days

       return
      end subroutine RK4SIR

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

      subroutine RK4SEIR(n,sh,school,ndata,ndays,nstep,tps,nparam,
     $     curpars,fMij,dsdt)

      implicit none
      integer iday_per_week
      parameter(iday_per_week=7)
      integer n, ndata, ndays, nstep, nparam
      real*8 sh(ndata,n), school(ndata,n)
      real*8 curpars(nparam,n)
      real*8 deltaR(n), Aparam(n), alpha(n)
      real*8 delta(n), ts(n), dur(n)
      real*8 dt
      real*8 t0(n), R0(n), Tg, sigma(n)
      real*8 fN(n), fMij(n,n)
      real*8 dsdt((ndays)*nstep,n)
      real*8 seed(n) !initial number of infectious
      real*8 tps(ndata) !time-series
      real*8 tps2(0:(ndata+1))
      real*8 y(5,n), tmpY(5,n)
      real*8 dy1(5,n), dy2(5,n), dy3(5,n), dy4(5,n)
      real*8 lambda(n)
      real*8 t_cur, Rt(n)
      real*8 p5, p3, p6
      integer np(n), index
      integer k, iweek, istep, iday
      integer icount, istart, iend

      dt = 1.0d0 / dble(nstep)

      p5 = 1.0d0 / 2.0d0
      p3 = 1.0d0 / 3.0d0
      p6 = 1.0d0 / 6.0d0

      fN = curpars(1,:)
      Tg  = curpars(2,1)
      R0 = curpars(3,:)
      sigma = curpars(4,:)
!pC = curpars(5,:)
      t0 = curpars(6,:)
      seed = curpars(7,:)
!e_bckgrnd = curpars(8,:)
      deltaR = curpars(9,:)
      Aparam = curpars(10,:)
      alpha  = curpars(11,:)
      delta = curpars(12,:)
      ts = curpars(13,:)
      dur = curpars(14,:)

! initialize things - this is the Y vector S-E-I-R-cummulative

      y(1,1:n) = fN - seed !fN * fraction-1.0
      y(2,1:n) = seed
      y(3,1:n) = 0.0d0   !* (1.0d0 - fraction)
      y(4,1:n) = 0.0d0
      y(5,1:n) = 0.0d0
      dsdt = 0.0d0

      tps2 = 0.0d0
      tps2(1:ndata) = tps
      tps2(0) = tps2(1) - (tps(2) - tps(1))
      tps2((ndata+1)) = tps2(ndata) +  (tps(ndata) - tps((ndata-1)))

c$$$      t_cur = tps(1)
      t_cur = tps2(0)

      icount = 0

!     Loop over ndays and in each day loop over nsteps
!     need to asign each day to a week because the sh and school data are weekly
!     tps has an arbitrary start (could be week 27 for example) but sh/school
!     go from week 1 to week ndata hence we substract the start

       istart = int(tps2(0))
       iend   = int(tps2((ndata+1)))

       do iday = istart, iend
          iweek = (iday-int(tps2(0)))/iday_per_week + 1
          iweek = min(iweek,ndata)
          iweek = max(1,iweek)

          do istep=1,nstep
             t_cur = t_cur+dt
             index = 0
             dY1 = 0.0d0
             dY2 = 0.0d0
             dY3 = 0.0d0
             dY4 = 0.0d0
             tmpY = 0.0d0
! calculate the time dependent R(t)
             Rt = R0
     $            * (1.0d0+deltaR*exp(-Aparam*sh(iweek,1:n)))
     $            * (1.0d0-alpha*school(iweek,1:n))

             do k = 1,n
                if (t_cur > t0(k)) then
                   index = index+1
                   np(index) = k
                end if

                if (t_cur >= ts(k) .and. t_cur <= (ts(k)+dur(k))) Then
                   Rt(k) = Rt(k) * (1.0d0 + delta(k))
                endif
             enddo

             Rt = Rt / Tg

             call calcForce(n,y(3,:),fN,fMij,Rt,lambda)

             call diffSEIR(n,np,index, lambda, sigma, Tg, Y, dY1)

             tmpY = Y + dt * dY1 * p5

             call calcForce(n,tmpY(3,:),fN,fMij,Rt,lambda)

             call diffSEIR(n,np,index, lambda, sigma, Tg, tmpY, dY2)

             tmpY = Y + dt * dY2 * p5

             call calcForce(n,tmpY(3,:),fN,fMij,Rt,lambda)

             call diffSEIR(n,np,index, lambda, sigma, Tg, tmpY, dY3)

             tmpY = Y + dt * dY3

             call calcForce(n,tmpY(3,:),fN,fMij,Rt,lambda)

             call diffSEIR(n,np,index, lambda, sigma, Tg, tmpY, dY4)

             tmpY = Y + dt * ( dY1 * p6 + dY2 * p3 +
     $              dY3 * p3 + dY4 * p6 )
             Y = tmpY

!     Update -dS/dt
             icount= icount + 1

             dsdt(icount,1:n) = y(5,1:n)
          enddo                 ! End of a single day
        enddo                    ! End of all days

       return
      end subroutine RK4SEIR

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

!-------------------------------------------------------
C Calcluate the force of infection for each compartment

      subroutine calcForce(n,fI,fN,fMij,beta,lambda)

      implicit none

      integer n
      real*8 fI(n), fN(n), fMij(n,n)
      real*8 beta(n), lambda(n)
      real*8 sum_numer,sum_denom, ratio(n)
      integer i, j, k

      lambda = 0.0d0

      do j=1,n
         sum_numer = 0.0d0
         sum_denom = 0.0d0

         do k=1,n

            sum_numer = sum_numer + fMij(k,j) * fI(k)
            sum_denom = sum_denom + fMij(k,j) * fN(k)
         enddo
         ratio(j) = sum_numer/sum_denom
      enddo

      do i = 1,n
         do j = 1,n
            lambda(i) = lambda(i) + beta(j) * fMij(i,j) *
     $                 ratio(j)
         enddo

      enddo

      return
      end subroutine calcForce

!------------------------------------
C The Main Differential Equations- SIR

      SUBROUTINE DiffSIR(n, np,index, lambda, Tg, Pop, dPop)

      implicit none

      integer n, np(n), index
      REAL*8 lambda(n), Tg
      REAL*8 Pop(4,n), dPop(4,n)
      integer i, j
C
C     The differential equations
C

C     dS/dt =
      do i = 1,index
         j = np(i)
         dPop(1,j) = - lambda(j)*Pop(1,j)
C     dI/dt =
         dPop(2,j) =   lambda(j)*Pop(1,j)
     $        - Pop(2,j)/Tg
C     dR/dt =
         dPop(3,j) =   Pop(2,j)/Tg
c cummulative dsdt
         dPop(4,j) = lambda(j)*Pop(1,j)
      enddo
      RETURN
      END subroutine DiffSIR

!------------------------------------
C The Main Differential Equations- SIR

      SUBROUTINE DiffSEIR(n,np,index, lambda, sigma, Tg, Pop, dPop)

      implicit none

      integer n, np(n), index
      REAL*8 lambda(n), sigma(n), Tg
      REAL*8 Pop(5,n), dPop(5,n)
      integer i, j
C
C     The differential equations
C

C     dS/dt =
      do i = 1,index
         j = np(i)
         dPop(1,j) = -lambda(j)*Pop(1,j)
C     dE/dt =
         dPop(2,J) = lambda(j)*Pop(1,j)
     $        - sigma(j) * Pop(2,j)
C     dI/dt =
         dPop(3,j) = sigma(j)*Pop(2,j)
     $        - Pop(3,j)/Tg
C     dR/dt =
         dPop(4,j) = Pop(3,j)/Tg
c cummulative dsdt
         dPop(5,j) = sigma(j) * Pop(2,j)
      enddo
      RETURN
      END subroutine DiffSEIR

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

        subroutine weekly(n,ndata,ndays,nstep,imid,dsdt,pC,e_bckgrnd,x)

        implicit none

        integer n,ndata, ndays, nstep
        real*8 dsdt((ndays)*nstep,n)
        real*8 x(ndata,n), pC(n), e_bckgrnd(n)
        integer imid(ndata+1)
        integer i, istart, iend
        integer iday_per_week, ishift, j

c$$$        do i=1,ndata
c$$$           iend   = imid(i+1)
c$$$           istart = imid(i)
c$$$
c$$$           x(i,1:n) = dsdt(iend,1:n) - dsdt(istart,1:n)
c$$$
c$$$           x(i,1:n) = x(i,1:n)* pC(1:n) +e_bckgrnd(1:n)
c$$$        enddo

        iday_per_week = 7

        ishift = nstep * iday_per_week * 0.50d0
        do i=1,ndata
           do j = 1,n
              x(i,j) = dsdt(i*nstep*iday_per_week+ishift,j) -
     $             dsdt(i*nstep*iday_per_week-ishift,j)
           enddo
           x(i,1:n) = x(i,1:n)* pC(1:n) +e_bckgrnd(1:n)
        enddo

        return
        end subroutine weekly

C--------------------------------------------------------------------------------
       subroutine calcFit(n,y,gamay,x,wght,ndata,curLLK)

        implicit none

        integer n, ndata, i, j

        real*8 y(ndata,n),x(ndata,n), gamay(ndata,n)
        real*8 wght(ndata,n)
        real*8 xi, yi,val,CurLLK(n)

c x is the simulated data
c y is the base profile

C calculate the P(yi,xi)

           do j = 1,n
              curLLK(j) = 0.0d0
              do i=1,ndata

                 yi = y(i,j)
                 xi = x(i,j)

                 val = yi * log(xi) - xi - gamay(i,j)
                 curLLK(j) = curLLK(j) + val * wght(i,j)

              enddo
              curLLK(j) = -curLLK(j) /sum(wght(:,j))
           enddo

        return

        end subroutine CalcFit

!
! -------------------------------------------------------------------
!
      subroutine MH(n,fnewLLK1D,curLLK1D,AUX,curMin1D,nparam,
     $     ndata,iseed,iaccept,iadapt,curpars,
     $     savepar,parBest,rtn,rtnBest,fnewLLK,curLLK,curMin)

      implicit none
      integer n,nparam,ndata,iseed,iaccept,iadapt
      real*8 fnewLLK1D, curLLK1D, AUX,curMin1D
      real*8 curMin(n), curLLK(n), fnewLLK(n)
      real*8 curpars(nparam,n),savepar(nparam,n)
      real*8 parBest(nparam,n)
      real*8 rtn(ndata,n), rtnBest(ndata,n)
      real*8 rnd,ran1,diff_LLK
      external ran1
      logical accept

      diff_LLK = fnewLLK1D - curLLK1D

      accept = .false.
      if (diff_LLK .le. 0.0d0 .and. AUX .ge. 1.0d0) Then
         accept = .true.
      else
         rnd = ran1(iseed)
         if ((exp(-diff_LLK)*AUX) .gt.rnd) Then
            accept = .true.
         else
            accept = .false.
         endif
      endif

! if step is accepted

      if (accept) Then
         iaccept = iaccept + 1
         iadapt= iadapt + 1
         savepar = curpars
         curLLK1D = fnewLLK1D
         curLLK   = fnewLLK

         if (curLLK1D .le. curMin1D) then !keep track of the best profile we have
            parBest = curpars
            curMin1D = curLLK1D
            curMin   = curLLK
            rtnBest = rtn
         endif
!     if step is rejected - restore saved values of parameters
      else

         curpars = savepar

      endif

      return
      end subroutine MH

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

      subroutine MH_CPL(n,fnewLLK1D,curLLK1D,AUX,curMin1D,nparam,
     $     ndata,iseed,iaccept,iadapt,curpars,savepar,
     $     parBest,curfMij,savefMij,rtn,rtnBest,fnewLLK,curLLK,curMin)

      implicit none
      integer n,nparam,ndata,iseed,iaccept,iadapt
      real*8 fnewLLK1D, curLLK1D, AUX,curMin1D
      real*8 fnewLLK(n), curLLK(n), curMin(n)
      real*8 curpars(nparam),savepar(nparam)
      real*8 parBest(nparam)
      real*8 curfMij(n,n), savefMij(n,n)
      real*8 rtn(ndata,n), rtnBest(ndata,n)
      real*8 rnd,ran1,diff_LLK
      external ran1
      logical accept

      diff_LLK = fnewLLK1D - curLLK1D

      accept = .false.
      if (diff_LLK .le. 0.0d0 .and. AUX .ge. 1.0d0) Then
         accept = .true.
      else
         rnd = ran1(iseed)
         if ((exp(-diff_LLK)*AUX) .gt.rnd) Then
            accept = .true.
         else
            accept = .false.
         endif
      endif

! if step is accepted

      if (accept) Then
         iaccept = iaccept + 1
         iadapt = iadapt + 1
         savepar = curpars
         curLLK1D = fnewLLK1D
         curLLK   = fnewLLK
         savefMij = curfMij

         if (curLLK1D .le. curMin1D) then !keep track of the best profile we have
            parBest= curpars
            curMin1D = curLLK1D
            curMin = fnewLLK
            rtnBest = rtn
         endif
!     if step is rejected - restore saved values of parameters
      else

         curpars = savepar
         curfMij = savefMij

      endif

      return
      end subroutine MH_CPL

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

      subroutine SR_from_unit_Block_Mult(n,nopt,x,ymin,ymax,
     $     logflag,rtn)

      implicit none

      integer n, nopt
      real*8 x(nopt,n),ymin(nopt,n),ymax(nopt,n),rtn(nopt,n)
      integer i, j
      logical logflag

! Block version
! Only allow for log base 10

      if (logflag) Then
            do i = 1,nopt
               do j = 1,n
                  rtn(i,j) = ymin(i,j) *
     $           10.0**(x(i,j)*(log10(ymax(i,j))-log10(ymin(i,j))))
               enddo
            enddo
      else
         do i=1,nopt
            do j = 1,n
               rtn(i,j) = ymin(i,j) + (ymax(i,j)-ymin(i,j))*x(i,j)
            enddo
         enddo
      endif

      return
      end subroutine SR_from_unit_Block_Mult

c----------------------------------------------------------------

      subroutine SR_to_unit_Block_Mult(n,nopt,y,ymin,ymax,
     $     logflag,rtn)

      implicit none
      integer n,nopt
      real*8 y(nopt,n), ymin(nopt,n),ymax(nopt,n),rtn(nopt,n)
      integer i,j
      logical logflag

! Block version
! Only allow for log base 10

      if (logflag) Then

         do i = 1,nopt
            do j = 1,n
               rtn(i,j) = (log10(y(i,j)) - log10(ymin(i,j))) /
     $              (log10(ymax(i,j))-log10(ymin(i,j)))
            enddo
         enddo
      else
         do i = 1, nopt
            do j = 1, n
               rtn(i,j) = (y(i,j) - ymin(i,j))/(ymax(i,j) - ymin(i,j))
            enddo
         enddo
      endif

      return
      end subroutine SR_to_unit_Block_Mult

c----------------------------------------------------------------

      subroutine fnProposeParamUpdatesBlock(n,nparam,curval,
     $     valmin,valmax,step,logflag,
     $     parupdt,iseed,nopt,iopt)

      implicit none
      integer n,nparam, nopt, iopt(nopt)
      real*8 curval(nparam,n),valmin(nparam,n),valmax(nparam,n)
      real*8 parupdt(nparam,n),step(nparam,n)
      real*8 x(nopt,n), rv(nopt,n), rtn(nopt,n)
      real*8 ran1
      integer iseed
      integer j,k

      logical logflag(nparam)
      external ran1

      parupdt = curval

      do j = 1, nopt
         do k = 1,n
            rv(j,k) = ran1(iseed)
            rv(j,k) = (rv(j,k) - 0.50d0) * step(iopt(j),k)
         enddo
      enddo

      call SR_to_unit_Block_Mult(n,nopt,curval(iopt,:),valmin(iopt,:),
     $         valmax(iopt,:),logflag(1),x)

      x = x + rv

      do j = 1, nopt
         do k = 1,n
            if (x(j,k) .le. 0.0d0) x(j,k) = 1.0d0 + x(j,k)
            if (x(j,k) .ge. 1.0d0) x(j,k) = x(j,k) - 1.0d0
         enddo
      enddo

      call SR_from_unit_Block_Mult(n,nopt,x,valmin(iopt,:),
     $     valmax(iopt,:),logflag(1),rtn)

      parupdt(iopt,:) = rtn

      return
      end subroutine fnProposeParamUpdatesBlock

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

      subroutine buildMij(n,nCPL,fN,fRij,parCPL,fMij)

      implicit none

      integer n,nCPL
      real*8 parCPL(nCPL)
      real*8 fMij(n,n),fN(n),fRij(n,n)
      real*8 tmp(n,n),sd, gamma, denom
      integer i, j

      sd    = parCPL(1)
      gamma = parCPL(2)

      tmp = 0.0d0
      tmp = fRij/sd
! The Mills & Riley Kernel
      do i=1,n
         do j=i,n
            denom = 1.0d0 + tmp(i,j)**gamma
            fMij(i,j) = fN(j)/denom
            fMij(j,i) = fN(i)/denom
         enddo
         fMij(i,:) = fMij(i,:) /sum(fMIJ(i,:))
      enddo

c$$$      fMij = 0.0d0
c$$$      do i = 1, n
c$$$         fMij(i,i) = 1.0d0
c$$$      enddo

      return
      end subroutine buildMij

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

      subroutine checkParam(n,nparam,pmin,pmax,pval)

      implicit none

      integer n, nparam
      real*8 pmin(nparam,n), pmax(nparam,n)
      real*8 pval(nparam,n)
      integer i, j

      do i=1,nparam
         do j = 1,n
            if (pval(i,j) .ge. pmax(i,j)) Then
               call dblepr('Greater than Pmax',-1,pval(i,j),1)
               call intpr('i',-1,i,1)
               call intpr('j',-1,j,1)
               stop
            endif

            if (pval(i,j) .le. pmin(i,j)) Then
               call dblepr('Smaller than Pmin',-1,pval(i,j),1)
               call intpr('i',-1,i,1)
               call intpr('j',-1,j,1)
               stop
            endif

         enddo
      enddo

      return
      end subroutine checkParam

ext: f


text:
      MODULE MOD_PRICORR

      CONTAINS

      SUBROUTINE PRICORR(POPNUM, POP1, LEVEL, N, ND, MODHEAD, LSPOP, CORMAX_FIN,
     $                   NCHARG, RNE, JOBNUM, REDUCE, GAMMAL, GAMMAR, T,
     $                   DELTAC, ELEMENT, NATOM, NFIRST, NLAST)

!     PRINTOUT OF CORRECTION FACTORS OF POPULATION NUMBERS
!     RELATIVE TO THE LAST ITERATION

      USE ABUNDANCES

      use mod_tictoc

      USE FILE_OPERATIONS
      USE COMMON_BLOCK

      IMPLICIT REAL*8(A - H, O - Z)

      INTEGER, INTENT(IN) :: JOBNUM

      DIMENSION POPNUM(ND, N), POP1(ND, N)
      DIMENSION NCHARG(N), RNE(ND), T(ND)
      DIMENSION NFIRST(NATOM), NLAST(NATOM)

      REAL*8, ALLOCATABLE :: ABXYZ_new(:)

      CHARACTER LEVEL(N)*10
      CHARACTER ELEMENT(NATOM)*10
      CHARACTER MODHEAD*100
      CHARACTER PRILINE*130, NUMBER*12

      REAL*8,  DIMENSION(NATOM) :: CORMAX_ELE
      REAL*8,  DIMENSION(N) ::     CORMAX_LEV

      INTEGER, DIMENSION(NATOM) :: DCM_ELE

      REAL*8,  INTENT(OUT) ::      CORMAX_FIN

      REAL*8 ::                    CORMAX_ALL
      REAL*8 ::                    CORMAX_NHE
      REAL*8 ::                    CORMAX_NHN

      INTEGER ::                   DCM_ALL
      INTEGER ::                   DCM_NHE
      INTEGER ::                   DCM_NHN

      REAL*8 ::                    CM

      INTEGER ::                   DCM, DCM_FIN

      CHARACTER(:), ALLOCATABLE :: CONV_FILE

      CHARACTER(LEN = 10) ::       ELE

      character(len = 8)  ::       time

      INTEGER, ALLOCATABLE, DIMENSION(:) :: DEPTHS

      INTEGER :: DEPTHS_DIM, DEPTH_IND

      LOGICAL :: ALLE, NOHE, NOHN

      LOGICAL NEGPOPL
C***  COMMON /GIIIERR/ COUNTS THE ERROR MESSAGES FROM SUBR. GAUNTFF
      COMMON /GIIIERR/  NTUP,NTLOW,NFUP,NFLOW,NHELP
C***  NEGINTL COUNTS THE ERROR MESSAGES  FROM SUBR. LINPOP
      COMMON / COMNEGI / NEGINTL

!      START = LAMBDA_ITER .EQ. 1! .AND. .NOT. OLDSTART .AND. .NOT. LBKG

      CORMAX_ELE(:) = 0.0D0
      CORMAX_LEV(:) = 0.0D0

      DCM_ELE(:) =    0

      IF (LSPOP .GE. 1) PRINT 30
   30 FORMAT (10H1$$$$$$$$$)
      PRINT 1,MODHEAD,JOBNUM
    1 FORMAT (              /,1X,  A  ,20X,'JOB NO.',I5,//,10X,
     $ 'RATIO OF POPULATION NUMBERS TO THOSE OF THE LAST ITERATION:')
      IF (LSPOP.LT.1) GOTO 8
      J1=1
    4 J2=MIN0(N,J1+9)
      PRINT 2,(LEVEL(J),J=J1,J2)
    2 FORMAT (//,' DEPTH',10(2X,A10))
      IF (N.LT.J2) PRINT 5
    5 FORMAT (1X)
      PRINT 5

      DO 3 L=1,ND
      IF(((L-1)/LSPOP)*LSPOP.NE.(L-1) .AND. L.NE.ND) GOTO 3
      ENCODE (130,6,PRILINE) L
    6 FORMAT (I6)

      DO 12 J=J1,J2
      IF (POP1(L,J) .NE. .0) THEN
            ENCODE (12,11,NUMBER) POPNUM(L,J)/POP1(L,J)
   11       FORMAT (F12.4)
            ELSE
            NUMBER='    INFINITE'
            ENDIF
      I=7+(J-J1)*12
      PRILINE(I:I+11)=NUMBER
   12 CONTINUE

      PRINT 13,PRILINE
   13 FORMAT (A)
    3 CONTINUE

      IF (J2.EQ.N) GOTO 8
      J1=J1+10
      GOTO 4
    8 CONTINUE

!----------------------------------------------------------------------------------
! CALCULATION BLOCK START
!----------------------------------------------------------------------------------

      DEPTHS_DIM = ND - 3

      ALLOCATE(DEPTHS(DEPTHS_DIM))

      DEPTH_IND = 1

      DO L = 3, ND - 1; DEPTHS(DEPTH_IND) = L; DEPTH_IND = DEPTH_IND + 1; ENDDO

      CORMAX_ALL = 0.0D0
      CORMAX_NHE = 0.0D0
      CORMAX_NHN = 0.0D0

      DCM_ALL = 0
      DCM_NHE = 0
      DCM_NHN = 0

      DO NA = 1, NATOM

         CALL CALC_CONV_IND(NFIRST(NA), NLAST(NA), ELEMENT(NA), DEPTHS, DEPTHS_DIM, POPNUM, POP1, LEVEL, N, ND, CM, DCM)

         CORMAX_ELE(NA) = CM

         DCM_ELE(NA) = DCM

         ALLE = CORMAX_ELE(NA) .GT. CORMAX_ALL
         NOHE = CORMAX_ELE(NA) .GT. CORMAX_NHE .AND. TRIM(ADJUSTL(ELEMENT(NA))) .NE. 'HELIUM'

         IF (ALLE) THEN; CORMAX_ALL = CORMAX_ELE(NA); DCM_ALL = DCM_ELE(NA); ENDIF

         IF (NOHE) THEN; CORMAX_NHE = CORMAX_ELE(NA); DCM_NHE = DCM_ELE(NA); ENDIF

         DO J = NFIRST(NA), NLAST(NA)

            CALL CALC_CONV_IND(J, J, LEVEL(J), DEPTHS, DEPTHS_DIM, POPNUM, POP1, LEVEL, N, ND, CM, DCM)

            CORMAX_LEV(J) = CM

            NOHN = CORMAX_LEV(J) .GT. CORMAX_NHN .AND.
     $             TRIM(ADJUSTL(ELEMENT(NA))) .NE. 'HELIUM' .AND.
     $             LEVEL(J) .NE. 'NeII-----1'

            IF (NOHN) THEN; CORMAX_NHN = CORMAX_LEV(J); DCM_NHN = DCM; ENDIF

         ENDDO

      ENDDO

!      CORMAX_FIN = CORMAX_ELE(1); DCM_FIN = DCM_ELE(1)
      CORMAX_FIN = CORMAX_NHN; DCM_FIN = DCM_NHN
!      CORMAX_FIN = CORMAX_ALL; DCM_FIN = DCM_ELE(1)

!----------------------------------------------------------------------------------
! PRINT OUT BLOCK START
!----------------------------------------------------------------------------------

      WRITE(*, '(/, A)') '---------------------------------------------------------
     $--------------------------------------------------------'

      WRITE(*, 10002) DEPTHS(1), DEPTHS(DEPTHS_DIM), CORMAX_ELEC, LOG10(CORMAX_ELEC)

      WRITE(*, '(/, A)') '---------------------------------------------------------
     $--------------------------------------------------------'

!==========================================================================================

      IF (LAMBDA_ITER .EQ. 1) THEN; CALL MKDIR(CONV_DIR); CALL CLEAN_DIR(CONV_DIR); ENDIF

      DO NA = 1, NATOM

         CONV_FILE = CONV_DIR//ELEMENT(NA)

         CALL OPEN_TO_APPEND(999 + NA, CONV_FILE)

         IF (LAMBDA_ITER .EQ. 1) THEN

             WRITE(999 + NA, '(3x, A, 5x, A, 6x, A, 1x, $)') 'LI', 'MAX', 'DI'

             DO J = NFIRST(NA), NLAST(NA); WRITE(999 + NA, '(2x, A10, $)') LEVEL(J); ENDDO

             WRITE(999 + NA, '(/)')

         ENDIF

         WRITE(999 + NA, '(I5, 2x, ES9.3, 2x, I3, $)') LAMBDA_ITER, CORMAX_ELE(NA), DCM_ELE(NA)

         DO J = NFIRST(NA), NLAST(NA)

            IF (J .NE. NLAST(NA)) WRITE(999 + NA, '(3x, ES9.3, $)') CORMAX_LEV(J)
            IF (J .EQ. NLAST(NA)) WRITE(999 + NA, '(3x, ES9.3)')    CORMAX_LEV(J)

         ENDDO

         CLOSE(UNIT = 999 + NA)

      ENDDO

!==========================================================================================

      CALL OPEN_TO_APPEND(1000 + NATOM, CONV_DIR//'ELEM')

      IF (LAMBDA_ITER .EQ. 1) THEN

         WRITE(1000 + NATOM, '(3x,A,$)') 'LI'

         DO NA = 1, NATOM

            ELE = ELEMENT(NA)

            IF (NA .EQ. 1) WRITE(1000 + NATOM, '(5x,A,$)') ELE(1 : 3)

            IF (NA .NE. 1) WRITE(1000 + NATOM, '(8x,A,$)') ELE(1 : 3)

         ENDDO

         WRITE(1000 + NATOM, '(8x,A,/)') 'ELE'

      ENDIF

      WRITE(1000 + NATOM, '(I5,$)') LAMBDA_ITER

      DO NA = 1, NATOM; WRITE(1000 + NATOM, '(2x,ES9.3,$)') CORMAX_ELE(NA); ENDDO

      WRITE(1000 + NATOM, '(2x,ES9.3)'), CORMAX_ELEC

      CLOSE(UNIT = 1000 + NATOM)

!==========================================================================================

      CALL OPEN_TO_APPEND(1001 + NATOM, CONV_DIR//'ALL')

      IF (LAMBDA_ITER .EQ. 1) THEN

         WRITE(1001 + NATOM, '(3x,A,7x,A,6x,A,3(8x,A),5x,A,/)') 'LI', 'TIME', 'ALL', 'NHE', 'NHN', 'HYD', 'DCMF'

      ENDIF

      time = trim(adjustl(writeTOC()))

      if (len_trim(time) == 5) time = '00:'//time

      WRITE(1001 + NATOM, '(I5,5x,A8,4(2x,ES9.3),2x,I3)') LAMBDA_ITER, time, CORMAX_ALL,
     $                                                    CORMAX_NHE, CORMAX_NHN, CORMAX_ELE(1), DCM_FIN

!----------------------------------------------------------------------------------
!PRINT OUT BLOCK END
!----------------------------------------------------------------------------------

      DEALLOCATE(DEPTHS)

      IF (GAMMAL.GT..0)
     $   PRINT 14, GAMMAL,DELTAC,GAMMAR
   14    FORMAT (/,10X,'SCHARMER-PARAMETERS:',  10X,
     $          '   GAMMAL=',F5.1,/,
     $        40X,'DELTAC=',F5.1,'   GAMMAR=',F7.1)

      IF (REDUCE .NE. 1.) PRINT 10,REDUCE
   10 FORMAT (/,10X,'CORRECTIONS REDUCED BY FACTOR',F5.2)

C***  INHIBIT NEGATIVE POP. NUMBERS
      NEGWARN=0
      DO 19 L=1,ND
      if (.NOT. allocated(ABXYZ_new)) allocate(ABXYZ_new(NATOM))
      ABXYZ_new(1:NATOM)=ABXYZn(1:NATOM,L)
      NEGPOPL=.FALSE.
      DO 16 J=1,N
      IF (POPNUM(L,J) .LE. .0) THEN
            NEGPOPL=.TRUE.
            NEGWARN=NEGWARN+1
            POPNUM(L,J)=0.5*ABS(POP1(L,J))
      ENDIF
   16 CONTINUE

      IF (NEGPOPL) THEN
C***     RENORMALIZATION OF THE POPULATION NUMBERS
         DO 15 NA=1,NATOM
         NFIRNA=NFIRST(NA)
         NLANA=NLAST(NA)
         SUM=0.0
         DO 18 J=NFIRNA,NLANA
   18    SUM=SUM+POPNUM(L,J)
         SUM=SUM/ABXYZ_new(NA)
         DO 15 J=NFIRNA,NLANA
         POPNUM(L,J)=POPNUM(L,J)/SUM
   15    CONTINUE
C***     RE-ADJUST THE ELECTRON DENSITY
         RNEL=0.0
         DO 20 J=1,N
         RNEL=RNEL+NCHARG(J)*POPNUM(L,J)
   20    CONTINUE
         RNE(L)=RNEL
      ENDIF
      if (allocated(ABXYZ_new)) deallocate(ABXYZ_new)
   19 CONTINUE

      IF (NEGWARN .GT. 0) PRINT 17,NEGWARN
   17 FORMAT (/,10X,'WARNING: NEGATIVE POP. NUMBERS ARE SET TO 0.5 ',
     $ 'TIMES THEIR OLD VALUE:',I4,' TIMES')
C***  PRINTOUT OF ERROR MESSAGES FROM SUBROUTINE GAUNTFF
      IF (LSPOP .EQ. 1) THEN
      IF (NTUP.GT.0)
     $PRINT  *, NTUP,' CALLS OF GAUNTFF BEYOND UPPER TEMPERATURE BOUND'
      IF (NTLOW.GT.0)
     $PRINT *,NTLOW,' CALLS OF GAUNTFF BEYOND LOWER TEMPERATURE BOUND'
      IF (NFUP.GT.0)
     $PRINT *,NFUP,' CALLS OF GAUNTFF BEYOND UPPER FREQUENCY BOUND'
      IF (NFLOW.GT.0)
     $PRINT *,NFLOW,' CALLS OF GAUNTFF BEYOND LOWER FREQUENCY BOUND'
      ENDIF

C***  PRINTOUT OF ERROR MESSAGES FROM SUBROUTINE LINPOP
      IF (NEGINTL .GT. 0)
     $  PRINT *,'           WARNING: NEGATIVE LINE INTENSITIES ',
     $      'REPLACED BY THEIR ABSOLUTE VALUE: ',NEGINTL,' TIMES'

      if (allocated(ABXYZ_new)) deallocate(ABXYZ_new)

      RETURN

10002 FORMAT(/,15x,'DEPTHS: ',I2,' TO ',I3,5X,'CORMAX= ',ES9.3,5X,'LOG=',F5.2,
     $         5X,'ELEMENT: ELECTRONS')

      END SUBROUTINE

      SUBROUTINE CALC_CONV_IND(NFIRNA, NLANA, ELEMENT, DEPTHS, DEPTHS_DIM,
     $                         POP, PREV_POP, LEVEL, N, ND, CORMAX, DCM)

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: DEPTHS_DIM, NFIRNA, NLANA

      INTEGER, DIMENSION(DEPTHS_DIM), INTENT(IN) :: DEPTHS

      CHARACTER(LEN = *), INTENT(IN) :: ELEMENT

      INTEGER, INTENT(IN) :: N, ND

      REAL*8, DIMENSION(ND, N), INTENT(IN) :: POP, PREV_POP

      CHARACTER (LEN = 10), DIMENSION(N), INTENT(IN) :: LEVEL

      REAL*8,  INTENT(OUT) :: CORMAX

      INTEGER, INTENT(OUT) :: DCM

      CHARACTER (LEN = 10) :: LEVMAX, LEVMIN, LEVMA2, LEVMI2

      REAL*8 :: Q, QMAX, QMA2, QMIN, QMI2

      INTEGER :: LMAX, LMA2, LMIN, LMI2

      REAL*8 :: CORLOG

      INTEGER :: I, J, DEPTH_IND

C***  CALCULATE MINIMUM AND MAXIMUM OF ALL CORRECTION FACTORS
      LMAX = 0
      LMA2 = 0
      LMIN = 0
      LMI2 = 0
      QMAX = 1.0D0
      QMA2 = 1.0D0
      QMIN = 1.0D0
      QMI2 = 1.0D0
      LEVMAX = 'NONE'
      LEVMA2 = 'NONE'
      LEVMIN = 'NONE'
      LEVMI2 = 'NONE'

      DCM = 0
      CORMAX = 0.0D0

      DEPTH: DO I = 1, DEPTHS_DIM

              DEPTH_IND = DEPTHS(I)

         LEVELS: DO J = NFIRNA, NLANA

                 IF (PREV_POP(DEPTH_IND, J) .EQ. 0.0D0) CYCLE LEVELS

                 Q = POP(DEPTH_IND, J) / PREV_POP(DEPTH_IND, J)

                 IF (Q .GT. QMA2) THEN

                    IF (Q .GT. QMAX) THEN

                       QMA2 = QMAX
                       QMAX = Q
                       LEVMA2 = LEVMAX
                       LEVMAX = LEVEL(J)
                       LMA2 = LMAX
                       LMAX = DEPTH_IND

                    ELSE

                       QMA2 = Q
                       LEVMA2 = LEVEL(J)
                       LMA2 = DEPTH_IND

                    ENDIF

                 ENDIF

                 IF (Q .LT. QMI2) THEN

                    IF (Q .LT. QMIN) THEN

                       QMI2 = QMIN
                       QMIN = Q
                       LEVMI2 = LEVMIN
                       LEVMIN = LEVEL(J)
                       LMI2 = LMIN
                       LMIN = DEPTH_IND

                    ELSE

                       QMI2 = Q
                       LEVMI2 = LEVEL(J)
                       LMI2 = DEPTH_IND

                    ENDIF

                 ENDIF

         ENDDO LEVELS

      ENDDO DEPTH

      CORMAX = DMAX1(QMAX - 1.0D0, 1.0D0 - QMIN)

      IF (ABS(CORMAX - QMAX + 1) .LT. 1D-15) DCM = LMAX

      IF (ABS(CORMAX + QMIN - 1) .LT. 1D-15) DCM = LMIN

      IF (CORMAX .GT. 0.0D0) CORLOG = LOG10(CORMAX)

      IF (NFIRNA .NE. NLANA) THEN

         WRITE(*, '(/,A)') '---------------------------------------------------------
     $--------------------------------------------------------'

         WRITE(*, 9) QMAX, LEVMAX, LMAX, QMA2, LEVMA2, LMA2, QMIN, LEVMIN, LMIN, QMI2, LEVMI2, LMI2,
     $            DEPTHS(1), DEPTHS(DEPTHS_DIM), CORMAX, CORLOG, ELEMENT

      ENDIF

    9 FORMAT(/,15X, 'QMAX1: ', ES9.3,'  (', A10,'  L=', I3,')',
     $         4X,  'QMAX2: ', ES9.3,'  (', A10,'  L=', I3,')',
     $       /,15X, 'QMIN1: ', ES9.3,'  (', A10,'  L=', I3,')',
     $         4X,  'QMIN2: ', ES9.3,'  (', A10,'  L=', I3,')', /,
     $       /,15X, 'DEPTHS: ', I2, ' TO ', I3,
     $         5X,  'CORMAX =',ES9.3,5X,'LOG(CORMAX) = ',F5.2,
     $         5X,  'ELEMENT: ', A)

      END SUBROUTINE

      END MODULE

ext: for

max_stars_repo_name: rtagirov/nessy

max_issues_repo_name: rtagirov/nessy

max_forks_repo_name: rtagirov/nessy-src


text:
      PROGRAM SCALAR01
      REAL X
      PRINT *,X
      END