1.Ice temperature profile output added 2.In profiles output, missing values...

1.Ice temperature profile output added 2.In profiles output, missing values included instead of extrapolation outside the grid being interpolated

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.

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

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  ')

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 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.
 
-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)
+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) .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) ) * &

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 &

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(:)

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, &