diff --git a/source/model/bathym_mod.f90 b/source/model/bathym_mod.f90
index c41531f5321491d581db93f7a480d724a0d4e165..461f3194ba1f4fbd3a566cd49e89981bd606c99b 100644
--- a/source/model/bathym_mod.f90
+++ b/source/model/bathym_mod.f90
@@ -170,6 +170,7 @@ else
enddo
enddo
bathymsoil(nsoilcols,ix,iy)%ibot = M+1
+ bathymsoil(:,ix,iy)%icent = 1
do i = 1, nsoilcols-1
do j = 1, M-1
if (work3(i) + dz > z_half(j) .and. &
diff --git a/source/model/lake_modules.f90 b/source/model/lake_modules.f90
index f834e22d854f785ef62203fe937bfae49fedbf47..ac83e7e50857a19c238f46aef37e416fda226d08 100644
--- a/source/model/lake_modules.f90
+++ b/source/model/lake_modules.f90
@@ -2025,8 +2025,8 @@ k_o2 = 0.672 / molm3tomgl_o2 ! (Lidstrom and Somers, 1984)
delta_Eq = 5.d+4 ! activation energy for methane oxidation, J/mol,
! after Arah&Stephen (1998), Dunfield et al. (1993),
! Nedwell&Watson (1995)
-Vmaxw = 2.5d-2/86400. ! reaction potential in oxygen-saturated Michaelis-Menten kinetics,
- ! after (Liikanen et al., 2002)
+Vmaxw = 1.d-2/86400. ! reaction potential in oxygen-saturated Michaelis-Menten kinetics,
+ ! 1.d-1/86400. after (Liikanen et al., 2002)
k_ch40 = 1.d-2 ! half-saturation constant in oxygen-saturated Michaelis-Menten kinetics,
! after (Liikanen et al., 2002)
diff --git a/source/model/oxygen_mod.f90 b/source/model/oxygen_mod.f90
index 28f48b513330f7213d72c791559821ddee52ab0a..87d69b5e04a1d19fb0f80d8278124eae9744c1b8 100644
--- a/source/model/oxygen_mod.f90
+++ b/source/model/oxygen_mod.f90
@@ -655,6 +655,10 @@ if (carbon_model%par == 2) then
Pd_DOC = rad_CDOM * &
& (1. - exp( - RadWater%extinct_CDOC(i)*gsp%ddz05(i-1)*ls%h1)) / &
& gsp%ddz05(i-1)*ls%h1 * AQY
+ !Pd_DOC is subtracted from DOC but not taken into account in other
+ !equations, as CO is the main product of photodegradation and there is
+ !no enough evidence to assume it is fully added to DIC and not emitted to the atmosphere
+ !(Zuo and Jones, Water Research, 1997)
!print*, rad_CDOM, &
!& (1. - exp( - RadWater%extinct_CDOC(i)*gsp%ddz05(i-1)*ls%h1))/ &
!& gsp%ddz05(i-1)*ls%h1, AQY
diff --git a/source/model/soil_mod.f90 b/source/model/soil_mod.f90
index 4b3f3f6d8ad5c5ddf74059fb6326402a79056a32..9c207b8eb13ae377b6285174297d057884d2724c 100644
--- a/source/model/soil_mod.f90
+++ b/source/model/soil_mod.f90
@@ -1,543 +1,543 @@
- 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
+MODULE SOIL_MOD
- !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
+use LAKE_DATATYPES, only : ireals, iintegers
+use NUMERICS, only : PROGONKA, STEP
+use NUMERIC_PARAMS, only : vector_length
- allocate (pressoil(1:ns))
- do i = 1, ns
- pressoil(i) = pressure + row0*g*(h1 + zsoil(i))
- enddo
+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
- if (tricemethhydr%par > 0.) then
- hicemelt => methhydrdiss
+!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
- 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)
+ 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
-
- 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
+ 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
-
- 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
+ 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
- filtr_low = 0
- wl3(ns) = wl2(ns) + dt*(-filtr_low+filtr(ns-1))/dzss(ns)
-
- do i=1,ns
- if (wl3(i)>WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)) then
- surplus=wl3(i)-WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)
- cy1:do j=1,ns
- if (WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par)-wl3(j)>0) then
- if (WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par) - wl3(j) > &
- & surplus*rosdry(i)*(1-por(i))*dzss(i) / &
- & (rosdry(j)*(1-por(j))*dzss(j))) then
- wl3(j)=wl3(j)+surplus*rosdry(i)*(1-por(i))*dzss(i) / &
- & (rosdry(j)*(1-por(j))*dzss(j))
- surplus=0
- exit cy1
- else
- surplus = surplus - &
- & (WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par)-wl3(j)) * &
- & (rosdry(j)*(1-por(j))*dzss(j)) / &
- & (rosdry(i)*(1-por(i))*dzss(i))
- wl3(j)=WL_MAX(por(j),rosdry(j),wi1(j,is),tricemethhydr%par)
- endif
- endif
- enddo cy1
- dhwfsoil3 = dhwfsoil3 + surplus*rosdry(i)*(1-por(i))*dzss(i)/row0
- wl3(i) = WL_MAX(por(i),rosdry(i),wi1(i,is),tricemethhydr%par)
- endif
- if(wl3(i)<0) then
- if (i/=1) then
- do j=i-1,1,-1
- if (wl3(j)>-wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
- & (rosdry(j)*(1-por(j))*dzss(j))) then
- wl3(j)=wl3(j)+wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
- & (rosdry(j)*(1-por(j))*dzss(j))
- goto 1
- endif
- enddo
- endif
- if (i/=ns) then
- do j=i+1,ns
- if (wl3(j)>-wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
- & (rosdry(j)*(1-por(j))*dzss(j))) then
- wl3(j)=wl3(j)+wl3(i)*rosdry(i)*(1-por(i))*dzss(i) / &
- & (rosdry(j)*(1-por(j))*dzss(j))
- goto 1
- endif
- enddo
- endif
- !dhwfsoil = dhwfsoil - abs(wl4(i))*rosdry(i)*(1-por(i))*dzss(i)/row
- 1 wl3(i)=0.
- endif
- enddo
+! STEP 3 OF SPLITTING-UP METHOD : PHASE PROCESSES
+mhsep = (1. + tricemethhydr%par*alphamh)
+20 c2: do i = 1, ns
+ work = MELTINGPOINT(Sals2(i,is)/wl3(i),pressoil(i),tricemethhydr%par)
+ if (wi1(i,is) > 0 .and. Tsoil2(i) > work + 0.01) then
+ Potphenergy = (Tsoil2(i) - (work + 0.01)) * &
+ & csoil(i)*rosoil(i)*dzss(i)
+ if (Potphenergy >= wi1(i,is)*rosdry(i)*dzss(i)*(1-por(i))*hicemelt) then
+ wl4(i,is) = wl3(i) + wi1(i,is)/mhsep
+ wi2(i,is) = 0.d0
+ Tsoil3(i,is) = work + 0.01 + &
+ & (Potphenergy-(wi1(i,is)*rosdry(i)*dzss(i)* &
+ & (1 - por(i))*hicemelt))/(csoil(i)*rosoil(i)*dzss(i))
+ else
+ wl4(i,is) = wl3(i) + Potphenergy / &
+ & (rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)/mhsep
+ wi2(i,is) = wi1(i,is) - Potphenergy / &
+ & (rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)
+ Tsoil3(i,is) = work + 0.01
+ endif
+ else
+ if (wl3(i) >= 0 .and. Tsoil2(i) < work - 0.01) then
+ Potphenergy = - (Tsoil2(i) - work + 0.01) * &
+ & csoil(i)*rosoil(i)*dzss(i)
+ watavailtofreeze = (wl3(i) - UNFRWAT(Tsoil2(i),i)) * &
+ & rosdry(i)*dzss(i)*(1 - por(i))*mhsep*hicemelt
+ if (watavailtofreeze < 0) then
+ !cc = csoil(i)*rosoil(i)*dzss(i)
+ !cc1 = rosdry(i)*dzss(i)*(1-por(i))*Lwi
+ !bb = (-Potphenergy-0.1*cc-wl3(i)*cc1)/cc
+ !aa = cc1/cc
+ !call phase_iter(aa,bb,i,Tsoil3(i))
+ Tsoil3(i,is) = Tsoil2(i) + watavailtofreeze/(csoil(i)*rosoil(i)*dzss(i))
+ wi2(i,is) = wi1(i,is) + (wl3(i) - UNFRWAT(Tsoil2(i),i))*mhsep
+ wl4(i,is) = UNFRWAT(Tsoil2(i),i)
+ cycle c2
+ endif
+ if (Potphenergy >= watavailtofreeze) then
+ Tsoil3(i,is) = work - 0.01 - &
+ & (Potphenergy - watavailtofreeze) / &
+ & (csoil(i)*rosoil(i)*dzss(i))
+ wi2(i,is) = wi1(i,is) + (wl3(i) - UNFRWAT(Tsoil2(i),i))*mhsep
+ wl4(i,is) = UNFRWAT(Tsoil2(i),i)
+ else
+ wl4(i,is) = wl3(i) - Potphenergy/(rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)/mhsep
+ wi2(i,is) = wi1(i,is) + Potphenergy/(rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)
+ Tsoil3(i,is) = work - 0.01
+ endif
+ else
+ wl4(i,is) = wl3(i)
+ wi2(i,is) = wi1(i,is)
+ Tsoil3(i,is) = Tsoil2(i)
+ endif
+ endif
+enddo c2
- ! STEP 3 OF SPLITTING-UP METHOD : PHASE PROCESSES
- mhsep = (1. + tricemethhydr%par*alphamh)
- 20 c2: do i = 1, ns
- work = MELTINGPOINT(Sals2(i,is)/wl3(i),pressoil(i),tricemethhydr%par)
- if (wi1(i,is) > 0 .and. Tsoil2(i) > work + 0.01) then
- Potphenergy = (Tsoil2(i) - (work + 0.01)) * &
- & csoil(i)*rosoil(i)*dzss(i)
- if (Potphenergy >= wi1(i,is)*rosdry(i)*dzss(i)*(1-por(i))*hicemelt) then
- wl4(i,is) = wl3(i) + wi1(i,is)/mhsep
- wi2(i,is) = 0.d0
- Tsoil3(i,is) = work + 0.01 + &
- & (Potphenergy-(wi1(i,is)*rosdry(i)*dzss(i)* &
- & (1 - por(i))*hicemelt))/(csoil(i)*rosoil(i)*dzss(i))
- else
- wl4(i,is) = wl3(i) + Potphenergy / &
- & (rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)/mhsep
- wi2(i,is) = wi1(i,is) - Potphenergy / &
- & (rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)
- Tsoil3(i,is) = work + 0.01
- endif
- else
- if (wl3(i) >= 0 .and. Tsoil2(i) < work - 0.01) then
- Potphenergy = - (Tsoil2(i) - work + 0.01) * &
- & csoil(i)*rosoil(i)*dzss(i)
- watavailtofreeze = (wl3(i) - UNFRWAT(Tsoil2(i),i)) * &
- & rosdry(i)*dzss(i)*(1 - por(i))*mhsep*hicemelt
- if (watavailtofreeze < 0) then
- !cc = csoil(i)*rosoil(i)*dzss(i)
- !cc1 = rosdry(i)*dzss(i)*(1-por(i))*Lwi
- !bb = (-Potphenergy-0.1*cc-wl3(i)*cc1)/cc
- !aa = cc1/cc
- !call phase_iter(aa,bb,i,Tsoil3(i))
- Tsoil3(i,is) = Tsoil2(i) + watavailtofreeze/(csoil(i)*rosoil(i)*dzss(i))
- wi2(i,is) = wi1(i,is) + (wl3(i) - UNFRWAT(Tsoil2(i),i))*mhsep
- wl4(i,is) = UNFRWAT(Tsoil2(i),i)
- cycle c2
- endif
- if (Potphenergy >= watavailtofreeze) then
- Tsoil3(i,is) = work - 0.01 - &
- & (Potphenergy - watavailtofreeze) / &
- & (csoil(i)*rosoil(i)*dzss(i))
- wi2(i,is) = wi1(i,is) + (wl3(i) - UNFRWAT(Tsoil2(i),i))*mhsep
- wl4(i,is) = UNFRWAT(Tsoil2(i),i)
- else
- wl4(i,is) = wl3(i) - Potphenergy/(rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)/mhsep
- wi2(i,is) = wi1(i,is) + Potphenergy/(rosdry(i)*dzss(i)*(1 - por(i))*hicemelt)
- Tsoil3(i,is) = work - 0.01
- endif
- else
- wl4(i,is) = wl3(i)
- wi2(i,is) = wi1(i,is)
- Tsoil3(i,is) = Tsoil2(i)
- endif
- endif
- enddo c2
-
- max1 = 0.
- do i = 1, ns
- if (Tsoil3(i,is) < &
- & MELTINGPOINT(Sals2(i,is)/wl4(i,is),pressoil(i),tricemethhydr%par) - 0.01 &
- &.and. abs(wl4(i,is) - UNFRWAT(Tsoil3(i,is),i)) > max1) then
- max1 = abs(wl4(i,is) - UNFRWAT(Tsoil3(i,is),i))
- endif
- enddo
+max1 = 0.
+do i = 1, ns
+ if (Tsoil3(i,is) < &
+ & MELTINGPOINT(Sals2(i,is)/wl4(i,is),pressoil(i),tricemethhydr%par) - 0.01 &
+ &.and. abs(wl4(i,is) - UNFRWAT(Tsoil3(i,is),i)) > max1) then
+ max1 = abs(wl4(i,is) - UNFRWAT(Tsoil3(i,is),i))
+ endif
+enddo
if (max1 > 0.01 .and. iter3 < 20) then
!wl3=(wl4+wl3)/2.
@@ -831,7 +831,7 @@ 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(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)
@@ -943,7 +943,8 @@ if (contr(2) == 1) then
& work, work, &
& wl4(1,is),wi2(1,is),wa,Sals2, &
& rootss,rosdry,por,veg,qsoilsurf(is)*por(1),TgrAnn, methgenmh, &
- & ddz,ddz05, wst, lammeth, 0.5*(zsoilcols(is) + zsoilcols(is+1)), &
+ & 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, &
@@ -979,7 +980,8 @@ if (contr(2) == 1) then
& work, work, &
& wl4(1,is),wi2(1,is),wa,Sals2, &
& rootss,rosdry,por,veg,fdiffbot(is),TgrAnn, methgenmh, &
- & ddz,ddz05, wst, lammeth, 0.5*(zsoilcols(is) + zsoilcols(is+1)), &
+ & 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, &
diff --git a/tools/plot_contour.py b/tools/plot_contour.py
index a9043df2d3565a5385f9c20da1d15ed47b69cf1b..c2f410ea6aa95405d7d7ac9121ba16cb80830de8 100644
--- a/tools/plot_contour.py
+++ b/tools/plot_contour.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
import numpy as nm
import matplotlib.pyplot as plt # for plotting
#import os
@@ -19,6 +20,15 @@ def runningMeanFast(x, N):
work1[i] = work1[len(x)-N]
return work1
+font = {'family': 'serif',
+ 'weight': 'normal'}
+rc('font', **font)
+rc('text', usetex=True)
+rc('text.latex',unicode=True)
+rc('text.latex',preamble='\usepackage[utf8]{inputenc}')
+rc('text.latex',preamble='\usepackage[russian]{babel}')
+#rc('text.latex',preamble='\usepackage{amsmath}')
+
varlist = \
['temp' ,\
@@ -127,8 +137,8 @@ if int(nplot) == varlist.index('temp') + 1: # Temperature
basename = 'water_temp 1 1'
fileout = 'temp_wat_cont'
nheader = 6
- labelx = 'Time' #, months
- labely = 'Depth, m'
+ labelx = u'Время' #'Time' #, months
+ labely = u'Глубина, м' #'Depth, m'
timescale = 1. #365/12 # 365/12*24 - number of hours in a month
time_unit = 1. #Time units in input data, expressed in days
ctime = 5
@@ -138,12 +148,12 @@ if int(nplot) == varlist.index('temp') + 1: # Temperature
#t0disp = 7.80 # Initial time for display
#t1disp = 7.95 # Final time for display
#nvardisp = nlevs_water
- title = "Temperature, $^\circ$C "
- #levmin = 6.
- #levmax = 27.
- nlevplot = 21 # Number of countour levels
+ title = u"Температура, $^\circ$C "#"Temperature, $^\circ$C "
+ levmin = 0.
+ levmax = 26.
+ nlevplot = 13 # Number of countour levels
pathsave = path
- nrunmean = 24
+ #nrunmean = 24
if int(nplot) == varlist.index('sal') + 1: # Salinity
basename = 'sal_water 1 1'
fileout = 'sal_wat_cont'
@@ -189,8 +199,8 @@ elif int(nplot) == varlist.index('ch4') + 1: # Methane
basename = 'methane_water 1 1'
fileout = 'meth_wat_cont'
nheader = 6 #Number of lines in the header
- labelx = 'Time, months'
- labely = 'Depth, m'
+ labelx = u'Время' #'Time, months'
+ labely = u'Глубина, м'#'Depth, m'
timescale = 365/12 # 365/12*24 - number of hours in a month
time_unit = 1. #Time units in input data, expressed in days
ctime = 5
@@ -199,19 +209,21 @@ elif int(nplot) == varlist.index('ch4') + 1: # Methane
#t0disp = 6 # Initial time for display
#t1disp = 11 # Final time for display
#nvardisp = nlevs_water
- title = "Methane, $mmol/m^3$"
- text = 'maxval'
+ title = u'CH$_4$, мкмоль/л' #"Methane, $mmol/m^3$"
+ #text = 'maxval'
#levmin = 0.
#levmax = 5000.
nlevplot = 20 # Number of countour levels
#stepticks = 60
pathsave = path
+ colormap = plt.cm.gnuplot
+ logscalebar = True
elif int(nplot) == varlist.index('o2') + 1: # Oxygen
basename = 'oxygen_water 1 1'
fileout = 'oxyg_wat_cont'
nheader = 6 #Number of lines in the header
- labelx = 'Time, months'
- labely = 'Depth, m'
+ labelx = u'Время' #'Time, months'
+ labely = u'Глубина, м' #'Depth, m'
timescale = 365/12 # 365/12*24 - number of hours in a month
time_unit = 1. #Time units in input data, expressed in days
ctime = 5
@@ -220,10 +232,10 @@ elif int(nplot) == varlist.index('o2') + 1: # Oxygen
#t0disp = 6 # Initial time for display
#t1disp = 10.3 # Final time for display
#nvardisp = nlevs_water
- title = "Oxygen, mg/l"
- #levmin = 0.
- #levmax = 10.
- nlevplot = 20 # Number of countour levels
+ title = u'Кислород, мг/л'#"Oxygen, mg/l"
+ levmin = 0.
+ levmax = 15.
+ nlevplot = 15 # Number of countour levels
pathsave = path
colormap = plt.cm.viridis
elif int(nplot) == varlist.index('co2') + 1: # Carbon dioxide
@@ -401,11 +413,12 @@ elif int(nplot) == varlist.index('ch4bub') + 1: # Methane bubble flux
ctime = 5
cstart = 7
time0 = 5 #Initial time
- t0disp = 5 # Initial time for display
- t1disp = 11 # Final time for display
+ #t0disp = 5 # Initial time for display
+ #t1disp = 11 # Final time for display
#nvardisp = nlevs_water
title = "Methane bubble flux, $mol/(m^2*s)$"
pathsave = path
+ logscalebar = True
elif int(nplot) == varlist.index('ri') + 1: # Richardson number
basename = 'Ri 1 1'
fileout = 'Ri_cont'
diff --git a/tools/plot_hourly.py b/tools/plot_hourly.py
index cf9c8e65bd3e992dc08cf396f3bc0809ce2d680f..2720b321a0c941786ced782557cc290747492bc5 100644
--- a/tools/plot_hourly.py
+++ b/tools/plot_hourly.py
@@ -27,7 +27,7 @@ plt.rc('text.latex',preamble='\usepackage[russian]{babel}')
#path = ['../results/windforc_base/hourly/', \
# '../results/windforc_kor/hourly/', \
# '../results/windforc_dynpgrad4/hourly/'] # path to files
-path = ['../results/Mojai2016-2017/hourly/','../results/Mojai2016-2017_damtribheat=0/hourly/']
+path = ['../results/Mojai2016-2017_damw5/hourly/']#,'../results/Mojai2016-2017_damtribheat=0/hourly/']
year = int(input("Enter year \n"))
nm0 = int(input("Enter month \n"))
@@ -35,6 +35,10 @@ nd0 = int(input("Enter first day \n"))
ndlast = int(input("Enter last day \n"))
nh0 = int(input("Enter first hour \n"))
nhlast = int(input("Enter last hour \n"))
+#figlabel = input("Enter figure label (put 0 for none) \n")
+
+figlabel = u'г)'
+
if nd0 == ndlast and nh0 == nhlast:
nhstep = 0
else:
@@ -85,6 +89,10 @@ ncols = []
nheader = []
scale = []
+varname_ = 'ch4'
+varnames_list = ['temp','o2','ch4','co2']
+nvarlist = varnames_list.index(varname_)
+
#varname = 'Absolute value of velocity'
#ylabel = 'Depth, m' #$\circ C$
#xlabel = varname + ', $m/s$'
@@ -93,50 +101,66 @@ scale = []
#ncols = 14 #Total number of columns
#nheader = 14 #Number of lines in the header
+if (nvarlist == varnames_list.index('temp')):
+ #varname.append(u'Temperature, model')
+ #ylabel = u'Depth, m' #$\circ C$
+ varname.append(u'T')
+ ylabel = u'Глубина, м'
+ xlabel = varname[len(varname)-1] + u',~$^{\circ}$C'
+ basefilename.append('Profiles')
+ ncol.append(2) #The column number to display, 2 for temperature
+ ncols.append(17) #Total number of columns
+ nheader.append(17) #Number of lines in the header
+ scale.append(1.)
+ logscale = False
+ xmin = 6.
+ xmax = 26.
+if (nvarlist == varnames_list.index('o2')):
+ #varname.append(u'O$_2$, model')
+ #ylabel = u'Depth, m' #$\circ C$
+ #xlabel = varname[len(varname)-1] + u',~mg/l'
+ varname.append(u'O$_2$')
+ ylabel = u'Глубина, м'
+ xlabel = varname[len(varname)-1] + u',~мг/л'
+ basefilename.append('Profiles')
+ ncol.append(15) #The column number to display, 2 for temperature
+ ncols.append(17) #Total number of columns
+ nheader.append(17) #Number of lines in the header
+ scale.append(1.)
+ logscale = False
+ xmin = 0.
+ xmax = 15.
+if (nvarlist == varnames_list.index('ch4')):
+ #varname.append(u'CH$_4$, model')
+ #ylabel = u'Depth, m' #$\circ C$
+ #xlabel = varname[len(varname)-1] + u',~mcmol/l'
+ varname.append(u'CH$_4$')
+ ylabel = u'Глубина, м'
+ xlabel = varname[len(varname)-1] + u',~ мкмоль/л'
+ basefilename.append('Profiles')
+ ncol.append(17) #The column number to display, 2 for temperature
+ ncols.append(17) #Total number of columns
+ nheader.append(17) #Number of lines in the header
+ scale.append(1.)
+ logscale = True
+ xmin = 8.E-2
+ xmax = 120.
+if (nvarlist == varnames_list.index('co2')):
+ #varname.append(u'CO$_2$, model')
+ #ylabel = u'Depth, m' #$\circ C$
+ #xlabel = varname[len(varname)-1] + u',~mcmol/l'
+ varname.append(u'CO$_2$')
+ ylabel = u'Глубина, м'
+ xlabel = varname[len(varname)-1] + u',~мкмоль/л'
+ basefilename.append('Profiles')
+ ncol.append(16) #The column number to display, 2 for temperature
+ ncols.append(17) #Total number of columns
+ nheader.append(17) #Number of lines in the header
+ scale.append(1.)
+ logscale = False
+ #xmin = -8.E-2
+ #xmax = 8.E-2
-#varname.append(u'Temperature, model')
-#ylabel = u'Depth, m' #$\circ C$
-#xlabel = varname[len(varname)-1] + u',~К'
-#basefilename.append('Profiles')
-#ncol.append(2) #The column number to display, 2 for temperature
-#ncols.append(17) #Total number of columns
-#nheader.append(17) #Number of lines in the header
-#scale.append(1.)
-##xmin = -8.E-2
-##xmax = 8.E-2
-
-varname.append(u'O$_2$, model')
-ylabel = u'Depth, m' #$\circ C$
-xlabel = varname[len(varname)-1] + u',~mg/l'
-basefilename.append('Profiles')
-ncol.append(15) #The column number to display, 2 for temperature
-ncols.append(17) #Total number of columns
-nheader.append(17) #Number of lines in the header
-scale.append(1.)
-#xmin = -8.E-2
-#xmax = 8.E-2
-
-#varname.append(u'CO$_2$, model')
-#ylabel = u'Depth, m' #$\circ C$
-#xlabel = varname[len(varname)-1] + u',~mcmol/l'
-#basefilename.append('Profiles')
-#ncol.append(16) #The column number to display, 2 for temperature
-#ncols.append(17) #Total number of columns
-#nheader.append(17) #Number of lines in the header
-#scale.append(1.)
-##xmin = -8.E-2
-##xmax = 8.E-2
-
-#varname.append(u'CH$_4$, model')
-#ylabel = u'Depth, m' #$\circ C$
-#xlabel = varname[len(varname)-1] + u',~mcmol/l'
-#basefilename.append('Profiles')
-#ncol.append(17) #The column number to display, 2 for temperature
-#ncols.append(17) #Total number of columns
-#nheader.append(17) #Number of lines in the header
-#scale.append(1.)
-##xmin = -8.E-2
-##xmax = 8.E-2
#varname.append(u'Компонента скорости по оси $x$')
#ylabel = u'Глубина, м' #$\circ C$
@@ -290,9 +314,11 @@ for np in range(npaths):
print k, line, ncol[k]
vals[j] = float(line[ncol[k]-1])
- legenditems.append(varname[k]+ ', ' +legendtime[i]+expnames[np])
+ #legenditems.append(varname[k] + u' (модель)') #+ ', ' +legendtime[i]+expnames[np])
+ legenditems.append(u'модель')
for m in range(len(vals)) : vals[m] = vals[m]*scale[k]
- plt.plot(vals,depths,ls=linestyles[np],color=linecolors[k+2],linewidth=3,\
+ plt.plot(vals,depths,ls=linestyles[np], \
+ color=linecolors[varnames_list.index(varname_)],linewidth=3,\
markersize=markersize)
#print(depths)
#print(vals)
@@ -307,10 +333,13 @@ for np in range(npaths):
if 'pathobs' in locals():
nprec = 3
nvars = 3
- legobs = [u'T obs',u'O$_2$ obs',u'CH$_4$ obs']
+ #legobs = [u'T obs',u'O$_2$ obs',u'CH$_4$ obs']
+ #legobs = [u'T (набл)',u'O$_2$ (набл)',u'CH$_4$ (набл)']
+ legobs = [u'набл',u'набл',u'набл']
nlocs = 6 #number of locations
- legobssuf = [u'section~V',u'section~IV',u'section~III',u'section~II',u'section~I',u'averaged']
- nvar = [2]
+ #legobssuf = [u'section~V',u'section~IV',u'section~III',u'section~II',u'section~I',u'averaged']
+ legobssuf = [u'ст.~V',u'ст.~IV',u'ст.~III',u'ст.~II',u'ст.~I',u'ср.']
+ nvar = [nvarlist+1]
nvarsd = len(nvar)
fobs = open(pathobs,'r')
#print(line)
@@ -333,11 +362,15 @@ if 'pathobs' in locals():
#print(z)
#print(varobs[:,2])
for i in range(nvarsd):
- plt.plot(varobs[:,i],z,'o', \
- color=linecolors[nvar[i]-1+k],linewidth=3,markersize=markersize)
+ #plt.plot(varobs[:,i],z,'o', \
+ #color=linecolors[nvar[i]-1+k],linewidth=3,markersize=markersize)
+ plt.plot(varobs[:,i],z,'o', \
+ color=linecolors[-1+k],linewidth=3,markersize=markersize)
+
plt.xlabel(xlabel,fontsize=fonts)
plt.ylabel(ylabel,fontsize=fonts)
+if (logscale) : plt.xscale('log')
plt.legend(legenditems,loc='best',fontsize=20)
plt.gca().tick_params(labelsize=20)
@@ -348,8 +381,12 @@ plt.gca().invert_yaxis()
if 'xmin' in locals(): plt.gca().set_xlim([xmin, xmax])
+print (figlabel)
+if figlabel != 0:
+ plt.figtext(0.035, 0.925, figlabel, \
+ backgroundcolor='w',size=20,fontdict=None)
-plt.savefig(path[np] + fileout + '_' + str(year) + str(nm).zfill(2) + str(ndlast).zfill(2) + str(nhlast).zfill(2) + '_hourly.png')
+#plt.savefig(path[np] + fileout + '_' + str(year) + str(nm).zfill(2) + str(ndlast).zfill(2) + str(nhlast).zfill(2) + '_hourly.png')
plt.show()
diff --git a/tools/plot_timser.py b/tools/plot_timser.py
index 6eaa00ce4c83515f0ab14bee2fef6efc6c59c088..a0180d2611cf1adcba50ab8515530fab6f26d238 100644
--- a/tools/plot_timser.py
+++ b/tools/plot_timser.py
@@ -256,7 +256,7 @@ def plot_timeser(plotid,nvar,legendcommon) :
#path_ = ['../results/IseoAugust2011_cdE-1/time_series/']#, \
# '../results/IseoAugust2011_cdE-1/time_series/']
#path_ = ['../results/Uvs1979-2017/time_series/']#, \
- path_ = ['../results/Mojai2016-2017_damw3/time_series/']
+ path_ = ['../results/Mojai2016-2017_damw5/time_series/']
#path_ = ['../results/IseoJuneOctober2011_cdE-1/time_series/']
#path_ = ['../results/IseoJuneOctober2011_botfric1backdiff1/time_series/']#, \
# '../results/IseoJuneOctober2011_botfric=1/time_series/'] #, \
@@ -280,7 +280,7 @@ def plot_timeser(plotid,nvar,legendcommon) :
nameexpobs = ['observ']
#nameexpobs = ['']
#nameexp_ = ['gradp+wind','gradp','wind']
- nameexp_ = ['model','corretahepi','tribheat0','netrad']
+ nameexp_ = ['','corretahepi','tribheat0','netrad']
#nameexp_ = ['dynpgrad4','dynpgrad0','cbotreduced']
#nameexp_ = [u'эксп. Като--Филлипса',u'эксп. Като--Филлипса + баротропные сейши', \
#u'эксп. Като--Филлипса + бароклинные сейши',u'эксп. Като--Филлипса + сила Кориолиса'] #
@@ -726,9 +726,9 @@ def plot_timeser(plotid,nvar,legendcommon) :
basename = 'layers 1 1'
fileout = 'layers_ser'
nheader = 19 #Number of lines in the header
- ncol = [15,16] #[6,15,16,17] #The number of columns in the file to display
- #legend = ['water thickness','ice thickness','snow thickness','deep ice thickness']
- legend = ['ice thickness','snow thickness']
+ ncol = [6,15,16,17] #The number of columns in the file to display
+ legend = ['water thickness','ice thickness','snow thickness','deep ice thickness']
+ #legend = ['ice thickness','snow thickness']
labelx = 'Time, months'
labely = 'Thickness, m'
ctime = 5
@@ -758,12 +758,14 @@ def plot_timeser(plotid,nvar,legendcommon) :
nheader = 33 #Number of lines in the header
ncol = [14,17,31] #The number of columns in the file to display
#ncol = [17] #The number of columns in the file to display
- legend = ['CH$_4$ ebullition flux', \
- 'CH$_4$ diffusion flux', \
- 'CH$_4$ flux through outlet']
- #legend = ['CH$_4$ diffusion flux']
- labelx = 'Time, months'
- labely = 'Methane flux, $\mu$mol/(m$^2$*s)'
+ #legend = ['CH$_4$ ebullition flux', \
+ # 'CH$_4$ diffusion flux', \
+ # 'CH$_4$ flux through outlet']
+ legend = [u'пузырьковый поток', \
+ u'диффузионный поток', \
+ u'поток через дамбу']
+ labelx = u'Время' #'Time, months'
+ labely = u'Поток метана, мкмоль/(м$^2$с)' #'Methane flux, $\mu$mol/(m$^2$*s)'
ctime = 5
timescale = 24 # 365/12*24 - number of hours in a month
time0 = 0. #month_doy(5,ny)/mdays #Initial time
@@ -771,7 +773,7 @@ def plot_timeser(plotid,nvar,legendcommon) :
path = path_
nameexp = nameexp_
npaths = len(path_)
- #logscale = 1
+ logscale = 1
elif int(nplot) == ch4flxMLmod_: # Methane fluxes
basename = 'methane_series 1 1'
fileout = 'methane_ML_budget'