MODULE fast_subs_univ

implicit none

CONTAINS

!*****************************************************
! Commonly used subroutines for fast method
! written by Norbert Schorghofer 2007-2011
!*****************************************************

pure function zint(y1,y2,z1,z2)
  ! interpolate between two values, y1*y2<0
  implicit none
  real(8), intent(IN) :: y1,y2,z1,z2
  real(8) zint
  zint = (y1*z2 - y2*z1)/(y1-y2)
end function zint



pure function equildepth(nz, z, rhosatav, rhosatav0, avrho1)
!***********************************************************************
!  returns equilibrium depth for given ice content
!  this is not the true (final) equilibrium depth
!***********************************************************************
  implicit none
  integer, intent(IN) :: nz
  real(8), intent(IN) :: z(nz), rhosatav(nz)
  real(8), intent(IN) :: rhosatav0, avrho1
  integer i, typeE
  real(8) equildepth  ! =zdepthE
  !real(8), external :: zint  ! defined in allinterfaces.mod
  
  typeE = -9; equildepth = -9999.
  do i=1,nz 
     if (rhosatav(i) <= avrho1) then
        typeE=i
        exit
     endif
  enddo
  if (typeE>1) then ! interpolate
     equildepth=zint(avrho1-rhosatav(typeE-1),avrho1-rhosatav(typeE),z(typeE-1),z(typeE))
  end if
  if (typeE==1) equildepth=zint(avrho1-rhosatav0,avrho1-rhosatav(1),0.d0,z(1))
  if (equildepth>z(nz)) equildepth=-9999.   ! happens very rarely
end function equildepth



subroutine depths_avmeth(typeT, nz, z, rhosatav, rhosatav0, rlow, avrho1,  &
     & porefill, typeF, zdepthF, B, ypp, typeG, zdepthG)
!***********************************************************************
!  returns interface depth and ypp
!  also returns lower geothermally limited boundary, if applicable
!
!  this is a crucial subroutine
!
!  B = Diff*bigstep/(porosity*icedensity)  [SI units]
!***********************************************************************
  use math_mod, only: deriv2_simple, deriv1_onesided, deriv1, colint
  use ice_table_mod, only: constriction
  implicit none
  integer, intent(IN) :: nz, typeT
  real(8), intent(IN), dimension(nz) :: z, rhosatav, porefill
  real(8), intent(IN) :: rhosatav0, rlow, avrho1
  integer, intent(OUT) :: typeF  ! index of depth below which filling occurs
  real(8), intent(INOUT) :: zdepthF
  real(8), intent(IN) :: B 
  real(8), intent(OUT) :: ypp(nz), zdepthG
  integer, intent(INOUT) :: typeG  ! positive on input when Fgeotherm>0

  integer i, typeP, nlb, newtypeG
  real(8) eta(nz), Jpump1, help(nz), yp(nz), zdepthFold, ap_one, ap(nz)
  real(8) leak, cumfill, cumfillabove

  if (typeT<0) then
     nlb = nz
     do i=1,nz
        eta(i) = constriction(porefill(i))
     enddo
  else
     !nlb = typeT-1
     nlb = typeT ! changed 2010-09-29
     do i=1,typeT-1
        eta(i) = constriction(porefill(i))
     enddo
     do i=typeT,nz
        eta(i)=0.
     enddo
  end if

!-fill depth
  zdepthFold = zdepthF
  typeF = -9;  zdepthF = -9999.
  call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp)  ! yp also used below
  do i=1,nlb
     Jpump1 = (rhosatav(i)-avrho1)/z(i)  ! <0 when stable
     ! yp is always <0
     help(i) = Jpump1 - eta(i)*yp(i)
     leak = porefill(i)/B*(z(i)-zdepthFold)/(18./8314.)
     !help(i) = help(i)-leak   ! optional
     if (help(i) <= 0.) then
        typeF=i
        !print *,'#',typeF,Jpump1,eta(typeF)*yp(typeF),leak
        exit
     endif
  enddo
  if (typeF>1) zdepthF = zint(help(typeF-1),help(typeF),z(typeF-1),z(typeF))
  if (typeF==1) zdepthF=z(1)


!-depth to shallowest perennial ice
  typeP = -9 
  do i=1,nz
     if (porefill(i)>0.) then
        typeP = i  ! first point with ice
        exit
     endif
  enddo

!-calculate ypp
  !call deriv1(z,nz,rhosatav,rhosatav0,rlow,yp)
  call deriv1(z,nz,eta(:),1.d0,eta(nz-1),ap)
  if (typeP>0 .and. typeP<nz-2) then
     call deriv1_onesided(typeP,z,nz,eta(:),ap_one)
     ! print *,typeP,ap(typeP),ap_one
     ap(typeP)=ap_one
  endif
  call deriv2_simple(z,nz,rhosatav(1:nz),rhosatav0,rlow,ypp(:))
  !call deriv2_full(z,nz,eta(:),rhosatav(:),1.d0,rhosatav0,rhosatav(nz-1),ypp(:))
  !write(40,*) rhosatav
  !write(41,*) yp
  !write(42,*) ypp
  ypp(:) = ap(:)*yp(1:)+eta(:)*ypp(:)
  !write(43,*) ypp
  !write(44,*) eta(1:nz)
  !write(45,*) (rhosatav(:)-avrho1)/z(:)
  ypp(:) = ypp(:)*18./8314.
  ! ypp values beyond nlb should not be used

!-geothermal stuff
  if (typeT>0) typeG=-9
  if (typeG<0) zdepthG=-9999.
  if (typeG>0 .and. typeT<0) then
     typeG=-9
     do i=2,nz
        if (yp(i)>0.) then  ! first point with reversed flux
           typeG=i
           zdepthG=zint(yp(i-1),yp(i),z(i-1),z(i))
           !zdepthG=z(typeG)
           exit
        endif
     enddo
  else
     typeG = -9
  endif
  if (typeG>0 .and. typeT<0) then
     call colint(porefill(:)/eta(:),z,nz,typeG-1,nz,cumfillabove) 
     newtypeG = -9
     do i=typeG,nz
        if (minval(eta(i:nz))<=0.) cycle  ! eta=0 means completely full
        call colint(porefill(:)/eta(:),z,nz,i,nz,cumfill)
        if (cumfill<yp(i)*18./8314.*B) then  ! usually executes on i=typeG
           if (i>typeG) then
              write(34,*) '# adjustment to geotherm depth by',i-typeG
              zdepthG = zint(yp(i-1)*18./8314.*B-cumfillabove, &  ! no good
                   &        yp(i)*18./8314.*B-cumfill,z(i-1),z(i))
              if (zdepthG>z(i) .or. zdepthG<z(i-1)) write(34,*) &
                   & '# WARNING: zdepthG interpolation failed',i,z(i-1),zdepthG,z(i)
              newtypeG=i
           endif
           ! otherwise leave zdepthG untouched
           exit
        endif
        cumfillabove = cumfill
     enddo
     if (newtypeG>0) typeG=newtypeG
  end if
  ! if typeG>0, then all ice at and below typeG should be erased 
end subroutine depths_avmeth



pure subroutine icechanges(nz,z,typeF,avdrho,avdrhoP,ypp, &
     & Diff,porosity,icefrac,bigstep,zdepthT,porefill,typeG)
!***********************************************************
! advances ice table, advances interface, grows pore ice
!
! a crucial subroutine
!***********************************************************
  use math_mod, only: colint
  use ice_table_mod, only: rho_ice
  implicit none
  integer, intent(IN) :: nz, typeF, typeG
  real(8), intent(IN) :: z(nz), ypp(nz), avdrho, avdrhoP
  real(8), intent(IN) :: Diff, porosity, icefrac, bigstep
  real(8), intent(INOUT) :: zdepthT, porefill(nz)
  integer j, erase, newtypeP, ub, typeP, typeT
  real(8) B, beta, integ

  B = Diff*bigstep*86400.*365.24/(porosity*927.)
  !B = Diff*bigstep*86400.*365.24/(porosity*rho_ice(T,'h2o'))

  ! advance ice table, avdrho>0 is retreat
  if (zdepthT>=0. .and. avdrho>0.) then 
     typeP=-9999; typeT=-9999
     do j=1,nz
        if (z(j)>zdepthT) then
           typeT=j
           exit
        endif
     enddo
     do j=1,nz
        if (porefill(j)>0.) then
           typeP=j
           exit
        endif
     enddo
     if (typeP==typeT) then   ! new 2011-09-01
        beta = (1-icefrac)/(1-porosity)/icefrac
        beta = Diff*bigstep*beta*86400*365.24/927.
        !beta = Diff*bigstep*beta*86400*365.24/rho_ice(T,'h2o')
        zdepthT = sqrt(2*beta*avdrho*18./8314. + zdepthT**2)
     endif
  endif
  if (zdepthT>z(nz)) zdepthT=-9999.
  
  ! advance interface, avdrhoP>0 is loss from zdepthP
  if (avdrhoP>0.) then
     erase=0
     do j=1,nz
        if (typeF>0 .and. j>=typeF) exit ! don't retreat beyond typeF
        if (zdepthT>=0. .and. z(j)>zdepthT) exit 
        call colint(porefill(1:nz)*z(1:nz),z(1:nz),nz,1,j,integ)
        erase = j
        if (integ>B*avdrhoP*18./8314.) exit
     end do
     if (erase>0) porefill(1:erase)=0.
  endif

  ! new depth
  newtypeP = -9 
  do j=1,nz
     if (zdepthT>=0. .and. z(j)>zdepthT) exit
     if (porefill(j)>0.) then
        newtypeP = j  ! first point with pore ice
        exit
     endif
  enddo

  ! diffusive filling
  ub = typeF
  if (newtypeP>0 .and. typeF>0 .and. newtypeP<ub) ub=newtypeP
  if (ub>0) then  
     do j=ub,nz
        porefill(j) = porefill(j) + B*ypp(j)
        if (porefill(j)<0.) porefill(j)=0.
        if (porefill(j)>1.) porefill(j)=1.
        if (zdepthT>=0. .and. z(j)>zdepthT) exit
     enddo
  end if

  ! below icetable
  if (zdepthT>=0.) then
     do j=1,nz
        if (z(j)>zdepthT) porefill(j) = icefrac/porosity
     enddo
  else
     ! geothermal lower boundary
     if (typeG>0) porefill(typeG:nz)=0.
  end if
end subroutine icechanges

END MODULE fast_subs_univ