MODULE INIT_VAR_MOD
use LAKE_DATATYPES, only : ireals, iintegers
contains
SUBROUTINE INIT_VAR &
( M, Mice, ns, ms, ml, nsoilcols, ndatamax, &
& init_T, skin, zero_model, &
& h10, l10, hs10, ls10, tempair, Ts0, Tb0, Tm, Tbb0, &
& Sals0, Salb0, us0, vs0, h_ML0, &
& rosoil, rosdry, por, depth_area, ip, isp, &
& flag_snow, flag_snow_init, itop, nstep, itherm, &
& h1, l1, hs1, ls1, &
& hx1, hx2, hy1, hy2, &
& hx1t, hx2t, hy1t, hy2t, &
& hx1ml, hx2ml, hy1ml, hy2ml, &
& veg, snmelt, snowmass, cdmw, velfrict_prev, &
& roughness, eflux0_kinem, Elatent, &
& totalevap, totalmelt, totalprecip, totalwat, totalpen, &
& time, dhwfsoil, dhw, dhw0, dhi, dhi0, dls0, &
& E1, eps1, zsoilcols, Tsoil1, wi1, wl1, Sals1, &
& rootss, qsoil, TgrAnn, qwater, &
& oxyg, DIC, phosph, oxygsoil, Sal1, u1, v1, Tw1, Tskin, T_0dim, Eseiches, &
& Ti1, salice, porice, Tis1, z_full, ddz, &
& dz, T, wl, dens, lamw, &
& dzeta_int, zsoil, pressure)
! Subroutine INIT_VAR initializes the prognostic variables of the model
use PHYS_FUNC, only : &
& MELTPNT, &
& UNFRWAT, &
& WL_MAX, &
& WI_MAX, &
& HENRY_CONST, &
& MELTINGPOINT
use PHYS_CONSTANTS, only : &
& row0, g, Kelvin0, R_univ, poricebot
use METH_OXYG_CONSTANTS, only : &
& ch4_atm0, o2_atm0, co2_atm0, &
& rel_conc_ebul_crit, &
& Henry_const0_ch4, &
& Henry_temp_dep_ch4, &
& Henry_temp_ref, &
& molmass_p
use DRIVING_PARAMS, only : &
& tricemethhydr, &
& nmeltpoint, &
& initprof ! derived datatype
use NUMERIC_PARAMS, only : &
& small_value
implicit none
real(kind=ireals), save :: E_init
real(kind=ireals), save :: eps_init
! Input variables
integer, intent(in) :: M, Mice, ns, ms, ml, nsoilcols, ndatamax ! Grid dimensions
integer, intent(in) :: init_T, skin, zero_model ! Driving parameters
real(kind=ireals), intent(in) :: h10, l10, hs10, ls10
real(kind=ireals), intent(in) :: tempair, Ts0, Tb0, Tm, Tbb0
real(kind=ireals), intent(in) :: Sals0, Salb0
real(kind=ireals), intent(in) :: us0, vs0
real(kind=ireals), intent(in) :: h_ML0
real(kind=ireals), intent(in) :: rosoil(1:ns), rosdry(1:ns), por(1:ns)
real(kind=ireals), intent(in) :: depth_area(1:ndatamax,1:2) ! Data for lake bathymetry
type(initprof), intent(in) :: ip, isp
real(kind=ireals), intent(in) :: zsoil(1:ns)
real(kind=ireals), intent(in) :: pressure
real(kind=ireals), intent(in) :: dzeta_int(1:M+1)
real(kind=ireals), intent(in) :: ddz(1:M)
! Output variables
integer(kind=iintegers), intent(out) :: flag_snow, flag_snow_init, itop
integer(kind=iintegers), intent(out) :: nstep
integer(kind=iintegers), intent(inout), allocatable :: itherm(:)
real(kind=ireals), intent(out) :: h1, l1, hs1, ls1
real(kind=ireals), intent(out) :: hx1, hx2, hy1, hy2
real(kind=ireals), intent(out) :: hx1t, hx2t, hy1t, hy2t
real(kind=ireals), intent(inout), allocatable :: hx1ml(:), hx2ml(:), hy1ml(:), hy2ml(:)
real(kind=ireals), intent(out) :: veg
real(kind=ireals), intent(out) :: snmelt, snowmass
real(kind=ireals), intent(out) :: cdmw, velfrict_prev
real(kind=ireals), intent(out) :: roughness
real(kind=ireals), intent(out) :: eflux0_kinem, Elatent
real(kind=ireals), intent(out) :: totalevap, totalmelt, totalprecip, totalwat, totalpen
real(kind=ireals), intent(out) :: time
real(kind=ireals), intent(out) :: dhwfsoil
real(kind=ireals), intent(out) :: dhw, dhw0, dhi, dhi0, dls0
real(kind=ireals), intent(out) :: E1(1:M+1), eps1(1:M+1)
real(kind=ireals), intent(out) :: zsoilcols(1:nsoilcols+1)
real(kind=ireals), intent(out) :: Tsoil1(1:ns,1:nsoilcols), Sals1(1:ns,1:nsoilcols)
real(kind=ireals), intent(out) :: wi1 (1:ns,1:nsoilcols), wl1 (1:ns,1:nsoilcols)
real(kind=ireals), intent(out) :: rootss(1:ns), qsoil(1:ns,1:nsoilcols), TgrAnn(1:ns), qwater(1:M+1)
real(kind=ireals), intent(out) :: oxyg(1:M+1), DIC(1:M+1), phosph(1:M+1)
real(kind=ireals), intent(out) :: oxygsoil(1:nsoilcols)
real(kind=ireals), intent(out) :: Sal1(1:M+1)
real(kind=ireals), intent(out) :: u1(1:M+1), v1(1:M+1)
real(kind=ireals), intent(out) :: Tw1(1:M+1), Tskin(1:2)
real(kind=ireals), intent(out) :: Ti1(1:Mice+1), Tis1(1:Mice+1)
real(kind=ireals), intent(out) :: salice(1:Mice+1), porice(1:Mice+1)
real(kind=ireals), intent(out) :: z_full(1:M+1)
real(kind=ireals), intent(out) :: dz(1:ms), T(1:ml), wl(1:ml), dens(1:ms)
real(kind=ireals), intent(out) :: lamw(1:M)
real(kind=ireals), intent(out) :: T_0dim
real(kind=ireals), intent(out) :: Eseiches
! Local variables
real(kind=ireals), parameter :: Ct_d_Ch = 2. ! Using results by West & Plug, 2007, may be estimated
! from 1.7 to 2.8
real(kind=ireals) :: work, work1
real(kind=ireals), allocatable :: pressoil(:), workarr(:)
real(kind=ireals) :: z, h_ML0zv, h_talik, ch4_bot
real(kind=ireals) :: Ti10, Ti11
integer(kind=iintegers) :: i, j, i_ML, i_talik, k
integer(kind=iintegers) :: n_1cm, n_5cm
logical :: flag
logical, save :: firstcall = .true.
logical, parameter :: soil_thermokarst_init = .false.
real(kind=ireals), parameter :: relch4_0 = 1.
if (firstcall) then
E_init = 10.d0**(-5.5)
eps_init = 1.d-9
endif
h1 = h10
hx1 = 0.; hx2 = 0.
hy1 = 0.; hy2 = 0.
hx1t = 0.; hx2t = 0.
hy1t = 0.; hy2t = 0.
if (allocated(hx1ml)) then
hx1ml(1:M+1) = 0.; hx2ml(1:M+1) = 0.
hy1ml(1:M+1) = 0.; hy2ml(1:M+1) = 0.
endif
if (allocated(itherm)) itherm(1:M+2) = -1
l1 = l10
hs1 = hs10
ls1 = ls10
do i = 1, M
E1(i) = E_init !+ float(i)/float(M)*(1.d-10 - 1.d-8)
eps1(i) = eps_init !+ float(i)/float(M)*(1.d-18 - 1.d-14)
enddo
do i = 1, M+1
z_full (i) = dzeta_int(i)*h1
enddo
i_ML = int((M+1)*h_ML0/h1) ! The index of level of mixed-layer depth
! Water temperature initialization
if (init_T == 2) then
h_ML0zv = (Tm - 0.5*(Tb0 + Ts0))/(Ts0/h1 - 0.5*(Ts0 + Tb0)/h1)
i_ML = int(M*h_ML0zv/h1)
endif
if (init_T == 1 .or. init_T == 2) then
Tw1(1:max(i_ML,1)) = Ts0
do i = max(i_ML+1,2), M+1
Tw1(i) = Ts0 + (Tb0-Ts0)*float(i-i_ML)/float(M+1-i_ML)
enddo
! The initial temperature profile is given from the input file
elseif (init_T == 3) then
call LININTERPOL (ip%zTinitprof,ip%Tinitprof,ip%lenprof, &
& z_full,Tw1,M+1,flag,.true.)
call LININTERPOL (ip%zTinitprof,ip%Sinitprof,ip%lenprof, &
& z_full,Sal1,M+1,flag,.true.)
if (.not.flag) then
print*, 'The error while interpolating the initial &
&temperature profile: terminating program'
STOP
endif
endif
! Lake depth intervals for the tops of soil columns.
! An alternative - to calculate from the vertical grid of depth-area data.
if (nsoilcols > 1 .and. depth_area(1,2) >= 0.) then
!work = h10/real(nsoilcols)
work1 = maxval(depth_area(:,1))
work = ( work1 - h10*0.5*ddz(M) )/real(nsoilcols-1)
zsoilcols(nsoilcols+1) = work1
do i = 1, nsoilcols
zsoilcols(i) = (i-1)*work
enddo
endif
Sals1(1:ns,1:nsoilcols) = Salb0*row0/( rosdry(1)*(1 - por(1)) ) ! assuming, that salinity in ground
! is the same
! as one in near bottom layer of water
allocate(pressoil(1:ns))
do i = 1, ns
pressoil(i) = pressure + row0*g*(h10 + zsoil(i))
enddo
! Temperature distribution in talik assuming homogeneous distribution of salinity in ground
flag = (h10 > 0. .and. Tb0 > MELTINGPOINT(Sals1(1,nsoilcols),pressoil(1),tricemethhydr%par) &
& .and. Tbb0 < MELTINGPOINT(Sals1(ns,nsoilcols),pressoil(ns),tricemethhydr%par))
if (flag .and. soil_thermokarst_init) then
! Initializing the talik depth, assuming thermokarst lake, following West and Plug, 2007
h_talik = min(Ct_d_Ch*h10,0.8*zsoil(ns))
do i = 1, ns
if (zsoil(i) >= h_talik) then
i_talik = i ! i_talik is the number of the first level below the estimated talik depth
exit
endif
enddo
! Temperature distribution in the talik
work = MELTINGPOINT(Sals1(i_talik,nsoilcols),pressoil(i_talik),tricemethhydr%par)
do i = 1, i_talik
Tsoil1(i,nsoilcols) = Tb0 + float(i-1)/float(i_talik-1) * (work - Tb0)
enddo
! Temperature distribution below the talik
do i = i_talik+1, ns
Tsoil1(i,nsoilcols) = work + float(i-i_talik)/float(ns-i_talik)*(Tbb0 - work)
enddo
else
h_talik = 0.
! Linear temperature profile in soil
if (init_T == 1 .or. init_T == 2) then
do i = 1, ns
Tsoil1(i,nsoilcols) = Tb0 + float(i-1)/float(ns-1)*(Tbb0-Tb0)
enddo
! The initial soil temperature profile is given from the input file
elseif (init_T == 3) then
call LININTERPOL (isp%zTinitprof,isp%Tinitprof,isp%lenprof, &
& zsoil,Tsoil1(1,nsoilcols),ns,flag,.true.)
if (.not.flag) then
print*, 'The error while interpolating the initial &
&soil temperature profile: terminating program'
STOP
endif
endif
endif
! Ice and liquid water content initialization in soil
do i = 1, ns
if (Tsoil1(i,nsoilcols) > MELTINGPOINT(Sals1(i,nsoilcols),pressoil(i),tricemethhydr%par)) then
wi1(i,nsoilcols) = 0.
wl1(i,nsoilcols) = WL_MAX(por(i),rosoil(i),0.e0_ireals,tricemethhydr%par) - 0.01
else
wl1(i,nsoilcols) = UNFRWAT(Tsoil1(i,nsoilcols),i)
wi1(i,nsoilcols) = WI_MAX(por(i),rosoil(i),tricemethhydr%par) - wl1(i,nsoilcols) - 0.01
endif
enddo
! Identical initialization of all soil columns
if (nsoilcols > 1 .and. depth_area(1,2) >= 0.) then
do i = 1, ns
wl1 (i,1:nsoilcols-1) = wl1 (i,nsoilcols)
wi1 (i,1:nsoilcols-1) = wi1 (i,nsoilcols)
enddo
!Temperature in side soil columns is initialized as a linear function from
!a water temperature at respective depth to Tbb0
allocate(workarr(1:M+1))
do j = 1, nsoilcols-1
work = 0.5*(zsoilcols(j) + zsoilcols(j+1))
workarr(1:M+1) = abs(z_full(1:M+1) - work)
k = minloc(workarr,1)
do i = 1, ns
!Tsoil1(i,j) = Tw1(k) + float(i-1)/float(ns-1)*(Tbb0 - Tw1(k))
Tsoil1(i,j) = Tw1(k) + (Tbb0 - Tw1(k))*(1. - exp( - 10.*float(i-1)/float(ns-1) ) )
enddo
enddo
deallocate(workarr)
endif
! Methane variables initialization
rootss(1:ns) = 0.e0_ireals ! No roots
veg = 0.e0_ireals ! No vegetation at the bottom
!i_ML = int((M+1)*h_ML0/h1)
if (init_T /= 3) then
ch4_bot = relch4_0*rel_conc_ebul_crit*(pressure + row0*g*h10)* &
& HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, &
& Henry_temp_ref, Tb0 + Kelvin0) ! Saturated concentration at the bottom
do i = 1, M+1
if (i <= i_ML) then
qwater(i) = ch4_atm0 ! Atmospheric concentration
else
qwater(i) = qwater(max(i_ML,1)) + &
& (ch4_bot - qwater(max(i_ML,1)))*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
endif
enddo
else
call LININTERPOL (ip%zTinitprof,ip%ch4initprof,ip%lenprof, &
& z_full,qwater,M+1,flag,.true.)
endif
! qwater2(1:M+1,1,1:2) = 0.5*ch4_atm0 ! two-meth
TgrAnn(1:ns) = Tsoil1(1:ns,nsoilcols)
! qsoil(1) = ch4_bot
! qsoil2(1,1:2) = 0.5*ch4_atm0 ! two-meth
! Initialization of methane in soil columns at the lake's slopes
if (nsoilcols > 1 .and. depth_area(1,2) >= 0.) then
do j = 1, nsoilcols-1
do i = 1, ns
qsoil(i,j) = relch4_0*rel_conc_ebul_crit*por(i) * &
& (pressure + row0*g*(0.5*(zsoilcols(j) + zsoilcols(j+1)) + zsoil(i))) * & ! assuming hydrostatic pressure in soil
& HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, &
& Henry_temp_ref, Tsoil1(i,j) + Kelvin0)
! qsoil2(i,1:2) = 0.5*qsoil(i) ! two-meth
enddo
enddo
endif
! Initialization of methane in deepest soil column
qsoil(1,nsoilcols) = qwater(M+1)*por(1) ! Boundary condition
do i = 2, ns
qsoil(i,nsoilcols) = relch4_0*rel_conc_ebul_crit*por(i) * &
& (pressure + row0*g*(h10 + zsoil(i))) * & ! assuming hydrostatic pressure in soil
& HENRY_CONST(Henry_const0_ch4, Henry_temp_dep_ch4, &
& Henry_temp_ref, Tsoil1(i,nsoilcols) + Kelvin0)
! qsoil2(i,1:2) = 0.5*qsoil(i) ! two-meth
enddo
!qwater(M+1) = qsoil(1,nsoilcols)/por(1) ! Boundary condition
! Dissolved inorganic carbon initialization
if (init_T /= 3) then
DIC(1:M+1) = 10.*co2_atm0 ! Atmospheric concentration
else
call LININTERPOL (ip%zTinitprof,ip%co2initprof,ip%lenprof, &
& z_full,DIC,M+1,flag,.true.)
endif
! Inorganic dissolved phosphorus
if (init_T /= 3) then
phosph(:) = 0.01/molmass_p
else
call LININTERPOL (ip%zTinitprof,ip%phosphinitprof,ip%lenprof, &
& z_full,phosph,M+1,flag,.true.)
endif
! Oxygen initialization
if (init_T /= 3) then
do i = 1, M+1
if (i <= i_ML) then
oxyg(i) = o2_atm0 ! Atmospheric concentration
else
oxyg(i) = oxyg(max(i_ML,1)) - oxyg(max(i_ML,1))*(h_ML0 - dzeta_int(i)*h10)/(h_ML0 - h10)
endif
oxyg(i) = max(oxyg(i),small_value)
enddo
else
call LININTERPOL (ip%zTinitprof,ip%o2initprof,ip%lenprof, &
& z_full,oxyg,M+1,flag,.true.)
endif
oxygsoil(1:nsoilcols) = 0. !Zero initial oxygen concentration in the aerobic soil layer
u1(1:max(i_ML,1)) = us0
v1(1:max(i_ML,1)) = vs0
do i = max(i_ML+1,2), M+1
u1(i) = us0 + (0.e0_ireals - us0)*float(i-i_ML)/float(M+1-i_ML)
v1(i) = vs0 + (0.e0_ireals - vs0)*float(i-i_ML)/float(M+1-i_ML)
enddo
if (init_T /= 3) then
Sal1(1:max(i_ML,1)) = Sals0
do i = max(i_ML+1,2), M+1
Sal1(i) = Sals0 + (Salb0-Sals0)*float(i-i_ML)/float(M+1-i_ML)
enddo
endif
if (skin == 0) then
Tskin(1:2) = 0.e0_ireals
else
Tskin(1) = Tw1(1)
endif
! Initial temperature for zero-dimensional model
if (zero_model == 1) then
T_0dim = Tw1(1)*0.5*ddz(1) + Tw1(M+1)*0.5*ddz(M)
do i = 2, M
T_0dim = T_0dim + Tw1(i)*0.5*(ddz(i-1) + ddz(i) )
enddo
else
T_0dim = 0.
endif
Ti1 = 0.
salice(:) = 0.
porice(:) = poricebot
Tis1 = 0.
Ti10 = min(tempair,-1.d-1) ! Initial ice surface temperature, Celsius
if (l1 /= 0) then
if (h1 > 0.) then
! Ice over water
Ti11 = MELTPNT(Sals0,0.e0_ireals,nmeltpoint%par)
else
! Ice over soil
Ti11 = Tb0
endif
do i = 1, Mice+1
Ti1(i)= Ti10 + (Ti11 - Ti10)*float(i-1)/float(Mice) ! Linear profile, if dzetai-grid is regular,
enddo ! water layer underneath is assumed to exist
endif
if (hs1 == 0.) then
flag_snow = 0
flag_snow_init = 1
itop = 1
else
flag_snow = 1
flag_snow_init = 0
hs1 = max(2.d-2, hs1)
hs1 = int(hs1/0.01)*0.01
n_5cm = int(hs1/0.05)
n_1cm = int((hs1-n_5cm*0.05)/0.01)
if (n_1cm == 0) then
n_1cm = 1
hs1 = hs1 + 0.01
endif
itop = ms - (n_1cm + n_5cm)
do i = itop, itop + n_1cm - 1
dz(i) = 0.01
T(i) = min(tempair,-5.d0)
wl(i) = 0
dens(i) = 150.
enddo
do i = itop + n_1cm, ms-1
dz(i) = 0.05
T(i) = min(tempair,-5.d0)
wl(i) = 0
dens(i) = 150.
enddo
dz(ms) = 0.5d0*dz(ms-1)
T(ms) = min(tempair,-5.d0)
wl(ms) = 0.e0_ireals
dens(ms) = 150.d0
endif
snmelt = 0.e0_ireals
snowmass = 0.e0_ireals
cdmw = 1.d-15
velfrict_prev = 1.d-2
roughness = 1.d-3
eflux0_kinem = 0.e0_ireals
Elatent = 0.
totalevap = 0
totalmelt = 0
totalprecip = 0.
totalwat = 0
totalpen = 0
time = 0
nstep = 0
dhwfsoil = 0.
dhw = 0.
dhw0 = 0.
dhi = 0.
dhi0 = 0.
dls0 = 0.
lamw(:) = 1.d+3
Eseiches = 0. !initial seiche energy
deallocate(pressoil)
if (firstcall) firstcall = .false.
END SUBROUTINE INIT_VAR
END MODULE INIT_VAR_MOD