Bugs fixed important for multi-lake simulations

source/model/init.f90

@@ -208,6 +208,7 @@ END MODULE GRID_OPERATIONS_LAKE
 
  allocate (dz_full(1:M), z_full(1:M+1), z_half (1:M), &
  & zsoilcols(1:nsoilcols+1,1:nx,1:ny))
+ zsoilcols(:,:,:) = missing_value
  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))

source/model/lake.f90

@@ -452,6 +452,7 @@ integer(kind=iintegers) :: layer_case
 integer(kind=iintegers) :: ndec
 integer(kind=iintegers) :: npoint_calibr(1:2) = 12
 integer(kind=iintegers) :: dnx, dny
+integer(kind=iintegers) :: iwork(1:2)
 
 !Logicals
 logical :: uniform_depth
@@ -648,7 +649,7 @@ dt_inv = 1.d0/dt
 dt_inv05 = 0.5d0*dt_inv
 
 ! Logical, indicating multiple soil columns configuration
-multsoil = (nsoilcols > 1 .and. depth_area(1,2) >= 0. .and. soilswitch%par == 1)
+multsoil = (nsoilcols > 1 .and. depth_area(1,2) > 0. .and. soilswitch%par == 1)
 
 if (init(ix,iy) == 0_iintegers) then
   ! Specification of initial profiles   
@@ -932,11 +933,12 @@ if (soilswitch%par == 1) then
 endif
 if (multsoil) then 
 ! Calculation of heat sources from heat fluxes from mutliple soil columns
-  call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+  iwork = (/1,0/)
+  call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols(1,ix,iy), &
                    & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                    & bathymwater,bathymice, &
                    & bathymdice,bathymsoil(1,ix,iy), &
-                   & soilflux(1,ix,iy),fdiffbot,(/1,0/))
+                   & soilflux(1,ix,iy),fdiffbot,iwork)
   call LATERHEAT(ix,iy,gs,ls, &
   & bathymwater,bathymice,bathymdice,bathymsoil, &
   & gsp,soilflux(1,ix,iy),.false.,lsh)
@@ -975,11 +977,12 @@ call OXYGEN_PRODCONS(gs, gsp, wst, bathymsoil, gas, area_int, area_half, ddz05,
 if (soilswitch%par == 1) then
   ! Calculation of methane in all soil columns, except for the lowest one
   if (multsoil) then
-    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+    iwork = (/1,0/)
+    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols(1,ix,iy), &
                      & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                      & bathymwater,bathymice, &
                      & bathymdice,bathymsoil(1,ix,iy), &
-                     & soilflux(1,ix,iy), fdiffbot, (/0,1/))
+                     & soilflux(1,ix,iy), fdiffbot, iwork)
     methsoilflux(1:nsoilcols-1) = fdiffbot(1:nsoilcols-1)
   endif
   call BUBBLE_BLOCK
@@ -1305,11 +1308,12 @@ if (flag_snow == 1) then
   endif
   ! Calculation of heat sources from heat fluxes from mutliple soil columns
   if (multsoil) then
-    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+    iwork = (/1,0/)
+    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols(1,ix,iy), &
                      & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                      & bathymwater,bathymice, &
                      & bathymdice,bathymsoil(1,ix,iy), &
-                     & soilflux(1,ix,iy), fdiffbot, (/1,0/))
+                     & soilflux(1,ix,iy), fdiffbot, iwork)
     call LATERHEAT(ix,iy,gs,ls, &
     & bathymwater,bathymice,bathymdice,bathymsoil, &
     & gsp,soilflux(1,ix,iy),.false.,lsh)
@@ -1334,11 +1338,12 @@ if (flag_snow == 1) then
   if (soilswitch%par == 1) then
     ! Calculation of methane in all soil columns, except for the lowest one
     if (multsoil) then
-      call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+      iwork = (/1,0/)
+      call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols(1,ix,iy), &
                        & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                        & bathymwater,bathymice, &
                        & bathymdice,bathymsoil(1,ix,iy), &
-                       & soilflux(1,ix,iy), fdiffbot, (/0,1/))
+                       & soilflux(1,ix,iy), fdiffbot, iwork)
       methsoilflux(1:nsoilcols-1) = fdiffbot(1:nsoilcols-1)
     endif
   endif
@@ -1561,11 +1566,12 @@ else
   endif
   ! Calculation of heat sources from heat fluxes from mutliple soil columns
   if (multsoil) then 
-    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+    iwork = (/1,0/)
+    call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols(1,ix,iy), &
                      & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                      & bathymwater,bathymice, &
                      & bathymdice,bathymsoil(1,ix,iy), &
-                     & soilflux(1,ix,iy), fdiffbot, (/1,0/))
+                     & soilflux(1,ix,iy), fdiffbot, iwork)
     call LATERHEAT(ix,iy,gs,ls, &
     & bathymwater,bathymice,bathymdice,bathymsoil, &
     & gsp,soilflux(1,ix,iy),.false.,lsh)
@@ -1598,11 +1604,12 @@ else
   if (soilswitch%par == 1) then
   ! Calculation of temperature in all soil columns, except for the lowest one
     if (multsoil) then
-      call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols, &
+      iwork = (/1,0/)
+      call SOILCOLSTEMP(nveclen,gs,gsp,dt,ls,ftot,ch4_pres_atm,ddz,ddzi,zsoilcols(1,ix,iy), &
                        & wst,RadWater%integr,a,b,c,d,add_to_winter, &
                        & bathymwater,bathymice, &
                        & bathymdice,bathymsoil(1,ix,iy), &
-                       & soilflux(1,ix,iy), fdiffbot, (/0,1/))
+                       & soilflux(1,ix,iy), fdiffbot, iwork)
       methsoilflux(1:nsoilcols-1) = fdiffbot(1:nsoilcols-1)
     endif
     ! Calling bubble model