diff --git a/setup/Mojai2016-2017_setup.dat b/setup/Mojai2016-2017_setup.dat
index 24454b663a4383880e56b73b6066fa2a4f09347f..3e3ec39620d9dae6ecb6befd102b0fe050acea75 100644
--- a/setup/Mojai2016-2017_setup.dat
+++ b/setup/Mojai2016-2017_setup.dat
@@ -146,6 +146,7 @@ sedim 0
salsoil 0
saltice 0
dyn_pgrad 0
+pgrad 0.
botfric 1
horvisc 0
backdiff 1
@@ -305,6 +306,13 @@ ngrid_out 26
23.
24.
#
+ngridice_out 5
+0.1
+0.2
+0.3
+0.4
+0.5
+#
ngridsoil_out 11
0.
1.
diff --git a/source/model/bathym_mod.f90 b/source/model/bathym_mod.f90
index 461f3194ba1f4fbd3a566cd49e89981bd606c99b..8079f2a80065e9491e77a6276ee2f7afb9173e65 100644
--- a/source/model/bathym_mod.f90
+++ b/source/model/bathym_mod.f90
@@ -30,7 +30,9 @@ use ARRAYS_GRID, only : &
& ddzi, nsoilcols, zsoilcols, &
& isoilcolc, ksoilcol, &
& z_full, z_half, &
-& dzeta_int, dzeta_05int, ddz05
+& z_full_ice, z_half_ice, &
+& dzeta_int, dzeta_05int, ddz05, &
+& dzetai_int, dzetai_05int
use ARRAYS_SOIL, only : dzs
@@ -80,7 +82,6 @@ logical :: flag
!real(kind=ireals), external :: ABSTR
-
do i = 1, M+1
z_full(i) = dzeta_int(i)*h1
preswat(i) = pressure + row0*g*z_full(i)
@@ -88,14 +89,8 @@ enddo
do i = 1, M
z_half(i) = dzeta_05int(i)*h1
enddo
-
-! Specification of profile output levels, if not given in setup file
-! In this case these levls coincide with numerical dzeta-levels
-! if (.not. allocated(zgrid_out)) then
-! ngrid_out%par = M+1
-! allocate (zgrid_out(1:ngrid_out%par))
-! endif
-! zgrid_out(:) = z_full(:)
+forall (i = 1:Mice+1) z_full_ice(i) = dzetai_int (i)*l1
+forall (i = 1:Mice ) z_half_ice(i) = dzetai_05int(i)*l1
! Constant cross-section
if (depth_area(1,2) < 0.) then ! No data on morphometry
diff --git a/source/model/driving_params_mod.f90 b/source/model/driving_params_mod.f90
index 4b7eb1dbd6b3f57fd5e6b2bba9de8af3a2d37e39..b12b43fa355f82c75672416b8671e46f3458034f 100644
--- a/source/model/driving_params_mod.f90
+++ b/source/model/driving_params_mod.f90
@@ -111,12 +111,9 @@ type(grarr1) :: rtemp
! Output grid
real(kind=ireals), allocatable :: zgrid_out(:)
+real(kind=ireals), allocatable :: zgridice_out(:)
real(kind=ireals), allocatable :: zgridsoil_out(:)
-
-! In perspective, it is expected, that area_lake be an array,
-! as the horizontal section of a lake varies with depth
-
integer(kind=iintegers) :: nunit = lake_subr_unit_min
type(intpar), target :: runmode
@@ -156,6 +153,7 @@ type(intpar), target :: everystep
type(intpar), target :: turb_out
type(intpar), target :: scale_output
type(intpar), target :: ngrid_out
+type(intpar), target :: ngridice_out
type(intpar), target :: ngridsoil_out
type(intpar), target :: zero_model
type(intpar), target :: nsoilcols_
@@ -401,13 +399,15 @@ if (firstcall) then
all_par_present = .false.
endif
if (.not.ngrid_out%ok) then
- write(*,*) 'The parameter ngrid_out is missing in setup file. &
- & Output grid for profiles will be created automatically'
- ! all_par_present = .false.
+ write(*,*) 'The parameter ngrid_out is missing in setup file.'
+ all_par_present = .false.
+ endif
+ if (.not.ngridice_out%ok) then
+ write(*,*) 'The parameter ngridice_out is missing in setup file.'
+ all_par_present = .false.
endif
if (.not.ngridsoil_out%ok) then
- write(*,*) 'The parameter ngridsoil_out is missing in &
- &setup file'
+ write(*,*) 'The parameter ngridsoil_out is missing in setup file.'
all_par_present = .false.
endif
if (.not.ok_assim) then
@@ -583,6 +583,18 @@ lineread: if (line(n1:n1) /= '#') then
zgrid_out(:) = work (:,1)
deallocate(work)
endif
+ CASE ('ngridice_out')
+ ngridice_out%par &
+ & = igetvarval(n1,n2,line,'Switch ice output regrid ')
+ allocate (zgridice_out (1:ngridice_out%par) )
+ ngridice_out%ok = .true.
+ if (ngridice_out%par > 0) then
+ ncol = 1
+ allocate (work(1:ngridice_out%par,ncol))
+ call READPROFILE(nunit,ngridice_out%par,ncol,work)
+ zgridice_out(:) = work (:,1)
+ deallocate(work)
+ endif
CASE ('ngridsoil_out')
ngridsoil_out%par &
& = igetvarval(n1,n2,line,'Switch soil out. regrid ')
diff --git a/source/model/init.f90 b/source/model/init.f90
index d22ca52213d673efb56ee5debf8f5dc946b161bb..019e588becdc3996febac8a7e8cda81f02b20321 100644
--- a/source/model/init.f90
+++ b/source/model/init.f90
@@ -33,6 +33,7 @@
allocate (dz_full(1:M), z_full(1:M+1), z_half (1:M), &
& zsoilcols(1:nsoilcols+1,1:nx,1:ny))
+ allocate (z_full_ice(1:Mice+1), z_half_ice(1:Mice))
allocate (isoilcolc(1:nsoilcols),ksoilcol(1:M+1))
allocate (ddz(1:M), ddz2(1:M-1), ddz05(0:M), ddz052(2:M-1), ddz054(2:M-1))
allocate (dzeta_int(1:M+1), dzeta_05int(1:M))
@@ -272,9 +273,10 @@
gs%nx => nx; gs%ny => ny
!Grid spacing group
- gsp%dz_full => dz_full; gsp%z_full => z_full; gsp%z_half => z_half;
+ gsp%dz_full => dz_full; gsp%z_full => z_full; gsp%z_half => z_half
+ gsp%z_full_ice => z_full_ice; gsp%z_half_ice => z_half_ice
gsp%zsoilcols => zsoilcols; gsp%isoilcolc => isoilcolc; gsp%ksoilcol => ksoilcol
- gsp%ddz => ddz; gsp%ddz2 => ddz2; gsp%ddz05 => ddz05;
+ gsp%ddz => ddz; gsp%ddz2 => ddz2; gsp%ddz05 => ddz05
gsp%ddz052 => ddz052; gsp%ddz054 => ddz054
gsp%ddzi => ddzi; gsp%ddzi05 => ddzi05
gsp%dzeta_int => dzeta_int; gsp%dzeta_05int => dzeta_05int
@@ -374,40 +376,45 @@
END SUBROUTINE GRID_CREATE
-SUBROUTINE LININTERPOL (z1,f1,N1,z2,f2,N2,flag)
+SUBROUTINE LININTERPOL (z1,f1,N1,z2,f2,N2,flag,trextr)
- use LAKE_DATATYPES, only : ireals, iintegers
+use LAKE_DATATYPES, only : ireals, iintegers
+use DRIVING_PARAMS, only : missing_value
implicit none
! Input variables
! Integers
integer(kind=iintegers), intent(in) :: N1, N2
-
! Reals
real(kind=ireals), intent(in) :: z1(1:N1), z2(1:N2)
-
real(kind=ireals), intent(in) :: f1(1:N1)
+logical, optional :: trextr
! Output variables
! Reals
real(kind=ireals), intent(out) :: f2(1:N2)
-
logical, intent(out) :: flag
! Local variables
integer(kind=iintegers) :: i, j ! Loop indices
+logical :: extrap
+
+extrap = .true.
+if (present(trextr)) extrap = trextr
flag = .true.
+if (extrap) then
+ forall(i = 1 : N2, z2(i) < z1(1)) f2(i) = f1(1) + &
+ & (f1(2) - f1(1)) / (z1(2)-z1(1)) * (z2(i) - z1(1))
+ forall(i = 1 : N2, z2(i) >= z1(N1)) f2(i) = f1(N1)
+else
+ forall(i = 1 : N2, z2(i) < z1(1 )) f2(i) = missing_value
+ forall(i = 1 : N2, z2(i) >= z1(N1)) f2(i) = missing_value
+endif
+
do i = 1, N2
- if (z2(i) < z1(1)) then
-! Linear extrapolation of the profile up to the upper surface
- f2(i) = f1(1) + (f1(2) - f1(1)) / (z1(2)-z1(1)) &
- & * (z2(i) - z1(1))
- elseif (z2(i) >= z1(N1)) then
-! The constant value beneath the lower point of the profile
- f2(i) = f1(N1)
- else
+ if (z2(i) >= z1(1) .and. z2(i) < z1(N1)) then
do j = 1, N1
if (z2(i) >= z1(j) .and. z2(i) < z1(j+1)) then
f2(i) = f1(j) + ( f1(j+1) - f1(j) ) * &
diff --git a/source/model/lake.f90 b/source/model/lake.f90
index 774f970c146320e46242645bc0571e31b292e6e1..b05781225a5bf0c08277248981f4af53442af03f 100644
--- a/source/model/lake.f90
+++ b/source/model/lake.f90
@@ -2559,6 +2559,15 @@ if (flag_print) then ! The output in ASCII files
& Tw1, z_full, M+1, &
& zgrid_out, ngrid_out%par, & ! ngrid_out%par
& outpath, 'water_temp', i, ndec, .false.) !, row, work1 - last optional argument!
+ ndec = -1
+ i = i + 2
+ call PROFILE_OUTPUT &
+ & (ix, iy, dnx, dny, &
+ & year, month, day, hour, &
+ & time, dt_out%par, &
+ & Ti1, z_full_ice, Mice+1, &
+ & zgridice_out, ngridice_out%par, & !ngridsoil_out%par, &
+ & outpath, 'ice_temp', i, ndec, .false.)
ndec = 4
i = i + 2
call PROFILE_OUTPUT &
diff --git a/source/model/lake_modules.f90 b/source/model/lake_modules.f90
index 40539137f433c6def6cbf7a5c5a834d30a0d6302..e992c8735f6f76960259b46e696ee0c2b0ca163c 100644
--- a/source/model/lake_modules.f90
+++ b/source/model/lake_modules.f90
@@ -43,7 +43,7 @@ real(kind=ireals), parameter :: min_ice_thick = 0.02
real(kind=ireals), parameter :: min_water_thick = 0.02
real(kind=ireals), parameter :: T_phase_threshold = 1.d-5
-real(kind=ireals), parameter :: minwind = 1.d-5
+real(kind=ireals), parameter :: minwind = 1.d-2
! ML --- total number of layers
! MS --- maximal number of layers in snow
@@ -447,6 +447,7 @@ type(gridsize_type) :: gs
! Grid spacing group
real(kind=ireals), allocatable, target :: dz_full(:), z_full(:), z_half(:), zsoilcols(:,:,:)
+real(kind=ireals), allocatable, target :: z_full_ice(:), z_half_ice(:)
real(kind=ireals), allocatable, target :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
real(kind=ireals), allocatable, target :: ddzi(:), ddzi05(:)
real(kind=ireals), allocatable, target :: dzeta_int(:), dzeta_05int(:)
@@ -455,6 +456,7 @@ integer(kind=iintegers), allocatable, target :: isoilcolc(:), ksoilcol(:)
type, public :: gridspacing_type
sequence
real(kind=ireals), pointer :: dz_full(:), z_full(:), z_half(:), zsoilcols(:,:,:)
+ real(kind=ireals), pointer :: z_full_ice(:), z_half_ice(:)
real(kind=ireals), pointer :: ddz(:), ddz2(:), ddz05(:), ddz052(:), ddz054(:)
real(kind=ireals), pointer :: ddzi(:), ddzi05(:)
real(kind=ireals), pointer :: dzeta_int(:), dzeta_05int(:)
diff --git a/source/model/out_mod.f90 b/source/model/out_mod.f90
index c1274be7a0d54b01fc4d08dbcd924d6029350fe4..a8258145690b06c19b9ddf8e6332111e4ccdf909 100644
--- a/source/model/out_mod.f90
+++ b/source/model/out_mod.f90
@@ -1367,9 +1367,9 @@ if (int(time/(dt_out*hour_sec))>count_out(ix,iy) .or. &
& abs(time - (count_out(ix,iy) + 1)*dt_out*hour_sec) < small_number) then
if (nout > 0) then
allocate (Profile_out(1:nout,1:nvars))
- call LININTERPOL (z_prof,Profile,nprof,z_out,Profile_out(1,1),nout,flag)
- if (nvars >= 2) call LININTERPOL (z_prof,Profile2,nprof,z_out,Profile_out(1,2),nout,flag)
- if (nvars == 3) call LININTERPOL (z_prof,Profile3,nprof,z_out,Profile_out(1,3),nout,flag)
+ call LININTERPOL (z_prof,Profile,nprof,z_out,Profile_out(1,1),nout,flag,.false.)
+ if (nvars >= 2) call LININTERPOL (z_prof,Profile2,nprof,z_out,Profile_out(1,2),nout,flag,.false.)
+ if (nvars == 3) call LININTERPOL (z_prof,Profile3,nprof,z_out,Profile_out(1,3),nout,flag,.false.)
k = nout
call DEFFORMAT
write (unit=n_unit(unitindex,ix,iy),fmt=formats(1)) year,month,day,hour, &