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), 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(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
max1 = 0.
do i = 1, ns
if (Tsoil3(i,is) < &
& 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
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,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
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%Tw1,wst%Ti1,wst%Tis1, &
& 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,1),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,Sals1, &
& 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, &
& 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,Sals1, &
& 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, &
& 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, &
& Tw1,Ti1,Tis1, &
& 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) :: Tw1(1:M+1)
real(kind=ireals), intent(in) :: Ti1(1:Mice+1), Tis1(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) = Ti1(1)
j = 2
! Ice interior
do i = 2, Mice
Temp(j) = Ti1(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) = Ti1(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) = Tw1(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) = Tw1(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) = Tis1(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) = Tis1(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) = Ti1(1)
j = 2
do i = 2, Mice
Temp(j) = Ti1(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) = Ti1(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) = Tw1(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) = Tw1(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) = Tw1(1)
j = 2
! Water interior
do i = 2, M
Temp(j) = Tw1(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) = Tw1(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) = Tis1(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) = Tis1(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) = Tw1(1)
j = 2
do i = 2, M
Temp(j) = Tw1(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) = Tw1(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
!print*, 'z0', z(0), ls%l1, gsp%ddzi05(gs%Mice), gsp%ddzi05(0)
do i = 1, gs%Mice+1
z(i) = z(i-1) + gsp%ddzi05(i-1)*ls%l1
enddo
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
do i = 1, gs%Mice+1
z(i) = z(i-1) + gsp%ddzi05(i-1)*ls%ls1
enddo
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