Skip to content
Snippets Groups Projects
soil_mod.f90 48.4 KiB
Newer Older
  • Learn to ignore specific revisions
  • Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
    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), 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
    
              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) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            endif
    
          enddo 
        endif
     !dhwfsoil = dhwfsoil - abs(wl4(i))*rosdry(i)*(1-por(i))*dzss(i)/row
    1   wl3(i)=0.
      endif
    enddo
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    ! STEP 3 OF SPLITTING-UP METHOD : PHASE PROCESSES
    mhsep = (1. + tricemethhydr%par*alphamh)
    20 c2: do i = 1, ns
    
      work = MELTINGPOINT(Sals1(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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
        & MELTINGPOINT(Sals1(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    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      
    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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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 / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      & (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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        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) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            & (rosdry(j)*(1-por(j))*dzss(j))) then
    
              wl4(j,is)=wl4(j,is)+surplus*rosdry(i)*(1-por(i))*dzss(i) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
              & (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)) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
              & (rosdry(i)*(1-por(i))*dzss(i))
    
              wl4(j,is) = &
              & POR(j)*ROW0/(rosdry(j)*(1-por(j))) - wi2(j,is)*row0_d_roi 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            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 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      endif
    
      if(wl4(i,is) < 0) then
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        if (i/=1) then 
          do j=i-1,1,-1
    
            if (wl4(j,is) > - wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            & (rosdry(j)*(1-por(j))*dzss(j))) then
    
              wl4(j,is) = wl4(j,is) + &
              & wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
              & (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) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
            & (rosdry(j)*(1-por(j))*dzss(j))) then
    
              wl4(j,is) = wl4(j,is) + &
              & wl4(i,is)*rosdry(i)*(1-por(i))*dzss(i) / &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
              & (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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    ! 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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    enddo
    
    deallocate(pressoil)
    
    if (firstcall) firstcall=.false.
    RETURN
    END SUBROUTINE SOILFORLAKE
    
    
    
    SUBROUTINE SOIL_COND_HEAT_COEF(is)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    implicit none
    
    
    integer(kind=iintegers), intent(in) :: is
    
    real(kind=ireals), parameter :: cr = 0.2*4180.d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    logical, save :: firstcall = .true.
    
    
    if (firstcall) then
      if (tricemethhydr%par > 0.) then
        cimh => cpmethhydr
      else
        cimh => ci
      endif
    endif
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
       
    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)) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !  ARG = (WL1(i)+WI1(i))/wlmax_i
    
    !  ARG = max(ARG,1.d-2)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !  PSI = PSIMAX(i)*(ARG)**(-BH(i))
    
    !  PF = log(-PSI*100.)/log(10.d0) 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !  if(PF>5.1) then
    !    lamsoil(i) = 4.1E-4*418.
    !  else
    
    !    lamsoil(i) = exp(-PF-2.7)*418.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    !  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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
        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)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      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) + &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      & filtr(i) ) *0.5*(rosdry(i)+rosdry(i+1))* &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    enddo 
    
    if (firstcall) firstcall = .false.
    END SUBROUTINE SOIL_COND_HEAT_COEF
    
    SUBROUTINE SOILCOLSTEMP(gs,gsp,dt,ls,ftot,ch4_pres_atm, &
    
                           & ddz,ddzi,zsoilcols, &
    
                           & bathymdice,bathymsoil, &
                           & soilflux,fdiffbot, contr) 
    
    
    ! Subroutine calculates temperature in soil columns,
    
    ! excepting for the lowest column, assuming 
    ! nsoilcols = M+1
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    use LAKE_DATATYPES, only : &
    & iintegers, ireals
    
    
    use DRIVING_PARAMS, only : soilcolconjtype
    
    
    use ATMOS, only : &
    & pressure, wind10
    
    
    & cw_m_row0, z0_bot
    
    & LOGFLUX, TEMPWATR
    
    use ARRAYS_BATHYM, only : bathym, dhw, dhw0, layers_type
    
    use ARRAYS_SOIL, only : csoil, rosoil, lamsoil, Tsoil3, Tsoil2, Tsoil1, &
    
    & wl4,wi2,wa,Sals1,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
    
    
    type(gridsize_type), intent(inout) :: gs
    
    type(gridspacing_type), intent(in) :: gsp
    
    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)
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    real(kind=ireals), intent(inout), dimension(1:vector_length) :: a, b, c, d
    
    logical, intent(inout) :: add_to_winter
    
    
    real(kind=ireals),      intent(inout) :: soilflux(1:gs%nsoilcols), fdiffbot(1:gs%nsoilcols)
    
    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(:)
    
    if ((ls%l1 /= 0. .or. ls%ls1 /= 0) .and. bcswitch == 2) then
    
      print*, "bcswitch == 2 is not operational when &
    
    
    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, &
    
                           & 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
    
          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)
    
    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, &
    
                           & bathymwater,bathymice,bathymdice,bathymsoil,&
                           & qsoilsurf)
          deallocate (work)
        endif
    
        allocate (work(0:gs%M+1))
        do is = 1, gs%nsoilcols-1
          gs%isoilcol = is
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,is),rosoil, &
    
          & rootss,rosdry,por,veg,qsoilsurf(is)*por(1),TgrAnn, methgenmh, &
    
          & ddz,ddz05, wst, lammeth, &
          & gsp%z_half(bathymsoil(is)%icent), & !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, &
    
          & 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
    
          ! 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
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          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), &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          & (gas,pressure,wind10,ch4_pres_atm,zsoil,Tsoil3(1,is),rosoil, &
    
          & rootss,rosdry,por,veg,fdiffbot(is),TgrAnn, methgenmh, &
    
          & ddz,ddz05, wst, lammeth, &
          & gsp%z_half(bathymsoil(is)%icent), & !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, &
    
          & 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)