        MODULE SOIL_MOD

        use LAKE_DATATYPES, only : ireals, iintegers
        use NUMERICS, only : PROGONKA, STEP
        use NUMERIC_PARAMS, only : vector_length

        contains
        SUBROUTINE COMSOILFORLAKE

        !COMSOILFORLAKE specifies parameters of soil according to soil type

        use NUMERIC_PARAMS
        use DRIVING_PARAMS
        use ARRAYS
        use ARRAYS_SOIL
        implicit none

        real(kind=ireals) :: b(1:Num_Soil)
        real(kind=ireals) :: psi_max(1:Num_Soil)
        real(kind=ireals) :: Porosity(1:Num_Soil)
        real(kind=ireals) :: gamma_max(1:Num_Soil)
        real(kind=ireals) :: lambda_max(1:Num_Soil)
        real(kind=ireals) :: W_0(1:Num_Soil)
        real(kind=ireals) :: W_m(1:Num_Soil)
        real(kind=ireals) :: pow
          
        integer(kind=iintegers) :: i, j
        logical :: firstcall

        data firstcall /.true./

        SAVE
        if (firstcall) then
        ! allocate (AL(1:ns+2),DLT(1:ns+2),DVT(1:ns+2),DL(1:ns+2),
        !& ALV(1:ns+2),DV(1:ns+2),Z(1:ns+2),T(1:ns+2),WL(1:ns+2),
        !& WV(1:ns+2),WI(1:ns+2),dens(1:ns+2))
        endif

        !I=1 ! SAND
        !I=2 ! LOAMY SAND
        !I=3 ! SANDY LOAM
        !I=4 ! LOAM
        !I=5 ! SILT LOAM
        !I=6 ! SANDY CLAY LOAM
        !I=7 ! CLAY LOAM
        !I=8 ! SILTY CLAY LOAM
        !I=9 ! SANDY CLAY
        !I=10! SILTY CLAY
        !I=11! CLAY
        !----------------------------------------------------------------
        !  NN    b      psi_max  Por   gamma_max lambda_max W_0     W_m
        !----------------------------------------------------------------
        !   1   4.05    3.50    0.395   0.01760   0.29800   0.02    0.01
        !   2   4.38    1.78    0.410   0.01560   0.13900   0.05    0.02
        !   3   4.90    7.18    0.435   0.00340   0.13700   0.08    0.03
        !   4   5.39    14.6    0.451   0.00069   0.06190   0.20    0.08
        !   5   5.30    56.6    0.485   0.00072   0.20400   0.18    0.07
        !   6   7.12    8.63    0.420   0.00063   0.03070   0.13    0.06
        !   7   8.52    36.1    0.476   0.00024   0.03720   0.27    0.12
        !   8   7.75    14.6    0.477   0.00017   0.01240   0.24    0.11
        !   9   10.4    6.16    0.426   0.00021   0.00641   0.23    0.10
        !  10   10.4    17.4    0.492   0.00010   0.00755   0.30    0.15
        !  11   11.4    18.6    0.482   0.00013   0.00926   0.40    0.20
        !-----------------------------------------------------------------

        zsoil(1) = 0.

        if (UpperLayer > 0.) then
          pow = log(UpperLayer/depth%par)/log(1.d0/float(ns-1))
          do i = 2, ns
            zsoil(i) = (1.*(i-1)/(ns-1))**pow*depth%par
          end do
        else
          do i = 2, ns
            zsoil(i) = (1.*(i-1)/(ns-1))*depth%par
          end do
        end if

        do i = 1, ns-1
          dzs(i) = zsoil(i+1) - zsoil(i)
        end do

        do i = 1, ns
          if (i == 1) then
            dzss(i) = 0.5d0*dzs(i)
          elseif (i == ns) then
            dzss(i) = 0.5d0*dzs(i-1)
          else
            dzss(i) = 0.5d0*(dzs(i-1) + dzs(i))
          endif 
        end do 

        !SoilType = 7

        b(1) = 4.05
        b(2) = 4.38
        b(3) = 4.90
        b(4) = 5.39
        b(5) = 5.30
        b(6) = 7.12
        b(7) = 8.52
        b(8) = 7.75
        b(9) = 10.4
        b(10)= 10.4
        b(11)= 11.4

        psi_max( 1) = 3.50/100. 
        psi_max( 2) = 1.78/100. 
        psi_max( 3) = 7.18/100. 
        psi_max( 4) = 14.6/100. 
        psi_max( 5) = 56.6/100. 
        psi_max( 6) = 8.63/100. 
        psi_max( 7) = 36.1/100. 
        psi_max( 8) = 14.6/100. 
        psi_max( 9) = 6.16/100. 
        psi_max(10) = 17.4/100. 
        psi_max(11) = 18.6/100. 

        Porosity( 1) = 0.395
        Porosity( 2) = 0.410
        Porosity( 3) = 0.435
        Porosity( 4) = 0.451
        Porosity( 5) = 0.485
        Porosity( 6) = 0.420
        Porosity( 7) = 0.476
        Porosity( 8) = 0.477
        Porosity( 9) = 0.426
        Porosity(10) = 0.492
        Porosity(11) = 0.482

        gamma_max( 1) = 0.01760/100. 
        gamma_max( 2) = 0.01560/100. 
        gamma_max( 3) = 0.00340/100. 
        gamma_max( 4) = 0.00069/100. 
        gamma_max( 5) = 0.00072/100. 
        gamma_max( 6) = 0.00063/100. 
        gamma_max( 7) = 0.00024/100. 
        gamma_max( 8) = 0.00017/100. 
        gamma_max( 9) = 0.00021/100. 
        gamma_max(10) = 0.00010/100. 
        gamma_max(11) = 0.00013/100. 

        lambda_max( 1) = 0.29800/10000. 
        lambda_max( 2) = 0.13900/10000.
        lambda_max( 3) = 0.13700/10000.
        lambda_max( 4) = 0.06190/10000.
        lambda_max( 5) = 0.20400/10000.
        lambda_max( 6) = 0.03070/10000.
        lambda_max( 7) = 0.03720/10000.
        lambda_max( 8) = 0.01240/10000.
        lambda_max( 9) = 0.00641/10000.
        lambda_max(10) = 0.00755/10000.
        lambda_max(11) = 0.00926/10000.

        W_0( 1) = 0.02
        W_0( 2) = 0.05
        W_0( 3) = 0.08
        W_0( 4) = 0.20
        W_0( 5) = 0.18
        W_0( 6) = 0.13
        W_0( 7) = 0.27
        W_0( 8) = 0.24
        W_0( 9) = 0.23
        W_0(10) = 0.30
        W_0(11) = 0.40

        W_m( 1) = 0.01
        W_m( 2) = 0.02
        W_m( 3) = 0.03
        W_m( 4) = 0.08
        W_m( 5) = 0.07
        W_m( 6) = 0.06
        W_m( 7) = 0.12
        W_m( 8) = 0.11
        W_m( 9) = 0.10
        W_m(10) = 0.15
        W_m(11) = 0.20

        if (SoilType%par > 0) then
          do j = 1, ns
            BH(j)= b(SoilType%par)   ! PARAMETER B, DIMENSOINLESS
            PSIMAX(j) = - psi_max(SoilType%par) ! SAT. WATER POTENTIAL, M.    
            POR(j)    = Porosity(SoilType%par)  ! POROSITY, DIMENSIONLESS
            FLWMAX(j) = gamma_max(SoilType%par) ! SAT. HYDR. CONDUCTIVITY, M/S
            DLMAX(j)  = lambda_max(SoilType%par)! SAT. WATER DIFFUSIVITY, ?KG/(M*SEC)
            WLM0(j)   = W_0(SoilType%par) ! MAXIMAL UNFREEZING WATER AT 0C
            WLM7(j)   = W_m(SoilType%par) ! MAXIMAL UNFREEZING WATER AT T<<0C
          end do
        else
          write(*,*) 'LAKE: Soil type is not given: STOP'
          STOP
        !do j = 1, ns
        ! BH(j)= BH_soil! PARAMETER B, DIMENSOINLESS
        ! PSIMAX(j) = PSIMAX_soil ! SAT. WATER POTENTIAL, M.    
        ! POR(j)    = POR_soil    ! POROSITY, DIMENSIONLESS
        ! FLWMAX(j) = FLWMAX_soil ! SAT. HYDR. CONDUCTIVITY, M/S
        ! DLMAX(j)  = DLMAX_soil  ! SAT. WATER DIFFUSIVITY, ?KG/(M*SEC)
        ! WLM0(j)   = WLM0_soil   ! MAXIMAL UNFREEZING WATER AT 0C
        ! WLM7(j)   = WLM7_soil   ! MAXIMAL UNFREEZING WATER AT T<<0C
        !end do
        end if

        rosdry(:) = 1.2E+3

        if (firstcall) firstcall=.false.   
        RETURN
        END SUBROUTINE COMSOILFORLAKE


        SUBROUTINE SOILFORLAKE(dt,a,b,c,d,is)

        !SOILFORLAKE calculates profiles of temperature, liquid and solid water 
        !content in the soil under a lake

        use PHYS_CONSTANTS 
        use METH_OXYG_CONSTANTS!, only : &
        !& methhydrdiss, &
        !& alphamh, &
        !& nch4_d_nh2o, &
        !& molmass_h2o
        use DRIVING_PARAMS
        use ARRAYS
        use ARRAYS_SOIL
        use ARRAYS_METHANE
        use ARRAYS_BATHYM
        use PHYS_FUNC!, only : &
        !& MELTPNT, &
        !& TEMPMHYDRDISS, &
        !& UNFRWAT, &
        !& WL_MAX, &
        !& MELTINGPOINT
        use ATMOS!, only : &
        !& pressure

        implicit none

        real(kind=ireals),    parameter :: pi = 3.141592653589d0

        real(kind=ireals), intent(in) :: dt

        !integer(kind=iintegers), intent(in) :: ix, iy
        !integer(kind=iintegers), intent(in) :: nx, ny
        integer(kind=iintegers), intent(in) :: is ! Number of soil column 

        !integer(kind=iintegers), intent(in) :: year, month, day

        !real(kind=ireals)   , intent(in) :: hour
        !real(kind=ireals)   , intent(in) :: phi
        !real(kind=ireals)   , intent(in) :: extwat, extice
        !real(kind=ireals)   , intent(in) :: fetch
            
        real(kind=ireals), intent(inout) :: a(1:vector_length)
        real(kind=ireals), intent(inout) :: b(1:vector_length)
        real(kind=ireals), intent(inout) :: c(1:vector_length)
        real(kind=ireals), intent(inout) :: d(1:vector_length)
        !real(kind=ireals) :: Temp(1:vector_length)

        real(kind=ireals) :: dzmean
        real(kind=ireals) :: ma, mi
        real(kind=ireals) :: wflow, wfhigh
        real(kind=ireals) :: lammoist1, lammoist2
        real(kind=ireals) :: surplus
        real(kind=ireals) :: Potphenergy
        real(kind=ireals) :: watavailtofreeze
        real(kind=ireals) :: filtr_low
        real(kind=ireals) :: dhwfsoil1, dhwfsoil2, dhwfsoil3, dhwfsoil4
        real(kind=ireals) :: max1
        real(kind=ireals) :: work, mhsep
        real(kind=ireals), pointer :: hicemelt

        real(kind=ireals) :: alp(10)
        real(kind=ireals) :: psiit(10)

        real(kind=ireals), allocatable :: pressoil(:)

        integer(kind=iintegers) :: i, j, k
        integer(kind=iintegers) :: iter, iter3

        logical :: firstcall

        data psiit/6,3,7,2,5,4,8,1,0,0/ !/3,2,4,1,0,0,0,0,0,0/  
        data firstcall /.true./

        SAVE

        !open (133, file=dir_out//'\dhwfsoil.dat', status = 'unknown')   

        !PARAMETERS OF ITERATIONAL PROCESS
        ma = 10.5 !5.5
        mi = 4.5
        do i = 1, 8 
          alp(i) = 2*(mi+ma-(ma-mi)*cos((2*psiit(i)-1)/16*pi))**(-1)
        enddo

        allocate (pressoil(1:ns))
        do i = 1, ns
          pressoil(i) = pressure + row0*g*(h1 + zsoil(i))
        enddo

        if (tricemethhydr%par > 0.) then
          hicemelt => methhydrdiss
        else
          hicemelt => Lwi
        endif
           
        !do i=1, ns
        ! wlmax_i = POR(i)*ROW0/(rosdry(i)*(1-por(i))) 
        ! BB1 = BH(i) + 2.
        ! csoil(i) = cr+WL1(i)*CW+WI1(i)*CI
        ! rosoil(i) = rosdry(i)*(1-por(i))*(1+wi1(i)+wl1(i)) 
        ! ARG = (WL1(i)+WI1(i))/wlmax_i
        ! ARG = max(ARG,1.d-2)
        ! PSI = PSIMAX(i)*(ARG)**(-BH(i))
        ! PF = log(-PSI*100.)/log(10.d0)
        ! IF(PF>5.1) THEN
        !  lamsoil(i) = 4.1E-4*418.
        ! ELSE
        !  lamsoil(i) = exp(-PF-2.7)*418.
        ! END IF
        ! wlmax_i = POR(i)*ROW0/(rosdry(i)*(1-por(i))) - wi1(i)*row0_d_roi
        ! ARG = (WL1(i))/wlmax_i
        ! ARG = max(ARG,1.d-2)
        ! lammoist(i) = dlmax(i)*ARG**BB1
        !end do


        !do i=1, ns-1
        ! wlmax_i=WL_MAX(por(i+1),rosdry(i+1),wi1(i+1))
        ! if (wl1(i+1)/wlmax_i>0.98.or.wlmax_i<0.01) then
        !  filtr(i) = 0
        ! else
        !  wlmax_i = (POR(i)+POR(i+1))*ROW0/((rosdry(i)+rosdry(i+1))*&
        !&  (1-(por(i)+por(i+1))/2)) - (wi1(i)+wi1(i+1))/2*row0_d_roi
        !  filtr(i) = (flwmax(i+1)+flwmax(i))/2*((wl1(i+1)+&
        !&  wl1(i))/2/wlmax_i)**(bh(i+1)+bh(i)+3)
        ! endif
        !enddo

        !SPLITTING-UP METHOD FOR TEMPERATURE EQUATION:
        !STEP 1: HEAT DIFFUSION

        !do i = 1, ns-1
        ! lammoist1 = 0.5*(lammoist(i)+lammoist(i+1))
        ! wsoil(i) = ( - lammoist1*(wl1(i+1)-wl1(i))/dzs(i) + &
        ! & filtr(i) ) *0.5*(rosdry(i)+rosdry(i+1))* &
        ! & (1-0.5*(por(i)+por(i+1)) )/row0
        !enddo


        !SPLITTING-UP METHOD FOR MOISTURE EQUATION:
        !STEP 1: MOISTURE DIFFUSION

        Wflow = 0.
        dhwfsoil=0.
        dhwfsoil1=0.
        dhwfsoil2=0.
        dhwfsoil3=0.
        dhwfsoil4=0.

        if ((l1 /= 0.and. h1 == 0.).or.ls1/=0.or. &
          & (hs1/=0.and.l1==0.and.h1==0)) then
          Wfhigh = 0
        !c(1)=1
        !b(1)=1
        !d(1)=0 
          c(1) = -1-dt*(lammoist(1)+lammoist(2))/(dzs(1)**2)
          b(1) =   -dt*(lammoist(1)+lammoist(2))/(dzs(1)**2)
          d(1) =   -wl1(1,is)
        endif
           
        if (h1 /= 0.and.ls1 == 0) then
          c(1)=1
          b(1)=0
          d(1)=WL_MAX(por(1),rosdry(1),wi1(1,is),tricemethhydr%par)
        !dhwfsoil1 = dt*(wl1(2)-wl1(1))*rosdry(1)*(1-por(1))/row0
        !& /dzs(1)*(lammoist(1)+lammoist(2))/2 !-
        !&  (WL_MAX(por(1),rosdry(1),wi1(1))-wl1(1)) VS,23.09.07
        !& *rosdry(1)*(1-por(1))*dzss(1)/row0  VS,23.09.07
        endif    
          
        lammoist1 = (lammoist(ns-1)+lammoist(ns))/2
        !c(ns)=-lammoist1/dzs(ns-1)
        !a(ns)=-lammoist1/dzs(ns-1)
        !d(ns)=Wflow
        c(ns) = -1-2*dt*lammoist1/dzs(ns-1)**2 
        a(ns) = -2*dt*lammoist1/dzs(ns-1)**2 
        d(ns) = -wl1(ns,is)

        do i=2,ns-1
          lammoist2 = (lammoist(i)+lammoist(i+1))/2
          lammoist1 = (lammoist(i-1)+lammoist(i))/2
          dzmean = (dzs(i-1)+dzs(i))/2
          a(i)=-dt/dzmean*lammoist1/dzs(i-1)
          b(i)=-dt/dzmean*lammoist2/dzs(i)
          c(i)=-dt/dzmean*(lammoist1/dzs(i-1) + &
          & lammoist2/dzs(i))-1
          d(i)=-wl1(i,is)
         !wl2(i) = wl1(i)+dt*(lammoist2/dzs(i)*(wl1(i+1)-wl1(i))-
        !& lammoist1/dzs(i-1)*(wl1(i)-wl1(i-1)))/dzmean
        enddo   
            
        call PROGONKA(a,b,c,d,wl2,1,ns)  
            
        if (h1 /= 0.and.ls1 == 0) then 
          lammoist1 = (lammoist(1)+lammoist(2))/2.
          dhwfsoil1 = -(wl2(1)*(dzs(1)/(2*dt)+lammoist1/dzs(1)) - &
          & wl2(2)*lammoist1/dzs(1)-wl1(1,is)*dzs(1)/(2*dt)) * &
          & rosdry(1)*(1-por(1))/row0*dt
        endif 
        !STEP 2 OF SPLITTING-UP METHOD FOR MOISTURE EQUATION: 
        !EVOLUTION OF MOISTURE DUE TO GRAVITATIONAL INFILTRATION  
         
        do i=2,ns-1
          wl3(i) = wl2(i) + dt*(filtr(i-1)-filtr(i))/dzss(i)
        enddo
         
        if ((h1==0.and.l1/=0).or.ls1 /= 0) then
          wl3(1) = wl2(1) + dt*(-filtr(1)+0)/dzss(1)
        endif

        if (h1/=0.and.ls1 == 0) then
          wl3(1) = WL_MAX(por(1),rosdry(1),wi1(1,is),tricemethhydr%par)
          dhwfsoil2 = - filtr(1)*rosdry(1)*(1-por(1))/row0*dt 
        endif

        filtr_low = 0
        wl3(ns) = wl2(ns) + dt*(-filtr_low+filtr(ns-1))/dzss(ns)

        do i=1,ns
          if (wl3(i)>WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)) then
            surplus=wl3(i)-WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)
            cy1:do j=1,ns
              if (WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par)-wl3(j)>0) then
                if (WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par) - wl3(j) > &
                & surplus*rosdry(i)*(1-por(i))*dzss(i) / &
                & (rosdry(j)*(1-por(j))*dzss(j))) then
                  wl3(j)=wl3(j)+surplus*rosdry(i)*(1-por(i))*dzss(i) / &
                  & (rosdry(j)*(1-por(j))*dzss(j))
                  surplus=0
                  exit cy1
                else
                  surplus = surplus - &
                  & (WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par)-wl3(j)) * &
                  & (rosdry(j)*(1-por(j))*dzss(j)) / &
                  & (rosdry(i)*(1-por(i))*dzss(i))
                  wl3(j)=WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par) 
                endif
              endif
            enddo cy1
            dhwfsoil3 = dhwfsoil3 + surplus*rosdry(i)*(1-por(i))*dzss(i)/row0
            wl3(i) = WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)
          endif
          if(wl3(i)<0) then
            if (i/=1) then 
              do j=i-1,1,-1
                if (wl3(j)>-wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
                & (rosdry(j)*(1-por(j))*dzss(j))) then
                  wl3(j)=wl3(j)+wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
                  & (rosdry(j)*(1-por(j))*dzss(j))
                  goto 1
                endif
              enddo
            endif
            if (i/=ns) then 
              do j=i+1,ns
                if (wl3(j)>-wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
                & (rosdry(j)*(1-por(j))*dzss(j))) then
                  wl3(j)=wl3(j)+wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
                  & (rosdry(j)*(1-por(j))*dzss(j))
                  goto 1
                endif
              enddo 
            endif
         !dhwfsoil = dhwfsoil - abs(wl4(i))*rosdry(i)*(1-por(i))*dzss(i)/row
        1   wl3(i)=0.
          endif
        enddo

        ! STEP 3 OF SPLITTING-UP METHOD : PHASE PROCESSES
        mhsep = (1. + tricemethhydr%par*alphamh)
        20 c2: do i = 1, ns
          work = MELTINGPOINT(Sals2(i,is)/wl3(i),pressoil(i),tricemethhydr%par)
          if (wi1(i,is) > 0 .and. Tsoil2(i) > work + 0.01) then
            Potphenergy = (Tsoil2(i) - (work + 0.01)) * &
            & csoil(i)*rosoil(i)*dzss(i)
            if (Potphenergy >= wi1(i,is)*rosdry(i)*dzss(i)*(1-por(i))*hicemelt) then
              wl4(i,is) = wl3(i) + wi1(i,is)/mhsep
              wi2(i,is) = 0.d0
              Tsoil3(i,is) = work + 0.01 + &
              & (Potphenergy-(wi1(i,is)*rosdry(i)*dzss(i)* &
              & (1 - por(i))*hicemelt))/(csoil(i)*rosoil(i)*dzss(i))
            else
              wl4(i,is) = wl3(i) + Potphenergy / &
              & (rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)/mhsep
              wi2(i,is) = wi1(i,is) - Potphenergy / &
              & (rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)
              Tsoil3(i,is) = work + 0.01
            endif
          else
            if (wl3(i) >= 0 .and. Tsoil2(i) < work - 0.01) then
              Potphenergy = - (Tsoil2(i) - work + 0.01) * &
              & csoil(i)*rosoil(i)*dzss(i)
              watavailtofreeze = (wl3(i) - UNFRWAT(Tsoil2(i),i)) * &
              & rosdry(i)*dzss(i)*(1 - por(i))*mhsep*hicemelt
              if (watavailtofreeze < 0) then 
                !cc = csoil(i)*rosoil(i)*dzss(i)
                !cc1 = rosdry(i)*dzss(i)*(1-por(i))*Lwi
                !bb = (-Potphenergy-0.1*cc-wl3(i)*cc1)/cc
                !aa = cc1/cc
                !call phase_iter(aa,bb,i,Tsoil3(i))
                Tsoil3(i,is) = Tsoil2(i) + watavailtofreeze/(csoil(i)*rosoil(i)*dzss(i))
                wi2(i,is) = wi1(i,is) + (wl3(i) - UNFRWAT(Tsoil2(i),i))*mhsep
                wl4(i,is) = UNFRWAT(Tsoil2(i),i)
                cycle c2  
              endif
              if (Potphenergy >= watavailtofreeze) then
                Tsoil3(i,is) = work - 0.01 - &
                & (Potphenergy - watavailtofreeze) / &
                & (csoil(i)*rosoil(i)*dzss(i))
                wi2(i,is) = wi1(i,is) + (wl3(i) - UNFRWAT(Tsoil2(i),i))*mhsep
                wl4(i,is) = UNFRWAT(Tsoil2(i),i)
              else
                wl4(i,is) = wl3(i) - Potphenergy/(rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)/mhsep
                wi2(i,is) = wi1(i,is) + Potphenergy/(rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)
                Tsoil3(i,is) = work - 0.01
              endif
            else
              wl4(i,is) = wl3(i)
              wi2(i,is) = wi1(i,is)
              Tsoil3(i,is) = Tsoil2(i)
            endif
          endif
        enddo c2

        max1 = 0.
        do i = 1, ns
          if (Tsoil3(i,is) < &
            & MELTINGPOINT(Sals2(i,is)/wl4(i,is),pressoil(i),tricemethhydr%par) - 0.01 &
            &.and. abs(wl4(i,is) - UNFRWAT(Tsoil3(i,is),i)) > max1) then
            max1 = abs(wl4(i,is) - UNFRWAT(Tsoil3(i,is),i))
          endif
        enddo    
  
if (max1 > 0.01 .and. iter3 < 20) then
  !wl3=(wl4+wl3)/2.
  !wi1=(wi2+wi1)/2.
  !Tsoil2=(Tsoil3+Tsoil2)/2.
  wl3 = wl3 + alp(iter+1)*(wl4(:,is)-wl3)
  wi1(:,is) = wi1(:,is) + &
  & alp(iter+1)*(wi2(:,is) - wi1(:,is))
  Tsoil2 = Tsoil2 + alp(iter+1)*(Tsoil3(:,is)-Tsoil2)
  iter = iter + 1
  iter3 = iter3 + 1
  if(iter == 8) iter = 0
  goto 20
else
  iter = 0
  iter3 = 0 
endif   

!evol=0

do i = 1, ns
  if (wi2(i,is) < 0) then
    wl4(i,is) = wl4(i,is) + wi2(i,is)/mhsep
    wi2(i,is) = 0
  endif
enddo

! Methane generation due to methane hydrate dissociation
do i = 1, ns
  methgenmh(i) = - tricemethhydr%par*(wi2(i,is) - wi1(i,is))*1.d+3 / &
  & (dt*mhsep*molmass_h2o)*nch4_d_nh2o*rosdry(i)*(1. - por(i))
enddo

do i=1,ns
  if (wl4(i,is) > &
  & POR(i)*ROW0/(rosdry(i)*(1-por(i))) - wi2(i,is)*row0_d_roi) then
    surplus = wl4(i,is) - &
    & (POR(i)*ROW0/(rosdry(i)*(1-por(i))) - wi2(i,is)*row0_d_roi)
    c3:do j=1,ns
      if (POR(j)*ROW0/(rosdry(j)*(1-por(j))) - wi2(j,is)*row0_d_roi - &
      & wl4(j,is) > 0) then
        if (POR(j)*ROW0/(rosdry(j)*(1-por(j))) - wi2(j,is)*row0_d_roi - &
        & wl4(j,is) > surplus*rosdry(i)*(1 - por(i))*dzss(i) / &
        & (rosdry(j)*(1-por(j))*dzss(j))) then
          wl4(j,is)=wl4(j,is)+surplus*rosdry(i)*(1-por(i))*dzss(i) / &
          & (rosdry(j)*(1-por(j))*dzss(j))
          surplus=0
          exit c3
        else
          surplus=surplus-(POR(j)*ROW0/(rosdry(j)*(1-por(j))) - &
          & wi2(j,is)*row0_d_roi-wl4(j,is)) * &
          & (rosdry(j)*(1-por(j))*dzss(j)) / &
          & (rosdry(i)*(1-por(i))*dzss(i))
          wl4(j,is) = &
          & POR(j)*ROW0/(rosdry(j)*(1-por(j))) - wi2(j,is)*row0_d_roi 
        endif
      endif
    enddo c3
    dhwfsoil4 = dhwfsoil4 + surplus*rosdry(i)*(1-por(i))*dzss(i)/row0
    wl4(i,is) = &
    & POR(i)*ROW0/(rosdry(i)*(1-por(i))) - wi2(i,is)*row0_d_roi 
  endif
  if(wl4(i,is) < 0) then
    if (i/=1) then 
      do j=i-1,1,-1
        if (wl4(j,is) > - wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
        & (rosdry(j)*(1-por(j))*dzss(j))) then
          wl4(j,is) = wl4(j,is) + &
          & wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
          & (rosdry(j)*(1-por(j))*dzss(j))
          goto 2
        endif
      enddo
    endif
    if (i/=ns) then 
      do j=i+1,ns
        if (wl4(j,is) > -wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
        & (rosdry(j)*(1-por(j))*dzss(j))) then
          wl4(j,is) = wl4(j,is) + &
          & wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
          & (rosdry(j)*(1-por(j))*dzss(j))
          goto 2
        endif
      enddo
    endif
    !dhwfsoil = dhwfsoil - abs(wl4(i))*rosdry(i)*(1-por(i))*dzss(i)/row
2   wl4(i,is) = 0
  endif
enddo

  !!!!! II. FREEZING AND MELTING IN CASE OF TSOIL<MELTING POINT (O C)!!!!!

    !do i=1, ns
    ! if (Tsoil3(i)<-1.) then
    !   wi3(i) = wi2(i) + (wl4(i)-unfrwat(Tsoil3(i),i))
  !   wl5(i) = unfrwat(Tsoil3(i),i)
    !   Tsoil4(i) = Tsoil3(i) + & 
    !   (wl4(i)-unfrwat(Tsoil3(i),i))*rosdry(i)*dzss(i)*(1-por(i))*Lwi/&
    !   (csoil(i)*rosoil(i)*dzss(i))
    ! else
    !   wi3(i) = wi2(i)
    !   wl5(i) = wl4(i)
    !   Tsoil4(i) = Tsoil3(i)
  ! end if
    ! if (wl5(i)>POR(i)*ROW/(rosdry(i)*(1-por(i))) - wi3(i)*row/roi) then
    !  surplus = wl5(i)-(POR(i)*ROW/(rosdry(i)*(1-por(i))) - wi3(i)*row/roi)
    !  dhwfsoil = dhwfsoil + surplus*rosdry(i)*(1-por(i))*dzss(i)/row
    !  wl5(i) = POR(i)*ROW/(rosdry(i)*(1-por(i))) - wi3(i)*row/roi 
  ! endif
    ! if(wl5(i)<0) then 
    !  dhwfsoil = dhwfsoil - abs(wl5(i))*rosdry(i)*(1-por(i))*dzss(i)/row
    !  wl5(i)=0
    ! endif
  !enddo
   
   ! evol=0
!   evol1=0

!   do i=1,ns
!    evol=evol+(Tsoil3(i)-Tsoil1(i))*rosoil(i)*csoil(i)*dzss(i)
!    evol1=evol1+(wi2(i)-wi1(i))*rosdry(i)*(1-por(i))*dzss(i)*Lwi
!   enddo

dhwfsoil = dhwfsoil1 + dhwfsoil2 + dhwfsoil3 + dhwfsoil4
  
Tsoil1(:,is) = Tsoil3(:,is)
wl1   (:,is) = wl4   (:,is)
wi1   (:,is) = wi2   (:,is)

! Diagnostic calculation of air content, kg/kg 
! (air mass per unit mass of the dry soil)
! assuming that the air occupies all the pore space 
! free from liquid water and ice
do i = 1, ns
  wa(i) = por(i)*roa0/((1.-por(i))*rosdry(i)) - &
  & roa0_d_row0*wl1(i,is) - roa0_d_roi*wi1(i,is)
  wa(i) = max(wa(i), 0.d0)
enddo

deallocate(pressoil)

if (firstcall) firstcall=.false.
RETURN
END SUBROUTINE SOILFORLAKE


SUBROUTINE SOIL_COND_HEAT_COEF(is)

use ARRAYS_SOIL, only : rosdry,csoil,rosoil,lamsoil, &
lammoist,filtr,wsoil,por,bh,wl1,wi1,psimax,flwmax,dlmax,dzs
use ARRAYS_GRID, only : nsoilcols
use DRIVING_PARAMS, only : ns, tricemethhydr
use PHYS_CONSTANTS, only : &
& row0, &
& cw, &
& ci, &
& roi, &
& row0_d_roi
use PHYS_FUNC, only : &
& WL_MAX, &
& SOIL_COND_JOHANSEN
use METH_OXYG_CONSTANTS, only : &
& cpmethhydr


implicit none

! Input variables
integer(kind=iintegers), intent(in) :: is

real(kind=ireals), parameter :: cr = 0.2*4180.d0

real(kind=ireals) :: wlmax_i
real(kind=ireals) :: bb1
real(kind=ireals) :: arg
real(kind=ireals) :: psi
real(kind=ireals) :: pf
real(kind=ireals) :: lammoist1
real(kind=ireals), save, pointer :: cimh

integer(kind=iintegers) :: i
logical, save :: firstcall = .true.


if (firstcall) then
  if (tricemethhydr%par > 0.) then
    cimh => cpmethhydr
  else
    cimh => ci
  endif
endif

!rosdry(:) = 1200. !2400. 
   
do i = 1, ns
  wlmax_i = POR(i)*ROW0/(rosdry(i)*(1-por(i))) 
  BB1 = BH(i) + 2.d0
  csoil(i) = (cr + wl1(i,is)*CW + wi1(i,is)*cimh) !*1.d-3 !1.d-3 - for sensitivity experiment
  rosoil(i) = rosdry(i)*(1 - por(i))*(1 + wi1(i,is) + wl1(i,is)) 
!  ARG = (WL1(i)+WI1(i))/wlmax_i
!  ARG = max(ARG,1.d-2)
!  PSI = PSIMAX(i)*(ARG)**(-BH(i))
!  PF = log(-PSI*100.)/log(10.d0) 
!  if(PF>5.1) then
!    lamsoil(i) = 4.1E-4*418.
!  else
!    lamsoil(i) = exp(-PF-2.7)*418.
!  end if
  lamsoil(i) = SOIL_COND_JOHANSEN(wl1(i,is),wi1(i,is),rosdry(i),por(i))
  wlmax_i = WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)
  ARG = (wl1(i,is))/wlmax_i
  ARG = max(ARG,1.d-2)
  lammoist(i) = dlmax(i)*ARG**BB1
end do

do i = 1, ns-1
  wlmax_i = WL_MAX(por(i+1),rosdry(i+1),wi1(i+1,is),tricemethhydr%par)
  if (wlmax_i < 0.01) then !wl1(i+1,is)/wlmax_i>0.98.or.
    filtr(i) = 0
  else
    wlmax_i = WL_MAX(0.5*(por(i)+por(i+1)),0.5*(rosdry(i)+rosdry(i+1)), &
    & 0.5*(wi1(i,is) + wi1(i+1,is)),tricemethhydr%par)
    filtr(i) = 0.5*(flwmax(i+1)+flwmax(i))*(0.5*(wl1(i+1,is) + &
    & wl1(i,is))/wlmax_i)**(bh(i+1)+bh(i)+3)
  endif
enddo

do i = 1, ns-1
  lammoist1 = 0.5*(lammoist(i)+lammoist(i+1))
  wsoil(i) = ( - lammoist1*(wl1(i+1,is) - wl1(i,is))/dzs(i) + &
  & filtr(i) ) *0.5*(rosdry(i)+rosdry(i+1))* &
  & (1. - 0.5*(por(i)+por(i+1)) )/row0
enddo 

if (firstcall) firstcall = .false.
END SUBROUTINE SOIL_COND_HEAT_COEF


SUBROUTINE SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm, &
                       & ddz,ddzi,zsoilcols, &
                       & wst,SR,a,b,c,d,add_to_winter, &
                       & bathymwater,bathymice, &
                       & bathymdice,bathymsoil, &
                       & soilflux,fdiffbot, contr) 

! Subroutine calculates temperature in soil columns,
! excepting for the lowest column, assuming 
! nsoilcols = M+1

use LAKE_DATATYPES, only : &
& iintegers, ireals

use DRIVING_PARAMS, only : soilcolconjtype

use T_SOLVER_MOD, only : &
& T_DIFF

use ATMOS, only : &
& pressure, wind10

use PHYS_CONSTANTS, only : &
& cw_m_row0, z0_bot

use PHYS_FUNC, only : &
& LOGFLUX, TEMPWATR

use ARRAYS_BATHYM, only : bathym, dhw, dhw0, layers_type
use ARRAYS_SOIL, only : csoil, rosoil, lamsoil, Tsoil3, Tsoil2, Tsoil1, &
& wl4,wi2,wa,Sals2,rosdry,por,lsm,dzs,zsoil,rosoil
use ARRAYS_GRID, only : ddz05, gridsize_type, gridspacing_type
use ARRAYS_METHANE, only : veg,TgrAnn,methgenmh,qwater,lammeth, &
& fplant,febul0, fdiff_lake_surf, &
& plant_sum,bull_sum,oxid_sum,rprod_sum, anox, &
& rprod_total_oldC,rprod_total_newC, &
& ice_meth_oxid_total, &
& h_talik,tot_ice_meth_bubbles,rootss 
use ARRAYS_WATERSTATE, only : Tw2, waterstate_type
use ARRAYS_TURB, only: eps1
use ARRAYS_OXYGEN, only : sodbot
use ARRAYS, only : gas, dt_inv

use METHANE_MOD, only : METHANE

use METH_OXYG_CONSTANTS, only : CH4_exp_growth

implicit none

! Input variables

type(gridsize_type), intent(inout) :: gs
type(gridspacing_type), intent(in) :: gsp
type(layers_type), intent(in) :: ls

real(kind=ireals),      intent(in) :: dt   ! timestep, st
real(kind=ireals),      intent(inout) :: ftot
real(kind=ireals),      intent(in) :: ch4_pres_atm
real(kind=ireals),      intent(in) :: ddz(1:gs%M), ddzi(1:gs%Mice), zsoilcols(1:gs%nsoilcols+1) 
type(waterstate_type),  intent(in) :: wst
real(kind=ireals),      intent(in) :: SR(0:gs%M+1) ! Shortwave radiation
type(bathym),           intent(in) :: bathymwater(1:gs%M+1), bathymice(1:gs%Mice+1) 
type(bathym),           intent(in) :: bathymdice(1:gs%Mice+1), bathymsoil(1:gs%nsoilcols+1)

integer(kind=iintegers), intent(in) :: contr(1:2)

real(kind=ireals), intent(inout), dimension(1:vector_length) :: a, b, c, d

logical, intent(inout) :: add_to_winter

! Output variables
real(kind=ireals),      intent(inout) :: soilflux(1:gs%nsoilcols), fdiffbot(1:gs%nsoilcols)

! Local variables
integer(kind=iintegers), parameter :: bcswitch = 1 !1 - continuity of flux and temperature
                                                   !2 - continuity of flux, calculated from surface layer theory

integer(kind=iintegers) :: is, i ! Number of soil column
real(kind=ireals) :: Tsoilsurf(1:gs%nsoilcols), qsoilsurf(1:gs%nsoilcols)
real(kind=ireals) :: exchcoef, Tws, qwaters
real(kind=ireals), allocatable :: work(:), work1(:)
integer(kind=iintegers), allocatable :: iwork(:)
real(kind=ireals) :: xx

if ((ls%l1 /= 0. .or. ls%ls1 /= 0) .and. bcswitch == 2) then
  print*, "bcswitch == 2 is not operational when &
  &l1 /= 0. .or. ls1 /= 0", ls%l1, ls%ls1
  STOP
endif


if (bcswitch == 2) then
  allocate(work1(1:gs%nsoilcols)) 
  allocate(iwork(1:gs%nsoilcols))
  if     (soilcolconjtype == 1) then
    forall(is=1:gs%nsoilcols-1) work1(is) = bathymsoil(is)%dzSLc
    forall(is=1:gs%nsoilcols-1) iwork(is) = bathymsoil(is)%icent
  elseif (soilcolconjtype == 2) then
    forall(is=1:gs%nsoilcols-1) work1(is) = bathymsoil(is)%dzSL
    forall(is=1:gs%nsoilcols-1) iwork(is) = bathymsoil(is)%itop
  endif
endif


if (contr(1) == 1) then
  ! Heat equation in soil columns
  if (bcswitch == 1) then ! Continuity of flux and temperature across soil-water interface
    if (soilcolconjtype == 1) then
      do is = 1, gs%nsoilcols-1
        i = bathymsoil(is)%icent
        Tsoilsurf(is) = wst%Tw1(i)
      enddo
    elseif (soilcolconjtype == 2) then
      call TSURFSOILCOL(gs%M,gs%Mice,gs%nsoilcols,ls%h1,ls%l1,ls%ls1, &
                       & ddz,ddzi,zsoilcols, &
                       & wst%Tw2,wst%Ti2,wst%Tis2, &
                       & bathymwater,bathymice,bathymdice,bathymsoil,&
                       & Tsoilsurf)
    endif
    
    do is = 1, gs%nsoilcols-1
      call SOIL_COND_HEAT_COEF(is)
      call T_DIFF(1,Tsoilsurf(is),dt,0,0,0,0,is,.false.)
      soilflux(is) = &
      & csoil(1)*rosoil(1)*dt_inv*0.5*dzs(1)*( Tsoil2(1) - Tsoil1(1,is) ) - &
      & 0.25* ( lamsoil(1) + lamsoil(2) ) * &
      & ( Tsoil2(2) + Tsoil1(2,is) - Tsoil2(1) - Tsoil1(1,is) )/dzs(1)
      call SOILFORLAKE(dt,a,b,c,d,is)
    enddo
  elseif (bcswitch == 2) then ! continuity of flux, calculated from surface layer theory above soil surface
    do is = 1, gs%nsoilcols-1
      call SOIL_COND_HEAT_COEF(is)
      i = iwork(is)
      soilflux(is) = LOGFLUX(sqrt(wst%u1(i)*wst%u1(i) + wst%v1(i)*wst%v1(i)), &
      & Tsoil1(1,is) - wst%Tw1(i), work1(is), z0_bot, z0_bot, cw_m_row0, 1) + SR(i)
      call T_DIFF(1,soilflux(is),dt,0,0,0,0,is,.true.)
      call SOILFORLAKE(dt,a,b,c,d,is)
    enddo
  endif
endif

if (contr(2) == 1) then
  ! Methane equation in soil columns
  if (bcswitch == 1) then ! Continuity of flux and concentration across soil-water interface
    if     (soilcolconjtype == 1) then
      do is = 1, gs%nsoilcols-1
        i = bathymsoil(is)%icent
        qsoilsurf(is) = qwater(i,1)
      enddo
    elseif (soilcolconjtype == 2) then
      allocate (work(1:gs%Mice+1)); work = 0.
      call TSURFSOILCOL(gs%M,gs%Mice,gs%nsoilcols,ls%h1,ls%l1,ls%ls1, &
                       & ddz,ddzi,zsoilcols, &
                       & qwater(1,2),work,work, &
                       & bathymwater,bathymice,bathymdice,bathymsoil,&
                       & qsoilsurf)
      deallocate (work)
    endif
    
    allocate (work(0:gs%M+1))
    do is = 1, gs%nsoilcols-1
      gs%isoilcol = is
      call METHANE &
      & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,is),rosoil, &
      & work, work, &
      & wl4(1,is),wi2(1,is),wa,Sals2, &
      & rootss,rosdry,por,veg,qsoilsurf(is)*por(1),TgrAnn, methgenmh, &
      & ddz,ddz05, wst, lammeth, 0.5*(zsoilcols(is) + zsoilcols(is+1)), &
      & ls, dhw, dhw0, .false., .false., lsm, bathymwater, &
      & fplant, febul0(is), fdiffbot(is), ftot, fdiff_lake_surf, &
      & plant_sum,bull_sum,oxid_sum,rprod_sum, &
      & anox,gs,gsp,dt,eps1(1),sodbot, &
      & rprod_total_oldC, rprod_total_newC, ice_meth_oxid_total, &
      & h_talik,tot_ice_meth_bubbles,add_to_winter) 
    enddo
    deallocate(work)
  elseif (bcswitch == 2) then! continuity of flux, calculated from surface layer theory above soil surface
    allocate (work(0:gs%M+1))
    do is = 1, gs%nsoilcols-1
      i = iwork(is)
      ! Exchange coefficient
      !exchcoef = LOGFLUX(sqrt(wst%u1(i)*wst%u1(i) + wst%v1(i)*wst%v1(i)), &
      !& qsoil(1,is)/por(1) - qwater(i,1), bathymsoil(is)%dzSL, z0_bot, z0_bot, 1._ireals,2)
      ! Methane concentration in water from similarity profile
      !qwaters = TEMPWATR(lammeth(i),bathymwater(i)%rad_int,qsoil(1,is)/por(1),exchcoef, &
      !& qwater(i,1), lammeth(i)*(qwater(i+1,1) - qwater(i,1))/(ddz(i)*h1))
      ! Methane flux at the bottom
      ! Taking into account depletion of methane in top sediments
      if (CH4_exp_growth /= 0.) then
        xx = CH4_exp_growth*0.5*dzs(1) / ( (exp(CH4_exp_growth*0.5*dzs(1)) - 1.) * por(1) )
      else
        xx = 1./por(1)
      endif
      fdiffbot(is) = LOGFLUX(sqrt(wst%u1(i)*wst%u1(i) + wst%v1(i)*wst%v1(i)), &
      & gas%qsoil(1,is)*xx - qwater(i,1), work1(is), &
      & z0_bot, z0_bot, 1._ireals,1)

      gs%isoilcol = is
      call METHANE &
      & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,is),rosoil, &
      & work, work, &
      & wl4(1,is),wi2(1,is),wa,Sals2, &
      & rootss,rosdry,por,veg,fdiffbot(is),TgrAnn, methgenmh, &
      & ddz,ddz05, wst, lammeth, 0.5*(zsoilcols(is) + zsoilcols(is+1)), &
      & ls, dhw, dhw0, .false., .true., lsm, bathymwater, &
      & fplant, febul0(is), fdiffbot(is), ftot, fdiff_lake_surf, &
      & plant_sum,bull_sum,oxid_sum,rprod_sum, &
      & anox,gs,gsp,dt,eps1(1),sodbot, &
      & rprod_total_oldC, rprod_total_newC, ice_meth_oxid_total, &
      & h_talik,tot_ice_meth_bubbles,add_to_winter) 
    enddo
    deallocate(work)
  endif

endif

if (bcswitch == 2) then
  deallocate(work1) 
  deallocate(iwork)
endif

END SUBROUTINE SOILCOLSTEMP


SUBROUTINE TSURFSOILCOL(M,Mice,nsoilcols,h1,l1,ls1,ddz,ddzi,zsoilcols, &
                       & Tw2,Ti2,Tis2, &
                       & bathymwater,bathymice,bathymdice,bathymsoil, &
                       & Tsoilsurf)

! Subroutine TSURFSOILCOL calculates mean temperature
! of a lake cross-section over specified depth intervals


use NUMERIC_PARAMS!, only : &
!& small_value, pi

use ARRAYS_BATHYM, bathymwater_=>bathymwater,bathymice_=>bathymice, &
& bathymsoil_=>bathymsoil,bathymdice_=>bathymdice,&
& h1_=>h1,l1_=>l1,ls1_=>ls1!, only : &
!& bathym

implicit none

!Input variables
integer(kind=iintegers), intent(in) :: M, Mice, nsoilcols
real(kind=ireals),    intent(in) :: h1, l1, ls1
real(kind=ireals),    intent(in) :: ddz(1:M), ddzi(1:Mice)
real(kind=ireals),    intent(in) :: zsoilcols(1:nsoilcols+1)
real(kind=ireals),    intent(in) :: Tw2(1:M+1)
real(kind=ireals),    intent(in) :: Ti2(1:Mice+1), Tis2(1:Mice+1)
type(bathym), intent(in) :: bathymwater(1:M+1)  , bathymice(1:Mice+1) 
type(bathym), intent(in) :: bathymdice(1:Mice+1), bathymsoil(1:nsoilcols+1)

!Output variables
real(kind=ireals), intent(out) :: Tsoilsurf(1:nsoilcols)

! Local variables
integer(kind=iintegers) :: i, j, j0, j1, nl 

real(kind=ireals), allocatable :: Temp(:), z(:), per(:), dlx(:), dly(:), aell(:), bell(:)
real(kind=ireals) :: pers,dz,aells,bells,areas,area,dlxs,dlys
!type(bathym),      allocatable :: bathymall(:)

logical, save :: firstcall = .true.


! Water, ice and deep ice
if (l1  > small_value .and. &
  & h1  > small_value .and. &
  & ls1 > small_value) then

  nl = 2*Mice + M
  allocate (Temp(1:nl+1), z(0:nl+1), per(0:nl+1), dlx(0:nl+1), dly(0:nl+1))
  allocate (aell(0:nl+1), bell(0:nl+1))
!  allocate (bathymall(1:nl+1))

  ! Ice top
  z(0)    = 0.                  ; z(1)   = 0.5*ddzi(1)*l1
  aell(0)  = 0.5*bathymice(1)%Lx ; aell(1) = 0.5*bathymice(1)%Lx_half
  bell(0)  = 0.5*bathymice(1)%Ly ; bell(1) = 0.5*bathymice(1)%Ly_half
  Temp(1) = Ti2(1)
  j = 2

  ! Ice interior
  do i = 2, Mice
    Temp(j) = Ti2(i)
    aell(j) = 0.5*bathymice(i)%Lx_half
    bell(j) = 0.5*bathymice(i)%Ly_half   
    z(j)  = z(j-1) + 0.5*(ddzi(i-1) + ddzi(i))*l1
    j = j + 1
  enddo

  ! Ice-water interface
  Temp(j) = Ti2(Mice+1)
  z(j) = z(j-1) + 0.5*(ddzi(Mice)*l1 + ddz(1)*h1)
  aell(j) = 0.5*bathymwater(1)%Lx_half
  bell(j) = 0.5*bathymwater(1)%Ly_half
  j = j + 1

  ! Water interior
  do i = 2, M
    Temp(j) = Tw2(i)
    aell(j) = 0.5*bathymwater(i)%Lx_half
    bell(j) = 0.5*bathymwater(i)%Ly_half
    z(j) = z(j-1) + 0.5*(ddz(i-1) + ddz(i))*h1
    j = j + 1
  enddo

  ! Water-deepice interface
  Temp(j) = Tw2(M+1)
  z(j) = z(j-1) + 0.5*(ddz(M)*h1 + ddzi(1)*ls1)
  aell(j) = 0.5*bathymdice(1)%Lx_half
  bell(j) = 0.5*bathymdice(1)%Ly_half
  j = j + 1

  ! Deepice interior
  do i = 2, Mice
    Temp(j) = Tis2(i)
    aell(j) = 0.5*bathymdice(i)%Lx_half
    bell(j) = 0.5*bathymdice(i)%Ly_half
    z(j) = z(j-1) + 0.5*(ddzi(i-1) + ddzi(i))*ls1
    j = j + 1
  enddo 

  ! Deepice bottom
  Temp(j) = Tis2(Mice+1)
  z(j) = z(j-1) + 0.5*ddzi(Mice)*ls1
  aell(j) = 0.5*bathymdice(Mice+1)%Lx
  bell(j) = 0.5*bathymdice(Mice+1)%Ly

  z(:) = z(:) + zsoilcols(nsoilcols+1) - h1 - l1 - ls1
endif

! Water and ice
if (l1  > small_value .and. &
  & h1  > small_value .and. &
  & ls1 < small_value) then

  nl = Mice + M
  allocate (Temp(1:nl+1), z(0:nl+1), per(0:nl+1), dlx(0:nl+1), dly(0:nl+1))
  allocate (aell(0:nl+1), bell(0:nl+1))
!  allocate (bathymall(1:nl+1))

  ! Ice top
  z(0)    = 0.                  ; z(1)   = 0.5*ddzi(1)*l1
  aell(0)  = 0.5*bathymice(1)%Lx ; aell(1) = 0.5*bathymice(1)%Lx_half
  bell(0)  = 0.5*bathymice(1)%Ly ; bell(1) = 0.5*bathymice(1)%Ly_half
  Temp(1) = Ti2(1)
  j = 2

  do i = 2, Mice
    Temp(j) = Ti2(i)
    aell(j) = 0.5*bathymice(i)%Lx_half
    bell(j) = 0.5*bathymice(i)%Ly_half   
    z(j)  = z(j-1) + 0.5*(ddzi(i-1) + ddzi(i))*l1
    j = j + 1
  enddo

  ! Ice-water interface
  Temp(j) = Ti2(Mice+1)
  z(j) = z(j-1) + 0.5*(ddzi(Mice)*l1 + ddz(1)*h1)
  aell(j) = 0.5*bathymwater(1)%Lx_half
  bell(j) = 0.5*bathymwater(1)%Ly_half
  j = j + 1

  ! Water interior
  do i = 2, M
    Temp(j) = Tw2(i)
    aell(j) = 0.5*bathymwater(i)%Lx_half
    bell(j) = 0.5*bathymwater(i)%Ly_half
    z(j) = z(j-1) + 0.5*(ddz(i-1) + ddz(i))*h1
    j = j + 1
  enddo

  ! Water bottom
  Temp(j) = Tw2(M+1)
  z(j) = z(j-1) + 0.5*ddz(M)*h1
  aell(j) = 0.5*bathymwater(M+1)%Lx
  bell(j) = 0.5*bathymwater(M+1)%Ly

  z(:) = z(:) + zsoilcols(nsoilcols+1) - h1 - l1
endif

! Water and deep ice
if (l1  < small_value .and. &
  & h1  > small_value .and. &
  & ls1 > small_value) then

  nl = Mice + M
  allocate (Temp(1:nl+1), z(0:nl+1), per(0:nl+1), dlx(0:nl+1), dly(0:nl+1))
  allocate (aell(0:nl+1), bell(0:nl+1))
!  allocate (bathymall(1:nl+1))

  ! Water top
  z(0)    = 0.                    ; z(1)   = 0.5*ddz(1)*h1
  aell(0)  = 0.5*bathymwater(1)%Lx ; aell(1) = 0.5*bathymwater(1)%Lx_half
  bell(0)  = 0.5*bathymwater(1)%Ly ; bell(1) = 0.5*bathymwater(1)%Ly_half
  Temp(1) = Tw2(1)
  j = 2

  ! Water interior
  do i = 2, M
    Temp(j) = Tw2(i)
    aell(j) = 0.5*bathymwater(i)%Lx_half
    bell(j) = 0.5*bathymwater(i)%Ly_half
    z(j) = z(j-1) + 0.5*(ddz(i-1) + ddz(i))*h1
    j = j + 1
  enddo

  ! Water-deepice interface
  Temp(j) = Tw2(M+1)
  z(j) = z(j-1) + 0.5*(ddz(M)*h1 + ddzi(1)*ls1)
  aell(j) = 0.5*bathymdice(1)%Lx_half
  bell(j) = 0.5*bathymdice(1)%Ly_half
  j = j + 1

  ! Deep ice interior
  do i = 2, Mice
    Temp(j) = Tis2(i)
    aell(j) = 0.5*bathymdice(i)%Lx_half
    bell(j) = 0.5*bathymdice(i)%Ly_half
    z(j) = z(j-1) + 0.5*(ddzi(i-1) + ddzi(i))*ls1
    j = j + 1
  enddo

  ! Deepice bottom
  Temp(j) = Tis2(Mice+1)
  z(j) = z(j-1) + 0.5*ddzi(Mice)*ls1
  aell(j) = 0.5*bathymdice(Mice+1)%Lx
  bell(j) = 0.5*bathymdice(Mice+1)%Ly

  z(:) = z(:) + zsoilcols(nsoilcols+1) - h1 - ls1
endif

! Water only
if (l1  < small_value .and. &
  & h1  > small_value .and. &
  & ls1 < small_value) then

  nl = M
  allocate (Temp(1:nl+1), z(0:nl+1), per(0:nl+1), dlx(0:nl+1), dly(0:nl+1))
  allocate (aell(0:nl+1), bell(0:nl+1))
!  allocate (bathymall(1:nl+1))

  ! Water top
  z(0)    = 0.                    ; z(1)   = 0.5*ddz(1)*h1
  aell(0)  = 0.5*bathymwater(1)%Lx ; aell(1) = 0.5*bathymwater(1)%Lx_half
  bell(0)  = 0.5*bathymwater(1)%Ly ; bell(1) = 0.5*bathymwater(1)%Ly_half
  Temp(1) = Tw2(1)
  j = 2

  do i = 2, M
    Temp(j) = Tw2(i)
    aell(j) = 0.5*bathymwater(i)%Lx_half
    bell(j) = 0.5*bathymwater(i)%Ly_half
    z(j) = z(j-1) + 0.5*(ddz(i-1) + ddz(i))*h1
    j = j + 1
  enddo

  ! Water bottom
  Temp(j) = Tw2(M+1)
  z(j) = z(j-1) + 0.5*ddz(M)*h1
  aell(j) = 0.5*bathymwater(M+1)%Lx
  bell(j) = 0.5*bathymwater(M+1)%Ly

  z(:) = z(:) + zsoilcols(nsoilcols+1) - h1
endif

!! Cross-section parameters calculation, assuming elliptic shape
!do i = 0, nl+1
!  per(i) = ELLAR(aell(i),bell(i))
!enddo
!! Assuming cross-section area reduces with depth
!do i = 0, nl
!  dlx(i) = aell(i) - aell(i+1) 
!  dly(i) = bell(i) - bell(i+1)
!  dlx(i) = sqrt(( z(i+1) - z(i) ) * ( z(i+1) - z(i) ) + dlx(i)*dlx(i))
!  dly(i) = sqrt(( z(i+1) - z(i) ) * ( z(i+1) - z(i) ) + dly(i)*dly(i))
!enddo

do i = 1, nsoilcols-1 ! Except for the lowest soil column
  ! Seeking the gridcells falling into the depth interval of current soil column
  if (zsoilcols(i) < z(0)) then
    j0 = 1
  else
    do j = 0, nl
      if (z(j) <= zsoilcols(i) .and. z(j+1) > zsoilcols(i)) then
        j0 = j+1
        exit
      endif
    enddo 
  endif
  do j = j0, nl
    if (z(j) <= zsoilcols(i+1) .and. z(j+1) > zsoilcols(i+1)) then
      j1 = j+1
      exit
    endif
  enddo

  ! Integrating temperature over the range of depths of current soil column,
  ! assuming at least two numerical layers overlay with this range.
  ! The bottom temperature is areally-averaged.
  areas = 0.
  Tsoilsurf(i) = 0.

  ! Top numerical layer, overlapping with the soil column depth interval
  if (zsoilcols(i) > z(j0-1)) then
    ! Soil column top is inside the numerical layer
    !dz = z(j0) - zsoilcols(i)
    aells = 0.5*bathymsoil(i)%Lx; bells = 0.5*bathymsoil(i)%Ly
    !pers  = ELLAR(aells,bells)
    !dlxs  = sqrt((aells - aell(j0))*(aells - aell(j0)) + dz*dz)
    !dlys  = sqrt((bells - bell(j0))*(bells - bell(j0)) + dz*dz)
    !area  = 0.25*(per(j0) + pers)*(dlxs + dlys)
    area  = pi*(aells*bells - aell(j0)*bell(j0))
    Tsoilsurf(i) = Tsoilsurf(i) + Temp(j0)*area !(zsoilcols(i+1) - z(j1-1))
    areas = areas + area
  else
    ! Soil column top is outside the numerical layer (j0 = 1)
    ! area = 0.25*(per(j0-1) + per(j0))*(dlx(j0-1) + dly(j0-1))
    area  = pi*(aell(j0-1)*bell(j0-1) - aell(j0)*bell(j0))
    Tsoilsurf(i) = Tsoilsurf(i) + Temp(j0)*area  !(z(j) - z(j-1))
    areas = areas + area
  endif

  ! Interior numerical layers
  if (j0+1 <= j1-1) then
    do j = j0+1, j1-1
      !area = 0.25*(per(j-1) + per(j))*(dlx(j-1) + dly(j-1))
      area  = pi*(aell(j-1)*bell(j-1) - aell(j)*bell(j))
      Tsoilsurf(i) = Tsoilsurf(i) + Temp(j)*area  !(z(j) - z(j-1))
      areas = areas + area
    enddo
  endif

 ! Lowest numerical layer, overlapping with the soil column depth interval
  !dz = zsoilcols(i+1) - z(j1-1)
  aells = 0.5*bathymsoil(i+1)%Lx; bells = 0.5*bathymsoil(i+1)%Ly
  !pers  = ELLAR(aells,bells)
  !dlxs  = sqrt((aell(j1-1) - aells)*(aell(j1-1) - aells) + dz*dz)
  !dlys  = sqrt((bell(j1-1) - bells)*(bell(j1-1) - bells) + dz*dz)
  !area  = 0.25*(per(j1-1) + pers)*(dlxs + dlys)
  area  = pi*(aell(j1-1)*bell(j1-1) - aells*bells)
  Tsoilsurf(i) = Tsoilsurf(i) + Temp(j1)*area !(zsoilcols(i+1) - z(j1-1))
  areas = areas + area

  ! Averaging
  Tsoilsurf(i) = Tsoilsurf(i)/areas

enddo

deallocate(Temp, z, per, dlx, dly, aell, bell)
!deallocate(bathymall)

if (firstcall) firstcall = .false.


contains
FUNCTION ELLAR(a,b)
! Approximate formula for the ellipse perimeter
implicit none

real(kind=ireals) :: ELLAR

real(kind=ireals), intent(in) :: a, b

ELLAR = pi*(3.*(a + b) - sqrt( (3.*a + b)*(a + 3.*b) ) )

END FUNCTION ELLAR


END SUBROUTINE TSURFSOILCOL


SUBROUTINE LATERHEAT(ix,iy,gs,ls,btmw,btmi,btmdi,btms,gsp,soilflux,onlywater,lsh)

!Subroutine LATERHEAT calculates heat sources from lateral heat exchange
! with soil columns for water, ice and deep ice layers

use ARRAYS_BATHYM,ls_=>ls!, only : bathym, layers_type
use ARRAYS_GRID,gs_=>gs,gsp_=>gsp!, only : gridsize_type, gridspacing_type
use ARRAYS_SOIL,lsh_=>lsh,soilflux_=>soilflux!, only : lsh_type

use NUMERIC_PARAMS!, only : &
!& small_value

implicit none

! Input variables
integer(kind=iintegers), intent(in) :: ix, iy
type(gridsize_type),     intent(in) :: gs ! Gridsizes
type(layers_type),       intent(in) :: ls ! Layers thicknesses
type(bathym),            intent(in) :: btmw(1:gs%M+1), btmi(1:gs%Mice+1) !Bathymetries
type(bathym),            intent(in) :: btmdi(1:gs%Mice+1), btms(1:gs%nsoilcols+1) ! Bathymetries 
type(gridspacing_type),  intent(in) :: gsp ! Grid spacing

real(kind=ireals), intent(in) :: soilflux(1:gs%nsoilcols) ! Heat flux into soil, W/m**2

logical, intent(in) :: onlywater ! Indicates whether the lateral heat influx is calculated for water column only

! Input/output variables

! Output variables
type(lsh_type), intent(inout) :: lsh ! Heat source from soil columns for water, ice and deep ice

! Local variables
real(kind=ireals) :: offset
real(kind=ireals), allocatable :: z(:), area(:)
integer(kind=iintegers) :: i, j, i0, i1


! Water layer exists
if (ls%h1 > small_value) then

  offset = maxval(gsp%zsoilcols) - ls%h1 - STEP(ls%ls1 - small_value)*ls%ls1
  allocate ( z(0:gs%M+1), area(0:gs%M+1) )
  z(0) = offset ; area(0) = btmw(1)%area_int
  z(1:gs%M) = offset + gsp%z_half(1:gs%M) ; area(1:gs%M) = btmw(1:gs%M)%area_half
  z(gs%M+1) = offset + ls%h1 ; area(gs%M+1) = btmw(gs%M+1)%area_int

  lsh%water(:) = 0.
  waterscan : do j = 1, gs%nsoilcols-1

    if ( OVERLAY(z(0),z(gs%M+1),gsp%zsoilcols(j,ix,iy),gsp%zsoilcols(j+1,ix,iy)) ) then

      if (z(0) <= gsp%zsoilcols(j,ix,iy)) then
        i0 = -1
        do i = 1, gs%M+1
          if ( z(i-1) <  gsp%zsoilcols(j,ix,iy) .and. &
          &    z(i)   >  gsp%zsoilcols(j,ix,iy) ) then
            i0 = i
            exit
          endif
        enddo
      else
        i0 = -1
      endif
      if (i0 /= -1) then
        lsh%water(i0) = lsh%water(i0) + 1./btmw(i0)%area_int * &
        & ( area(i0) - btms(j)%area_int ) / (ls%h1*gsp%ddz05(i0-1)) * soilflux(j)
      endif 

      if (z(gs%M+1) >= gsp%zsoilcols(j+1,ix,iy)) then
        i1 = -1
        do i = 1, gs%M+1
          if ( z(i-1) <  gsp%zsoilcols(j+1,ix,iy) .and. &
          &    z(i)   >  gsp%zsoilcols(j+1,ix,iy) ) then
            i1 = i
            exit
          endif
        enddo
      else
        i1 = -1
      endif
      if (i1 /= -1) then
        lsh%water(i1) = + 1./btmw(i1)%area_int * &
        & ( btms(j+1)%area_int - area(i1-1) ) / (ls%h1*gsp%ddz05(i1-1)) * soilflux(j)
      endif

      forall ( i = 1:gs%M+1, z(i-1) >= gsp%zsoilcols(j,ix,iy) &
      & .and. z(i) <= gsp%zsoilcols(j+1,iy,iy) ) &
      & lsh%water(i) = + 1./btmw(i)%area_int * &
      & ( area(i) - area(i-1) )/ (ls%h1*gsp%ddz05(i-1)) * soilflux(j)

      if (z(gs%M+1) <= gsp%zsoilcols(j,ix,iy)) exit waterscan

    endif
  enddo waterscan

  deallocate(z,area)
endif

! Ice layer exists
if (ls%l1 > small_value .and. (.not. onlywater) ) then

  offset = maxval(gsp%zsoilcols) - ls%l1 &
  & - STEP(ls%ls1 - small_value)*ls%ls1 - STEP(ls%h1 - small_value)*ls%h1
  allocate ( z(0:gs%Mice+1), area(0:gs%Mice+1) )
  z(0) = offset ; area(0) = btmi(1)%area_int
  forall (i=1:gs%Mice+1) z(i) = z(i-1) + gsp%ddzi05(i-1)*ls%l1
  area(1:gs%Mice) = btmi(1:gs%Mice)%area_half  
  area(gs%Mice+1) = btmi(gs%Mice+1)%area_int

  lsh%ice(:) = 0.
  icescan : do j = 1, gs%nsoilcols-1

    if ( OVERLAY(z(0),z(gs%Mice+1),gsp%zsoilcols(j,ix,iy),gsp%zsoilcols(j+1,ix,iy)) ) then

      if (z(0) <= gsp%zsoilcols(j,ix,iy)) then
        i0 = -1
        do i = 1, gs%Mice+1
          if ( z(i-1) <  gsp%zsoilcols(j,ix,iy) .and. &
          &    z(i)   >  gsp%zsoilcols(j,ix,iy) ) then
            i0 = i
            exit
          endif
        enddo
      else
        i0 = -1
      endif
      if (i0 /= -1) then
        lsh%ice(i0) = lsh%ice(i0) + 1./btmi(i0)%area_int * &
        & ( area(i0) - btms(j)%area_int ) / (ls%l1*gsp%ddzi05(i0-1)) * soilflux(j)
      endif 

      if (z(gs%Mice+1) >= gsp%zsoilcols(j+1,ix,iy)) then
        i1 = -1
        do i = 1, gs%Mice+1
          if ( z(i-1) <  gsp%zsoilcols(j+1,ix,iy) .and. &
          &    z(i)   >  gsp%zsoilcols(j+1,ix,iy) ) then
            i1 = i
            exit
          endif
        enddo
      else
        i1 = -1
      endif
      if (i1 /= -1) then
        lsh%ice(i1) = + 1./btmi(i1)%area_int * &
        & ( btms(j+1)%area_int - area(i1-1) ) / (ls%l1*gsp%ddzi05(i1-1)) * soilflux(j)
      endif

      forall ( i = 1:gs%Mice+1, z(i-1) >= gsp%zsoilcols(j,ix,iy) &
      & .and. z(i) <= gsp%zsoilcols(j+1,iy,iy) ) &
      & lsh%ice(i) = + 1./btmi(i)%area_int * &
      & ( area(i) - area(i-1) )/ (ls%l1*gsp%ddzi05(i-1)) * soilflux(j)

      if ( z(gs%Mice+1) <= gsp%zsoilcols(j,ix,iy) ) exit icescan

    endif
  enddo icescan

  deallocate(z,area)
endif

! Deep ice layer exists
if (ls%ls1 > small_value .and. (.not. onlywater) ) then

  offset = maxval(gsp%zsoilcols) - ls%ls1
  allocate ( z(0:gs%Mice+1), area(0:gs%Mice+1) )
  z(0) = offset ; area(0) = btmdi(1)%area_int
  forall (i=1:gs%Mice+1) z(i) = z(i-1) + gsp%ddzi05(i-1)*ls%ls1
  area(1:gs%Mice) = btmdi(1:gs%Mice)%area_half  
  area(gs%Mice+1) = btmdi(gs%Mice+1)%area_int

  lsh%dice(:) = 0.
  dicescan : do j = 1, gs%nsoilcols-1

    if ( OVERLAY(z(0),z(gs%Mice+1),gsp%zsoilcols(j,ix,iy),gsp%zsoilcols(j+1,ix,iy)) ) then

      if (z(0) <= gsp%zsoilcols(j,ix,iy)) then
        i0 = -1
        do i = 1, gs%Mice+1
          if ( z(i-1) <  gsp%zsoilcols(j,ix,iy) .and. &
          &    z(i)   >  gsp%zsoilcols(j,ix,iy) ) then
            i0 = i
            exit
          endif
        enddo
      else
        i0 = -1
      endif
      if (i0 /= -1) then
        lsh%dice(i0) = lsh%dice(i0) + 1./btmdi(i0)%area_int * &
        & ( area(i0) - btms(j)%area_int ) / (ls%ls1*gsp%ddzi05(i0-1)) * soilflux(j)
      endif 

      if (z(gs%Mice+1) >= gsp%zsoilcols(j+1,ix,iy)) then
        i1 = -1
        do i = 1, gs%Mice+1
          if ( z(i-1) <  gsp%zsoilcols(j+1,ix,iy) .and. &
          &    z(i)   >  gsp%zsoilcols(j+1,ix,iy) ) then
            i1 = i
            exit
          endif
        enddo
      else
        i1 = -1
      endif
      if (i1 /= -1) then
        lsh%dice(i1) = + 1./btmdi(i1)%area_int * &
        & ( btms(j+1)%area_int - area(i1-1) ) / (ls%ls1*gsp%ddzi05(i1-1)) * soilflux(j)
      endif

      forall ( i = 1:gs%Mice+1, z(i-1) >= gsp%zsoilcols(j,ix,iy) &
      & .and. z(i) <= gsp%zsoilcols(j+1,iy,iy) ) &
      & lsh%dice(i) = + 1./btmdi(i)%area_int * &
      & ( area(i) - area(i-1) )/ (ls%ls1*gsp%ddzi05(i-1)) * soilflux(j)

      if ( z(gs%Mice+1) <= gsp%zsoilcols(j,ix,iy) ) exit dicescan

    endif
  enddo dicescan

  deallocate(z,area)
endif

contains
FUNCTION OVERLAY(x1,x2,y1,y2)
! Function OVERLAY checks, if two intervals overlay
! x2>x1, y2>y1
implicit none

logical :: OVERLAY
real(kind=ireals), intent(in) :: x1, x2, y1, y2 

OVERLAY = (y1 >= x1 .and. y2 <= x2) .or. &
        & (y1 <= x1 .and. y2 >= x2) .or. &
        & (y1 <= x1 .and. y2 >  x1) .or. &
        & (y1 <  x2 .and. y2 >= x2)

END FUNCTION OVERLAY
END SUBROUTINE LATERHEAT


END MODULE SOIL_MOD

