diff --git a/source/Makefile b/source/Makefile
index 097572a9e06b5940d96a88f2103f10e70532f005..857cdba7516c1a5a66831148aa725c3cb8f3b57d 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -45,11 +45,14 @@
  $(objfiles_path)lake_modules.o \
  $(objfiles_path)dzeta_mod.o \
  $(objfiles_path)numerics_mod.o \
+ $(objfiles_path)time.o \
  $(objfiles_path)phys_func.o \
+ $(objfiles_path)init.o \
  $(objfiles_path)out_mod.o \
  $(objfiles_path)surf_scheme1_mod.o \
  $(objfiles_path)surf_scheme2_mod.o \
  $(objfiles_path)surf_scheme_inm_mod.o \
+ $(objfiles_path)surf_scheme3.o \
  $(objfiles_path)heatbalsurf_mod.o \
  $(objfiles_path)t_solver_mod.o \
  $(objfiles_path)driver_datatypes_mod.o \
@@ -75,14 +78,13 @@
  \
  $(objfiles_path)massflux_convection_v10.o \
  $(objfiles_path)convectpar.o \
+ $(objfiles_path)snowcalc.o \
+ $(objfiles_path)snowtemp.o \
+ $(objfiles_path)t_0dim.o \
  $(objfiles_path)lake.o \
- $(objfiles_path)init.o \
  $(objfiles_path)netcdf_inout.o \
- $(objfiles_path)time.o \
  $(objfiles_path)lakinteract.o \
  $(objfiles_path)driver.o \
- $(objfiles_path)snowtemp.o \
- $(objfiles_path)snowcalc.o \
  $(objfiles_path)flake_albedo_ref.o \
  $(objfiles_path)flake_derivedtypes.o \
  $(objfiles_path)flake_paramoptic_ref.o \
@@ -92,8 +94,6 @@
  $(objfiles_path)src_flake_interface_1D.o \
  $(objfiles_path)postprocessing.o \
  $(objfiles_path)postprocess_surf.o \
- $(objfiles_path)surf_scheme3.o \
- $(objfiles_path)t_0dim.o \
  \
  $(objfiles_path)sgetrf.o \
  $(objfiles_path)sgetrs.o \
diff --git a/source/allocarrays_mod.mod b/source/allocarrays_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..cbe3cd5df5bc381452a54dedcbef8823434bce32
Binary files /dev/null and b/source/allocarrays_mod.mod differ
diff --git a/source/arrays.mod b/source/arrays.mod
new file mode 100644
index 0000000000000000000000000000000000000000..92377da5dbecd0c623f0c8f57fc12156bbf565b8
Binary files /dev/null and b/source/arrays.mod differ
diff --git a/source/arrays_bathym.mod b/source/arrays_bathym.mod
new file mode 100644
index 0000000000000000000000000000000000000000..7bc6f516a84ab2a396d44df2bbb3a26bb4135099
Binary files /dev/null and b/source/arrays_bathym.mod differ
diff --git a/source/arrays_grid.mod b/source/arrays_grid.mod
new file mode 100644
index 0000000000000000000000000000000000000000..68ed382b4530da1f8ce6f03b107a2bdfca3760c9
Binary files /dev/null and b/source/arrays_grid.mod differ
diff --git a/source/arrays_methane.mod b/source/arrays_methane.mod
new file mode 100644
index 0000000000000000000000000000000000000000..543c8b3f265f9c6c1e1e80d1ec77f193644b8df2
Binary files /dev/null and b/source/arrays_methane.mod differ
diff --git a/source/arrays_oxygen.mod b/source/arrays_oxygen.mod
new file mode 100644
index 0000000000000000000000000000000000000000..eff7fd534bce8caadffefb0866c802b38a1e9db9
Binary files /dev/null and b/source/arrays_oxygen.mod differ
diff --git a/source/arrays_soil.mod b/source/arrays_soil.mod
new file mode 100644
index 0000000000000000000000000000000000000000..bc58d3d0b57e3612b044aeeccf91a8a80c25f215
Binary files /dev/null and b/source/arrays_soil.mod differ
diff --git a/source/arrays_turb.mod b/source/arrays_turb.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1332fa9140794ba8ab43faf37da69429993e4bf0
Binary files /dev/null and b/source/arrays_turb.mod differ
diff --git a/source/arrays_waterstate.mod b/source/arrays_waterstate.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d59f2a8d8ff91b7c162d852100c1b5076f26f5a6
Binary files /dev/null and b/source/arrays_waterstate.mod differ
diff --git a/source/atmos.mod b/source/atmos.mod
new file mode 100644
index 0000000000000000000000000000000000000000..4af56232526954edc19c1e1019fdecdc90fadd8b
Binary files /dev/null and b/source/atmos.mod differ
diff --git a/source/bathymsubr.mod b/source/bathymsubr.mod
new file mode 100644
index 0000000000000000000000000000000000000000..7efd873238696ee65d2dca983838cad32db1832f
Binary files /dev/null and b/source/bathymsubr.mod differ
diff --git a/source/bl_mod_lake.mod b/source/bl_mod_lake.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1e7ed2cffff82fc95b4ea469ae98d232ba11b381
Binary files /dev/null and b/source/bl_mod_lake.mod differ
diff --git a/source/bubble_lake_mod.mod b/source/bubble_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5918bca7ad219c3751ad68e0eb789c6db69d895e
Binary files /dev/null and b/source/bubble_lake_mod.mod differ
diff --git a/source/carbon_dioxide_lake_mod.mod b/source/carbon_dioxide_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5717bf65c51d680a64414edac0b8c2dec30b6528
Binary files /dev/null and b/source/carbon_dioxide_lake_mod.mod differ
diff --git a/source/comparams.mod b/source/comparams.mod
new file mode 100644
index 0000000000000000000000000000000000000000..670c603c872bbecc5150bd962b26939a742b361a
Binary files /dev/null and b/source/comparams.mod differ
diff --git a/source/control_point_lake_mod.mod b/source/control_point_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5565a2f05b82bb021b095ce0933c6200df30db85
Binary files /dev/null and b/source/control_point_lake_mod.mod differ
diff --git a/source/convectpar_lake_mod.mod b/source/convectpar_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..e5234e216f30b3c8fe5c8519d13e314416abe4c9
Binary files /dev/null and b/source/convectpar_lake_mod.mod differ
diff --git a/source/data_parameters.mod b/source/data_parameters.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a49acec8c35713ee2d55ed01044675b5644cd617
Binary files /dev/null and b/source/data_parameters.mod differ
diff --git a/source/diffusion_lake_mod.mod b/source/diffusion_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3afc04af2fc2489d1c86c4a28f4327efa7b7384a
Binary files /dev/null and b/source/diffusion_lake_mod.mod differ
diff --git a/source/driver/driver.f90 b/source/driver/driver.f90
index b004d351b69f7f2e8d60d001dbccef53621a5670..1e9ba302ee4caae2400aa7d22b762e7384a70ae8 100644
--- a/source/driver/driver.f90
+++ b/source/driver/driver.f90
@@ -29,6 +29,9 @@
   use PHYS_FUNC, only : &
   & NETLWRAD
 
+  use LAKE_INTERACTION, only : LAKETRIBINT, OUTFLOW_DISCHARGE
+  use TIME_LAKE_MOD, only : JULIAN_DATE
+
 #ifdef MPI
       use MPI
 #endif
@@ -246,8 +249,6 @@
   character(len=100) :: format_screen ! Work variable
   character(len=60) :: work_character
 
-  real(kind=ireals), external :: outflow_discharge
-
   data time    /0./
   data i_year  /0/
   data nstep /0/
@@ -745,7 +746,7 @@
      iycyc : do iy = 1, ny
 
        tribinfl = trib_inflow_2d(ix,iy) - &
-       & OUTFLOW_DISCHARGE(ndatamax,effl_outflow_2d,h_2d) + &
+       & OUTFLOW_DISCHARGE(ndatamax,effl_outflow_2d(1,ix,iy),h_2d(ix,iy)) + &
        & dhdt(iy) ! Adding interaction of two lakes
 
        call LAKE(ta(ix),qa(ix),pa(ix),ua(ix),va(ix),atm_rad(ix),Rad_bal(ix), &
diff --git a/source/driver/lakinteract.f90 b/source/driver/lakinteract.f90
index 6742509251fa5c6801a926db74f7228a54793a6c..1a29857a8689d6c715173a46885d1c3634905db8 100644
--- a/source/driver/lakinteract.f90
+++ b/source/driver/lakinteract.f90
@@ -1,3 +1,5 @@
+MODULE LAKE_INTERACTION
+contains
 SUBROUTINE LAKETRIBINT(larea,h,hbot,hbottrib,width,length,nlakes, &
 & dhdt)
 
@@ -34,8 +36,6 @@ real(kind=ireals) :: sgn(1:nlakes)
 integer(kind=iintegers) :: botfric = 1 ! Switch between Shezy (1) and Mannings (2) bottom friction
 integer(kind=iintegers) :: i ! Loop index
 
-real(kind=ireals), external :: HYDRAULIC_RADIUS, COEF_SHEZY, SCROSS_AREA
-
 ! It is assumed that the channel flow is directed from the edge (i.e. lake)
 ! with higher channel bottom and has uniform depth and width
 select case (hbot(1) + hbottrib(1) > hbot(2) + hbottrib(2))
@@ -178,3 +178,4 @@ HYDRAULIC_RADIUS = max(HYDRAULIC_RADIUS,small_number)
 
 END FUNCTION HYDRAULIC_RADIUS
 
+END MODULE LAKE_INTERACTION
diff --git a/source/driver/netcdf_inout.f90 b/source/driver/netcdf_inout.f90
index b5e9954ca88d6180ed21c641963c6860e570cf78..69713629a60bfb34b834b29dc639fcda00eb9382 100644
--- a/source/driver/netcdf_inout.f90
+++ b/source/driver/netcdf_inout.f90
@@ -1,4 +1,6 @@
 #ifdef NETCDF_LIB
+MODULE NETCDF_LAKE_DRIVER
+contains
 SUBROUTINE NETCDF_READ_FORCING &
  & (filename1, rewind_netcdf, npoints, &
  & lat, lon, &
@@ -366,5 +368,6 @@ ENDIF
 END SUBROUTINE NCERROR
 
 END SUBROUTINE NETCDF_WRITE_SURF
+END MODULE NETCDF_LAKE_DRIVER
 #endif
 
diff --git a/source/driver/postprocess_surf.f90 b/source/driver/postprocess_surf.f90
index 2077c481dbcc32deb999cd9781d9dc1d93bc85dd..821fd6119f44a50f5158c4aaa687db7f23e0a866 100644
--- a/source/driver/postprocess_surf.f90
+++ b/source/driver/postprocess_surf.f90
@@ -1,4 +1,6 @@
 #ifdef NETCDF_LIB
+MODULE NETCDF_POSTPROCESS_LAKE
+contains
 SUBROUTINE NETCDF_POSTPROCESS_SURF &
 & (filename1,moving_average_window,mean_cycle_period)
 
@@ -355,7 +357,6 @@ call NCERROR(NF90_CLOSE(ncid))
 
 END SUBROUTINE NETCDF_WRITE_MEAN_CYCLE
 
-
 END SUBROUTINE NETCDF_POSTPROCESS_SURF
+END MODULE NETCDF_POSTPROCESS_LAKE
 #endif
-
diff --git a/source/driver/postprocessing.f90 b/source/driver/postprocessing.f90
index 82bc350f3910a3df94c63134e002d9b8b5d1f155..54239523fcffb4ab979da7540b83ba15829c3819 100644
--- a/source/driver/postprocessing.f90
+++ b/source/driver/postprocessing.f90
@@ -1,5 +1,6 @@
+MODULE POSTPROCESSING_LAKE
+contains
 ! Postprocessing procedures
-
 SUBROUTINE SIMPLE_MOVING_AVERAGE(series_input, series_averaged, &
 & series_length, average_window)
 
@@ -92,4 +93,4 @@ mean_cycle = mean_cycle/real(number_periods)
 
 END SUBROUTINE MEAN_CYCLE_CALC
 
- 
+END MODULE POSTPROCESSING_LAKE
diff --git a/source/driver_datatypes.mod b/source/driver_datatypes.mod
new file mode 100644
index 0000000000000000000000000000000000000000..befc27a796175f49dba9054422edc1846a169b35
Binary files /dev/null and b/source/driver_datatypes.mod differ
diff --git a/source/driver_interfaces_mod.mod b/source/driver_interfaces_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..2be2776e4a7c1b7920602497ddcd5dd024af4d0c
Binary files /dev/null and b/source/driver_interfaces_mod.mod differ
diff --git a/source/driver_parameters.mod b/source/driver_parameters.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d303f709f241b21cc8d2672eb5652d1c53572bed
Binary files /dev/null and b/source/driver_parameters.mod differ
diff --git a/source/driving_params.mod b/source/driving_params.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a1039719b9b354ff692bcc14c1efb465bce69c7a
Binary files /dev/null and b/source/driving_params.mod differ
diff --git a/source/dzeta_mod.mod b/source/dzeta_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5e17b5f080d7300b156a434fd793988947c58e8c
Binary files /dev/null and b/source/dzeta_mod.mod differ
diff --git a/source/evolution_variables.mod b/source/evolution_variables.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d36eaedcbec3472e15e2cb9e1c86de3dd4733c6f
Binary files /dev/null and b/source/evolution_variables.mod differ
diff --git a/source/flake.mod b/source/flake.mod
new file mode 100644
index 0000000000000000000000000000000000000000..cb58ca14c59c637ad3edf32d2944ab8ecfb2d2b5
Binary files /dev/null and b/source/flake.mod differ
diff --git a/source/flake_albedo_ref.mod b/source/flake_albedo_ref.mod
new file mode 100644
index 0000000000000000000000000000000000000000..7b9883fa5deca6df56a455de3c9f6cbf74c514fc
Binary files /dev/null and b/source/flake_albedo_ref.mod differ
diff --git a/source/flake_configure.mod b/source/flake_configure.mod
new file mode 100644
index 0000000000000000000000000000000000000000..8afe5bf2351f4f647b2c4f8aacdacfdd872a7633
Binary files /dev/null and b/source/flake_configure.mod differ
diff --git a/source/flake_derivedtypes.mod b/source/flake_derivedtypes.mod
new file mode 100644
index 0000000000000000000000000000000000000000..8ce5ca1a7b09ebb03c7481d8803ed828ab060f0e
Binary files /dev/null and b/source/flake_derivedtypes.mod differ
diff --git a/source/flake_parameters.mod b/source/flake_parameters.mod
new file mode 100644
index 0000000000000000000000000000000000000000..df0867ea7914b431b2d90df1c31010fdd4bb3745
Binary files /dev/null and b/source/flake_parameters.mod differ
diff --git a/source/flake_paramoptic_ref.mod b/source/flake_paramoptic_ref.mod
new file mode 100644
index 0000000000000000000000000000000000000000..ed4444eb2e63a77627d0193697c25291fa143ad7
Binary files /dev/null and b/source/flake_paramoptic_ref.mod differ
diff --git a/source/grid_operations_lake.mod b/source/grid_operations_lake.mod
new file mode 100644
index 0000000000000000000000000000000000000000..64e58e25f90e3ebf12a74926dde35693cdd310ec
Binary files /dev/null and b/source/grid_operations_lake.mod differ
diff --git a/source/heatbalance_mod.mod b/source/heatbalance_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..b959c2edadf79a93747d6a11756ea140e5df5395
Binary files /dev/null and b/source/heatbalance_mod.mod differ
diff --git a/source/init_var_mod.mod b/source/init_var_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..9cd3366e994482c74e839c95db4b0c8f43441fcf
Binary files /dev/null and b/source/init_var_mod.mod differ
diff --git a/source/inout.mod b/source/inout.mod
new file mode 100644
index 0000000000000000000000000000000000000000..cf1b7ea0854156aa7358fe0220167bf545e852a8
Binary files /dev/null and b/source/inout.mod differ
diff --git a/source/inout_driver_parameters.mod b/source/inout_driver_parameters.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5f3f3657acced031ef350451dd21bf02ff404e08
Binary files /dev/null and b/source/inout_driver_parameters.mod differ
diff --git a/source/inout_parameters.mod b/source/inout_parameters.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d2c2c32f5fd67fbff191e081d143c1969cf602fd
Binary files /dev/null and b/source/inout_parameters.mod differ
diff --git a/source/lake_datatypes.mod b/source/lake_datatypes.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3bba1749f937f1f7910345b9cc11c1f8a4610e33
Binary files /dev/null and b/source/lake_datatypes.mod differ
diff --git a/source/lake_interaction.mod b/source/lake_interaction.mod
new file mode 100644
index 0000000000000000000000000000000000000000..4590d6a44880254d554e05476ca6c0d63ddff5b8
Binary files /dev/null and b/source/lake_interaction.mod differ
diff --git a/source/lake_mod.mod b/source/lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..06ec6ba90f14d9f65eda118c56d6f2c4492f9fb1
Binary files /dev/null and b/source/lake_mod.mod differ
diff --git a/source/meth_oxyg_constants.mod b/source/meth_oxyg_constants.mod
new file mode 100644
index 0000000000000000000000000000000000000000..4ed57f390d0e6b53109d4a9b8c5a3940dd8e4279
Binary files /dev/null and b/source/meth_oxyg_constants.mod differ
diff --git a/source/methane_lake_mod.mod b/source/methane_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..2a1567c5c8197a5b3bbf5cb913bad4f006ed2114
Binary files /dev/null and b/source/methane_lake_mod.mod differ
diff --git a/source/methane_simple.mod b/source/methane_simple.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f8cd695bd9159b772d4f0278b4fdca60eeb5a356
Binary files /dev/null and b/source/methane_simple.mod differ
diff --git a/source/model/bathym_mod.f90 b/source/model/bathym_mod.f90
index 5ca2124dd07d8dad8a3effa25498cd1c0e060a95..6b4070320f7f588861e0984ddc934e98de45c1cd 100644
--- a/source/model/bathym_mod.f90
+++ b/source/model/bathym_mod.f90
@@ -1,5 +1,7 @@
 MODULE BATHYMSUBR
 
+use GRID_OPERATIONS_LAKE, only : LININTERPOL, PIECEWISECONSTINT
+
 contains
 !> Subroutine BATHYMETRY calculates lake bathymetry characteristics
 SUBROUTINE BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
diff --git a/source/model/bubble_mod.f90 b/source/model/bubble_mod.f90
index 3f8558e0334b996b6559ede74b15fe20cb0a56ea..6b2b1aa6fa7cf6718a0af8c51fbb821c468e0822 100644
--- a/source/model/bubble_mod.f90
+++ b/source/model/bubble_mod.f90
@@ -1,4 +1,4 @@
-MODULE BUBBLE_MOD
+MODULE BUBBLE_LAKE_MOD
 
 use LAKE_DATATYPES, only : ireals, iintegers
 
@@ -411,4 +411,4 @@ enddo
 END SUBROUTINE BUBBLEFLUXAVER
 
 
-END MODULE BUBBLE_MOD
+END MODULE BUBBLE_LAKE_MOD
diff --git a/source/model/carbon_dioxide.f90 b/source/model/carbon_dioxide.f90
index c22f3f8a94b36ec3d83ddc5403775ca56117a65b..58733c2a6218fdccd3acf8c4dc4756e6ed7758b0 100644
--- a/source/model/carbon_dioxide.f90
+++ b/source/model/carbon_dioxide.f90
@@ -1,4 +1,4 @@
-MODULE CARBON_DIOXIDE_MOD
+MODULE CARBON_DIOXIDE_LAKE_MOD
 
 contains
 !> Subroutine CARBON_DIOXIDE calculates :
@@ -274,4 +274,4 @@ deallocate(conc)
 
 END SUBROUTINE CARBON_DIOXIDE
 
-END MODULE CARBON_DIOXIDE_MOD
+END MODULE CARBON_DIOXIDE_LAKE_MOD
diff --git a/source/model/control_point_mod.f90 b/source/model/control_point_mod.f90
index 4bc128f326aeef922d6f82785cf37e9ec749db2c..37f42f6c033045c1f895bb570d136ba2c69152eb 100644
--- a/source/model/control_point_mod.f90
+++ b/source/model/control_point_mod.f90
@@ -1,4 +1,4 @@
-MODULE CONTROL_POINT
+MODULE CONTROL_POINT_LAKE_MOD
 
 use LAKE_DATATYPES, only : ireals, iintegers
 use COMPARAMS, only : parallel_params
@@ -756,4 +756,4 @@ endif
 END SUBROUTINE SHIFTVI
 
 
-END MODULE CONTROL_POINT
+END MODULE CONTROL_POINT_LAKE_MOD
diff --git a/source/model/convectpar.f90 b/source/model/convectpar.f90
index 2f571e5c8ceca55408cbbd171fbeb3790ff73bb5..16e4321235322ff18c2edcc7fa3a904064beaf60 100644
--- a/source/model/convectpar.f90
+++ b/source/model/convectpar.f90
@@ -1,5 +1,4 @@
-MODULE CONVECTPAR
-
+MODULE CONVECTPAR_LAKE_MOD
 
 !INTERFACE
 !SUBROUTINE WaterMix(M,flag, T1, T2, ro1)
@@ -248,4 +247,4 @@ implicit none
   !print *, 'Nr, mix ', Nr, mix
 END SUBROUTINE UnDenGrMix
 
-END MODULE CONVECTPAR
+END MODULE CONVECTPAR_LAKE_MOD
diff --git a/source/model/diffusion_mod.f90 b/source/model/diffusion_mod.f90
index 754f2765aed13a64e01cea5546d90257c53c19a7..2a38f8be78b2c0b3f62d9a78167e470f8012acf3 100644
--- a/source/model/diffusion_mod.f90
+++ b/source/model/diffusion_mod.f90
@@ -1,4 +1,4 @@
-MODULE DIFFUSION_MOD
+MODULE DIFFUSION_LAKE_MOD
 
 use LAKE_DATATYPES, only : ireals, iintegers
 
@@ -26,7 +26,7 @@ real(kind=ireals)      , intent(in) :: dh0         !> Increment of the layer at
 real(kind=ireals)      , intent(in) :: bcs(1:2)    !> Value of variable or flux at boundaries
 real(kind=ireals)      , intent(in) :: diff(1:N-1) !> Diffusion coefficient
 real(kind=ireals)      , intent(in) :: ddz(1:N-1)  !> Grid spacing, n/d
-real(kind=ireals)      , intent(in) :: ddz05(0:N-1) !> Distance between cells centers, n/d
+real(kind=ireals)      , intent(in) :: ddz05(0:N-1) !> Distance between cells centers, n/dKE
 real(kind=ireals)      , intent(in) :: dzeta_int(1:N) !> Non-dimensional coordinate of cell interfaces
 real(kind=ireals)      , intent(in) :: lsh(1:N)    !> Exchange with sediments
 type(bathym)           , intent(in) :: bath(1:N)   !> Bathymetry
@@ -97,4 +97,4 @@ call PROGONKA (a,b,c,f,var2,1_iintegers,N)
 deallocate (a,b,c,f,y)
 END SUBROUTINE DIFFUSION
 
-END MODULE DIFFUSION_MOD
+END MODULE DIFFUSION_LAKE_MOD
diff --git a/source/model/heatbalsurf_mod.f90 b/source/model/heatbalsurf_mod.f90
index 6cddd7fa3692ef2a0deffb1ee93de5d80cb80eaf..26ad8e54e7627d60c7f51b2fdb04b39a0e19a6a7 100644
--- a/source/model/heatbalsurf_mod.f90
+++ b/source/model/heatbalsurf_mod.f90
@@ -217,7 +217,9 @@ FUNCTION HEATBALANCE(Tsurf,surftyp,RadWater,RadIce,fetch,dt)
  use ARRAYS_BATHYM, only : h1,l1
 
  use SURF_SCHEME1_MOD, only : DRAGVL
+ use SURF_SCHEME3_LAKE_MOD, only : SURF_SCHEME3
  use SURF_SCHEME_INM, only : DRAG3_LAKE
+ use RICHPBL_LAKE_MOD, only: RichPBL
 
  use SFCFLX, only : SFCFLX_MOMSENLAT
 
diff --git a/source/model/init.f90 b/source/model/init.f90
index 0d9888c40ef009dfcbb2bb0fa9ed30b7c912fa49..b5fc89c141f6210f5bf7350d546030a5658de221 100644
--- a/source/model/init.f90
+++ b/source/model/init.f90
@@ -1,3 +1,177 @@
+ MODULE GRID_OPERATIONS_LAKE
+
+ contains
+ SUBROUTINE GRID_CREATE
+
+ use LAKE_DATATYPES, only : ireals, iintegers
+
+ use ARRAYS_GRID
+ use DRIVING_PARAMS
+ use DZETA_MOD, only : DZETA, DZETAI
+
+ implicit none
+
+ integer(kind=iintegers) :: i
+ real(kind=ireals) :: dzeta1, dzeta2
+
+!The grid according to Burchard, 2002, pp. 97-98
+ do i = 1, M
+   dzeta1 = &
+   & (tanh(d_surf%par)+tanh((d_surf%par+d_bot%par)*real(i-1)/real(M)-d_surf%par))/ &
+   & (tanh(d_surf%par)+tanh(d_bot%par))
+   dzeta2 = &
+   & (tanh(d_surf%par)+tanh((d_surf%par+d_bot%par)*real(i)/real(M)-d_surf%par))/ &
+   & (tanh(d_surf%par)+tanh(d_bot%par))
+   ddz(i) = dzeta2 - dzeta1
+ enddo
+ 
+ ddz05(0) = 0.5d0*ddz(1)
+ do i = 1, M-1
+   ddz2(i) = ddz(i) + ddz(i+1)
+   ddz05(i) = 0.5d0*(ddz(i) + ddz(i+1))   
+ enddo
+ ddz05(M) = 0.5d0*ddz(M)
+ 
+ do i = 2, M-1
+   ddz052(i) = 0.5d0*ddz(i-1) + ddz(i) + 0.5d0*ddz(i+1)
+   ddz054(i) = 2.d0*ddz052(i)
+ enddo
+ 
+ do i = 1, M+1 
+   dzeta_int(i) = DZETA(real(i))
+ enddo
+ 
+ do i = 1, M
+   dzeta_05int(i) = DZETA(real(i)+0.5)
+ enddo
+ 
+ ddzi(:) = 1.d0/float(Mice)
+
+ ddzi05(0) = 0.5d0*ddzi(1)
+ do i = 1, Mice-1
+   ddzi05(i) = 0.5d0*(ddzi(i) + ddzi(i+1))
+ enddo
+ ddzi05(Mice) = 0.5d0*ddzi(Mice)
+
+ do i = 1, Mice+1 
+   dzetai_int(i) = DZETAI(real(i))
+ enddo
+ 
+ do i = 1, Mice
+   dzetai_05int(i) = DZETAI(real(i)+0.5)
+ enddo
+
+ 
+ END SUBROUTINE GRID_CREATE
+
+
+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.
+
+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) ) * &
+        & ( z2(i) - z1(j) )/( z1(j+1) - z1(j) )
+        exit
+      endif
+    enddo
+  endif  
+enddo 
+
+END SUBROUTINE LININTERPOL
+
+
+SUBROUTINE PIECEWISECONSTINT (z1,f1,N1,z2,f2,N2,flag)
+
+ use LAKE_DATATYPES, only : ireals, iintegers
+
+ 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)
+
+! Output variables
+! Reals
+real(kind=ireals), intent(out) :: f2(1:N2)
+
+integer(kind=iintegers), intent(in) :: flag
+
+! Local variables
+integer(kind=iintegers) :: i, j, k ! Loop indices
+
+
+select case (flag)
+case(1_iintegers)
+  k = 1
+  cyc: do i = 2, N1
+    do while (z2(k) <= z1(i))
+      f2(k) = f1(i-1)
+      k = k + 1
+      if (k > N2) exit cyc
+    enddo
+  enddo cyc 
+case(2_iintegers)
+  k = 1
+  cyc1: do i = 2, N1
+    do while (z2(k) <= 0.5*(z1(i) + z1(i-1)) )
+      f2(k) = f1(i-1)
+      k = k + 1
+      if (k > N2) exit cyc1
+    enddo
+  enddo cyc1
+  if (k <= N2) then
+    f2(k:N2) = f1(N1)
+  endif
+end select
+
+END SUBROUTINE PIECEWISECONSTINT
+END MODULE GRID_OPERATIONS_LAKE
+
+ MODULE ALLOCARRAYS_MOD
+ contains
+
+
  !> Subroutine ALLOCARRAYS allocates arrays, sets arrays to reference values, 
  !! and assigns pointers
  SUBROUTINE ALLOCARRAYS(nx,ny)
@@ -18,6 +192,7 @@
  & ngasb
  use RADIATION, only : RadWater, RadIce, RadDeepIce, &
  & fracbands, nbands, RAD_POINTERS
+ use GRID_OPERATIONS_LAKE, only : GRID_CREATE
 
  implicit none
 
@@ -247,8 +422,6 @@
 ! Specification of dzeta-grid (liquid water layers)
 ! and of the dzetai-grid (deep ice and upper ice layers)
  call GRID_CREATE
- !print*, ddz
- !stop
 
  ! Associating pointers to groups of variables
  ! Layers' thicknesses
@@ -311,170 +484,6 @@
  gas%phosph   => phosph
 
  END SUBROUTINE ALLOCARRAYS
+ END MODULE ALLOCARRAYS_MOD
 
 
- SUBROUTINE GRID_CREATE
-
- use LAKE_DATATYPES, only : ireals, iintegers
-
- use ARRAYS_GRID
- use DRIVING_PARAMS
- use DZETA_MOD, only : DZETA, DZETAI
-
- implicit none
-
- integer(kind=iintegers) :: i
- real(kind=ireals) :: dzeta1, dzeta2
-
-!The grid according to Burchard, 2002, pp. 97-98
- do i = 1, M
-   dzeta1 = &
-   & (tanh(d_surf%par)+tanh((d_surf%par+d_bot%par)*real(i-1)/real(M)-d_surf%par))/ &
-   & (tanh(d_surf%par)+tanh(d_bot%par))
-   dzeta2 = &
-   & (tanh(d_surf%par)+tanh((d_surf%par+d_bot%par)*real(i)/real(M)-d_surf%par))/ &
-   & (tanh(d_surf%par)+tanh(d_bot%par))
-   ddz(i) = dzeta2 - dzeta1
- enddo
- 
- ddz05(0) = 0.5d0*ddz(1)
- do i = 1, M-1
-   ddz2(i) = ddz(i) + ddz(i+1)
-   ddz05(i) = 0.5d0*(ddz(i) + ddz(i+1))   
- enddo
- ddz05(M) = 0.5d0*ddz(M)
- 
- do i = 2, M-1
-   ddz052(i) = 0.5d0*ddz(i-1) + ddz(i) + 0.5d0*ddz(i+1)
-   ddz054(i) = 2.d0*ddz052(i)
- enddo
- 
- do i = 1, M+1 
-   dzeta_int(i) = DZETA(real(i))
- enddo
- 
- do i = 1, M
-   dzeta_05int(i) = DZETA(real(i)+0.5)
- enddo
- 
- ddzi(:) = 1.d0/float(Mice)
-
- ddzi05(0) = 0.5d0*ddzi(1)
- do i = 1, Mice-1
-   ddzi05(i) = 0.5d0*(ddzi(i) + ddzi(i+1))
- enddo
- ddzi05(Mice) = 0.5d0*ddzi(Mice)
-
- do i = 1, Mice+1 
-   dzetai_int(i) = DZETAI(real(i))
- enddo
- 
- do i = 1, Mice
-   dzetai_05int(i) = DZETAI(real(i)+0.5)
- enddo
-
- 
- END SUBROUTINE GRID_CREATE
-
-
-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.
-
-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) ) * &
-        & ( z2(i) - z1(j) )/( z1(j+1) - z1(j) )
-        exit
-      endif
-    enddo
-  endif  
-enddo 
-
-END SUBROUTINE LININTERPOL
-
-
-SUBROUTINE PIECEWISECONSTINT (z1,f1,N1,z2,f2,N2,flag)
-
- use LAKE_DATATYPES, only : ireals, iintegers
-
- 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)
-
-! Output variables
-! Reals
-real(kind=ireals), intent(out) :: f2(1:N2)
-
-integer(kind=iintegers), intent(in) :: flag
-
-! Local variables
-integer(kind=iintegers) :: i, j, k ! Loop indices
-
-
-select case (flag)
-case(1_iintegers)
-  k = 1
-  cyc: do i = 2, N1
-    do while (z2(k) <= z1(i))
-      f2(k) = f1(i-1)
-      k = k + 1
-      if (k > N2) exit cyc
-    enddo
-  enddo cyc 
-case(2_iintegers)
-  k = 1
-  cyc1: do i = 2, N1
-    do while (z2(k) <= 0.5*(z1(i) + z1(i-1)) )
-      f2(k) = f1(i-1)
-      k = k + 1
-      if (k > N2) exit cyc1
-    enddo
-  enddo cyc1
-  if (k <= N2) then
-    f2(k:N2) = f1(N1)
-  endif
-end select
-
-END SUBROUTINE PIECEWISECONSTINT
diff --git a/source/model/init_var_mod.f90 b/source/model/init_var_mod.f90
index b29d1b072f3b843041b52400a5afa49e8056a296..82cd5564a4c43a329e3d3d2606b6113be048d378 100644
--- a/source/model/init_var_mod.f90
+++ b/source/model/init_var_mod.f90
@@ -1,6 +1,7 @@
   MODULE INIT_VAR_MOD
 
   use LAKE_DATATYPES, only : ireals, iintegers
+  use GRID_OPERATIONS_LAKE, only : LININTERPOL
 
   contains
   SUBROUTINE INIT_VAR &
diff --git a/source/model/lake.f90 b/source/model/lake.f90
index 168d8de489cc19113968654d009f1279b62ba65f..8769ccfb10088f219f4f3cbc276be043847b1a5d 100644
--- a/source/model/lake.f90
+++ b/source/model/lake.f90
@@ -19,13 +19,14 @@ use INOUT_PARAMETERS, only : &
 & lake_subr_unit_min, &
 & lake_subr_unit_max
 use WATER_DENSITY
-use SOIL_MOD, only : COMSOILFORLAKE
+use SOIL_LAKE_MOD, only : COMSOILFORLAKE
 use SURF_SCHEME1_MOD, only : PBLDAT
 use SURF_SCHEME_INM, only : PBLDAT_INM
 use INOUT, only : fext2, readgrd_lake
 use BL_MOD_LAKE, only : BL_MOD_LAKE_ALLOC
 use SNOWSOIL_MOD, only : SNOWSOIL_MOD_ALLOC
 use METH_OXYG_CONSTANTS, only : SET_METH_OXYG_CONSTANTS
+use ALLOCARRAYS_MOD, only : ALLOCARRAYS
 
 implicit none
 
@@ -175,22 +176,23 @@ use ARRAYS_WATERSTATE
 use ARRAYS_METHANE
 use ARRAYS_OXYGEN
 use MODI_MASSFLUX_CONVECTION
-use CONVECTPAR, only : UnDenGrMix
-use SALINITY, only : S_DIFF, UPDATE_ICE_SALINITY
-use TURB, only : K_EPSILON, ED_TEMP_HOSTETLER, ED_TEMP_HOSTETLER2, ZIRIVMODEL
+use CONVECTPAR_LAKE_MOD, only : UnDenGrMix
+use SALINITY_LAKE_MOD, only : S_DIFF, UPDATE_ICE_SALINITY
+use TURB_LAKE_MOD, only : K_EPSILON, ED_TEMP_HOSTETLER, &
+&                         ED_TEMP_HOSTETLER2, ZIRIVMODEL
 use TURB_CONST, only : min_diff
 use DZETA_MOD, only : VARMEAN
 
-use MOMENTUM, only : &
+use MOMENTUM_LAKE_MOD, only : &
 & MOMENTUM_SOLVER
 
 use VERTADV, only : &
 & VERTADV_UPDATE
 
-use METHANE_MOD, only : &
+use METHANE_LAKE_MOD, only : &
 & METHANE, METHANE_OXIDATION
 
-use SOIL_MOD, only : &
+use SOIL_LAKE_MOD, only : &
 & SOILFORLAKE, &
 & SOIL_COND_HEAT_COEF, &
 & SOILCOLSTEMP, &
@@ -243,7 +245,7 @@ use OXYGEN_MOD, only : &
 & ADDOXPRODCONS, &
 & CHLOROPHYLLA
 
-use CARBON_DIOXIDE_MOD, only : &
+use CARBON_DIOXIDE_LAKE_MOD, only : &
 & CARBON_DIOXIDE
 
 use PHOSPHORUS_MOD, only : &
@@ -257,9 +259,9 @@ use TIMEVAR, only : &
 
 use OUT_MOD
 
-use BUBBLE_MOD, only : BUBBLE, BUBBLEFLUXAVER
+use BUBBLE_LAKE_MOD, only : BUBBLE, BUBBLEFLUXAVER
 
-use CONTROL_POINT, only : CONTROL_POINT_OUT, CONTROL_POINT_IN
+use CONTROL_POINT_LAKE_MOD, only : CONTROL_POINT_OUT, CONTROL_POINT_IN
 
 use TRIBUTARIES, only : TRIBTEMP
 
@@ -270,11 +272,14 @@ use NUMERICS, only : ACCUMM
 use BL_MOD_LAKE
 
 use SNOWSOIL_MOD
+use SNOWTEMP_LAKE, only : SNOWTEMP, SNOW_COND_HEAT_COEF
 
 use RADIATION, only : &
 & RadWater, RadIce, RadDeepIce, &
 & fracbands, nbands, extwatbands
 
+use ZERODIM_MODEL_LAKE_MOD, only : ZERODIM_MODEL
+
 implicit none
 
 !Input variables
diff --git a/source/model/methane_mod.f90 b/source/model/methane_mod.f90
index b5bf9e44a7aa84ab852a808ba256c37042d32f85..a12ef0ea475bb1d798f42be4aa1e45aea37451e5 100644
--- a/source/model/methane_mod.f90
+++ b/source/model/methane_mod.f90
@@ -1,4 +1,4 @@
-MODULE METHANE_MOD
+MODULE METHANE_LAKE_MOD
 
 use NUMERICS, only : PROGONKA
 
@@ -47,7 +47,7 @@ use PHYS_FUNC, only : &
 
 use T_SOLVER_MOD, only : DIFF_COEF   
 
-use BUBBLE_MOD, only : BUBBLEFLUXAVER
+use BUBBLE_LAKE_MOD, only : BUBBLEFLUXAVER
 
 use ARRAYS_SOIL, only : lsh_type
 use ARRAYS_GRID, only : gridsize_type, gridspacing_type
@@ -1464,4 +1464,4 @@ endif ifdeep
  END FUNCTION MEAN_RPROD_NUM_OLDC2    
 
 
- END MODULE METHANE_MOD
+ END MODULE METHANE_LAKE_MOD
diff --git a/source/model/momentum_mod.f90 b/source/model/momentum_mod.f90
index 297946502383fb4240a99941456cbd1a90f33159..62b6d26a73585bdbc93c04e476aa4c1a458ea793 100644
--- a/source/model/momentum_mod.f90
+++ b/source/model/momentum_mod.f90
@@ -1,7 +1,7 @@
-MODULE MOMENTUM
+MODULE MOMENTUM_LAKE_MOD
 
 use NUMERICS, only : PROGONKA, KRON, MATRIXPROGONKA
-use TURB, only : CE_CANUTO, SMOMENT_GALPERIN
+use TURB_LAKE_MOD, only : CE_CANUTO, SMOMENT_GALPERIN
 use NUMERIC_PARAMS, only : vector_length, pi, small_value, minwind
 
 use LAKE_DATATYPES, only : ireals, iintegers
@@ -94,7 +94,7 @@ use WATER_DENSITY, only: &
 & DDENS_DTEMP0, DDENS_DSAL0, &
 & SET_DDENS_DTEMP, SET_DDENS_DSAL
 
-use TURB, only : &
+use TURB_LAKE_MOD, only : &
 & VSMOP_LAKE
 
 implicit none
@@ -2250,8 +2250,7 @@ deallocate(hx1, hx2, hy1, hy2)
 
 END SUBROUTINE RELAYER
 
-
-END MODULE MOMENTUM
+END MODULE MOMENTUM_LAKE_MOD
 
 
 
diff --git a/source/model/out_mod.f90 b/source/model/out_mod.f90
index a8258145690b06c19b9ddf8e6332111e4ccdf909..3a609f02b2eb61a5015b93049f697a178763ac99 100644
--- a/source/model/out_mod.f90
+++ b/source/model/out_mod.f90
@@ -3,6 +3,8 @@ MODULE OUT_MOD
 use LAKE_DATATYPES, only : ireals, iintegers
 use INOUT, only : CHECK_UNIT
 use DRIVING_PARAMS, only : missing_value
+use GRID_OPERATIONS_LAKE, only : LININTERPOL
+  use TIME_LAKE_MOD, only : DATEMINUS, TIMESTR
 
 character, save :: outpath*60 
 
diff --git a/source/model/oxygen_mod.f90 b/source/model/oxygen_mod.f90
index f56db4569d2e1dcc39bd0829d4a3386b398be75c..08031111ac94d5520bb739a6b042eb89bbb75bb0 100644
--- a/source/model/oxygen_mod.f90
+++ b/source/model/oxygen_mod.f90
@@ -10,7 +10,7 @@ use ARRAYS_SOIL, only : lsh_type
 use METH_OXYG_CONSTANTS
 use DRIVING_PARAMS, only : carbon_model
 use UNITS, only : kg_to_g
-use DIFFUSION_MOD, only : DIFFUSION
+use DIFFUSION_LAKE_MOD, only : DIFFUSION
 
 real(4), parameter :: k_PA = 1.7 ! n/d, coefficient in Poole&Atkins formula
 real(4), parameter :: extwat1 = k_PA/2., extwat2 = k_PA/3.5 ! Extinction coefficients
diff --git a/source/model/phosphorus_mod.f90 b/source/model/phosphorus_mod.f90
index e4f1f4257269561b7d2821bbf309685a677680dc..f2552bbc9714c8030bcdce061d9338aa3f87311b 100644
--- a/source/model/phosphorus_mod.f90
+++ b/source/model/phosphorus_mod.f90
@@ -10,7 +10,7 @@ contains
 !> Subroutine performs a time step for phosphorus 
 SUBROUTINE PHOSPHORUS(M,dt,ls,gsp,bath,gas,lsp,diff)
 
-use DIFFUSION_MOD, only : DIFFUSION
+use DIFFUSION_LAKE_MOD, only : DIFFUSION
 use ARRAYS       , only : gas_type
 use ARRAYS_BATHYM, only : bathym, layers_type
 use ARRAYS_GRID  , only : gridspacing_type
diff --git a/source/model/phys_func.f90 b/source/model/phys_func.f90
index 98ba30a46d50f43d6fcc7b62dfee677051035af7..b64b795429f964cc4b0cf3f3e657143ae2fc1cba 100644
--- a/source/model/phys_func.f90
+++ b/source/model/phys_func.f90
@@ -579,7 +579,7 @@ END FUNCTION GSW_RHO
   FUNCTION SINH0(year,month,day,hour,phi)
 
 ! SINH0 is sine of solar angle ( = cosine of zenith angle)   
-  
+  use TIME_LAKE_MOD, only : JULIAN_DAY  
   implicit none
   
   real(kind=ireals) :: sinh0
@@ -597,8 +597,6 @@ END FUNCTION GSW_RHO
    real(kind=ireals) phi_rad
 
   integer(kind=iintegers) nday
-
-  integer(kind=iintegers), external:: JULIAN_DAY
  
   pi=4.*datan(1.d0)
 
diff --git a/source/model/salinity_mod.f90 b/source/model/salinity_mod.f90
index c108749c3387be3982e47b6d6d5fafc097efe3fd..20d3510220781196bcac4157295094067bb27ecc 100644
--- a/source/model/salinity_mod.f90
+++ b/source/model/salinity_mod.f90
@@ -1,4 +1,4 @@
-MODULE SALINITY
+MODULE SALINITY_LAKE_MOD
 
 use NUMERICS, only : PROGONKA, CHECK_PROGONKA, IND_STAB_FACT_DB
 use NUMERIC_PARAMS, only : vector_length, small_value
@@ -503,4 +503,4 @@ END FUNCTION RESID
 END SUBROUTINE UPDATE_ICE_SALINITY
 
 
-END MODULE SALINITY
+END MODULE SALINITY_LAKE_MOD
diff --git a/source/model/snowcalc.f90 b/source/model/snowcalc.f90
index 445735551c44be62a83f20b20f4e972a7ba31e0e..fa1298caa96201e8d170b2952bba0f29ffa43485 100644
--- a/source/model/snowcalc.f90
+++ b/source/model/snowcalc.f90
@@ -1,3 +1,6 @@
+MODULE SNOW_LAKE
+
+contains
       SUBROUTINE SNOW_CALC(t2,Erad_sc,ts,snow_sc,iyear,imonth, &
       & iday,CCT,pmelt,pheat,dt)
       
@@ -703,3 +706,5 @@
       
       return
       END SUBROUTINE addPrecipToSnow
+
+END MODULE SNOW_LAKE
diff --git a/source/model/snowtemp.f90 b/source/model/snowtemp.f90
index 45a2961c1831c5ff015442778a427dc11921da10..936e8f9a32c7ed1b17dd60c31b679efd8040cd61 100644
--- a/source/model/snowtemp.f90
+++ b/source/model/snowtemp.f90
@@ -1,3 +1,8 @@
+MODULE SNOWTEMP_LAKE
+
+use SNOW_LAKE, only : SNOW_CALC, ADDPRECIPTOSNOW
+
+contains
 SUBROUTINE SNOWTEMP(ix,iy,nx,ny,year,month,day,hour,snowmass, &
 & snowmass_init,a,b,c,d,Temp,phi,fetch,dt)
 
@@ -198,3 +203,5 @@ end do
 !print*, minval(cs), maxval(cs)
 
 END SUBROUTINE SNOW_COND_HEAT_COEF
+
+END MODULE SNOWTEMP_LAKE
diff --git a/source/model/soil_mod.f90 b/source/model/soil_mod.f90
index 275f0831b8e08abac7966c81c013aa718351a703..e16063cb5f64621c1a92371d24a3430d0afae364 100644
--- a/source/model/soil_mod.f90
+++ b/source/model/soil_mod.f90
@@ -1,4 +1,4 @@
-MODULE SOIL_MOD
+MODULE SOIL_LAKE_MOD
 
 use LAKE_DATATYPES, only : ireals, iintegers
 use NUMERICS, only : PROGONKA, STEP
@@ -820,7 +820,7 @@ use ARRAYS_TURB, only: eps1
 use ARRAYS_OXYGEN, only : sodbot
 use ARRAYS, only : gas, dt_inv
 
-use METHANE_MOD, only : METHANE
+use METHANE_LAKE_MOD, only : METHANE
 
 use METH_OXYG_CONSTANTS, only : CH4_exp_growth
 
@@ -1598,5 +1598,5 @@ END FUNCTION OVERLAY
 END SUBROUTINE LATERHEAT
 
 
-END MODULE SOIL_MOD
+END MODULE SOIL_LAKE_MOD
 
diff --git a/source/model/surf_scheme2_mod.f90 b/source/model/surf_scheme2_mod.f90
index ef65ff41d0ae2be48c80e9ec68ea809a7806b527..44c379ef2ca4710d4b133b298e987ffe7c414618 100644
--- a/source/model/surf_scheme2_mod.f90
+++ b/source/model/surf_scheme2_mod.f90
@@ -1,3 +1,6 @@
+MODULE RICHPBL_LAKE_MOD
+
+contains
       SUBROUTINE RichPBL(h,le,tau,Tsurf1,humsurf,z0)
       
 !     RichPBL calculates heat fluxes according to Louis parameterization (Louis, 1079)    
@@ -48,8 +51,8 @@
 
       real(kind=ireals) z0hz0,Tsurf1,humsurf,z,xmu,xfh,ts,h,le,leg, &
      & ta,qa,ps,ua,z0h,z0,hu,veg,ztvi,rhoa,zdepl,zua,ztvis,Ri,&
-     & valnon,zeps,zch,zsta,iyra,zdi,rra,zdim,cdm,zds,chstar2,qs,ph2, &
-     & cdh,cmstar2,pm2,zdsm,xhu,tau
+     & valnon,zeps,zch,zsta,iyra,zdi,rra,zdim,cdm,zds,qs, &
+     & cdh,zdsm,xhu,tau
 
       SAVE
       
@@ -252,3 +255,6 @@
       pm2    =0.5233-0.0815*xmu+.0135*xmu*xmu-.001*xmu*xmu*xmu
       return
       end
+
+
+END MODULE RICHPBL_LAKE_MOD
diff --git a/source/model/surf_scheme3.f90 b/source/model/surf_scheme3.f90
index 3e3bd0e24ca020d254f5cbb3350af66bc4e0c7a7..08c3831595412dd67edcbdd6a3472366309b6233 100644
--- a/source/model/surf_scheme3.f90
+++ b/source/model/surf_scheme3.f90
@@ -1,3 +1,6 @@
+MODULE SURF_SCHEME3_LAKE_MOD
+
+contains
       SUBROUTINE SURF_SCHEME3(u2,t2,t1,q2,q1,z1,z2,ro,l1,   &
       & hflux, Elatent,tau)
       
@@ -486,3 +489,5 @@
           return
           
           END SUBROUTINE SURF_SCHEME3
+
+END MODULE SURF_SCHEME3_LAKE_MOD
diff --git a/source/model/t_0dim.f90 b/source/model/t_0dim.f90
index 69367198466457cb3084f544ce8db424d48af9b7..bca68a13c560d22e6744bcc4d2578c6cded962b0 100644
--- a/source/model/t_0dim.f90
+++ b/source/model/t_0dim.f90
@@ -1,3 +1,6 @@
+MODULE ZERODIM_MODEL_LAKE_MOD
+
+contains
 SUBROUTINE ZERODIM_MODEL &
 & (h, dt, &
 & shortwave, longwave, tempair, humair, pressure, &
@@ -108,3 +111,5 @@ Temp0 = Temp0 + dt*(eflux - botflux)/(cw_m_row0*h)
 !Temp0 = Temp0 + dt*eflux/(cw_m_row0*h)
 
 END SUBROUTINE ZERODIM_MODEL
+
+END MODULE ZERODIM_MODEL_LAKE_MOD
diff --git a/source/model/trib.f90 b/source/model/trib.f90
index 0e76961acaf99e9768192594cc995cabc813c05f..f5c2734b411f0e5e7d8e71be3b93c9063fcafce4 100644
--- a/source/model/trib.f90
+++ b/source/model/trib.f90
@@ -8,6 +8,8 @@ use ARRAYS_WATERSTATE, only : waterstate_type
 use TIMEVAR, only : day_sec
 use NUMERIC_PARAMS, only : small_value
 
+use GRID_OPERATIONS_LAKE, only : PIECEWISECONSTINT
+
 integer(kind=iintegers), save :: itriblev
 
 contains
diff --git a/source/model/turb_mod.f90 b/source/model/turb_mod.f90
index 6ad5857b63c3f9a86ff7dba997e8287d017f2d12..500c4d6d63d17dc0302fde58e16f6aa7c851d1a7 100644
--- a/source/model/turb_mod.f90
+++ b/source/model/turb_mod.f90
@@ -1,4 +1,4 @@
-MODULE TURB
+MODULE TURB_LAKE_MOD
 
 use NUMERICS, only : PROGONKA, IND_STAB_FACT_DB
 use NUMERIC_PARAMS, only : vector_length, small_value
@@ -1887,4 +1887,4 @@ endif
 
 END SUBROUTINE ZIRIVMODEL
 
-END MODULE TURB
+END MODULE TURB_LAKE_MOD
diff --git a/source/modi_massflux_convection.mod b/source/modi_massflux_convection.mod
new file mode 100644
index 0000000000000000000000000000000000000000..2493a0e573205876c928bd9c710ec8e99363dd77
Binary files /dev/null and b/source/modi_massflux_convection.mod differ
diff --git a/source/momentum_lake_mod.mod b/source/momentum_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..dc7e085adffd60faa6a90abb5b2e7db9192f3620
Binary files /dev/null and b/source/momentum_lake_mod.mod differ
diff --git a/source/numeric_params.mod b/source/numeric_params.mod
new file mode 100644
index 0000000000000000000000000000000000000000..a3d9888ddada4783713f93c08fbb15419fae97d8
Binary files /dev/null and b/source/numeric_params.mod differ
diff --git a/source/numerics.mod b/source/numerics.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f42436e6417779c7ca9031faabb028d02fc490e6
Binary files /dev/null and b/source/numerics.mod differ
diff --git a/source/out_mod.mod b/source/out_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1fd4ea97a61ad04f0b8f6b831db08f965f09c817
Binary files /dev/null and b/source/out_mod.mod differ
diff --git a/source/oxygen_mod.mod b/source/oxygen_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3d8ac8eac5aa3f46ee186f7fabf2d54a3b3b0787
Binary files /dev/null and b/source/oxygen_mod.mod differ
diff --git a/source/phosphorus_mod.mod b/source/phosphorus_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1ee91941764aa71959d45dbd56120d6db1879875
Binary files /dev/null and b/source/phosphorus_mod.mod differ
diff --git a/source/phys_constants.mod b/source/phys_constants.mod
new file mode 100644
index 0000000000000000000000000000000000000000..7da5bbd928482a874de11c32b9a66c2d4f3491a7
Binary files /dev/null and b/source/phys_constants.mod differ
diff --git a/source/phys_func.mod b/source/phys_func.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d512d2902bc0101c108bd7fd6f37db96914c9993
Binary files /dev/null and b/source/phys_func.mod differ
diff --git a/source/phys_parameters.mod b/source/phys_parameters.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3902399839e946a998b275da1f6843f5de650489
Binary files /dev/null and b/source/phys_parameters.mod differ
diff --git a/source/postprocessing_lake.mod b/source/postprocessing_lake.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1d1faf3e0d40932337be1f77481dd5ec3454052c
Binary files /dev/null and b/source/postprocessing_lake.mod differ
diff --git a/source/radiation.mod b/source/radiation.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1980310ac0980d03f01d5936ea76c2963431d4be
Binary files /dev/null and b/source/radiation.mod differ
diff --git a/source/richpbl_lake_mod.mod b/source/richpbl_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..4621bb8b6b5da4affb38b691417a3177c91897ba
Binary files /dev/null and b/source/richpbl_lake_mod.mod differ
diff --git a/source/salinity_lake_mod.mod b/source/salinity_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..1534b716403996f0e49b766f1c0d2017de3f5b09
Binary files /dev/null and b/source/salinity_lake_mod.mod differ
diff --git a/source/seiches_param.mod b/source/seiches_param.mod
new file mode 100644
index 0000000000000000000000000000000000000000..3ea27e800e73a4175ff0c010ac13619303df0f3f
Binary files /dev/null and b/source/seiches_param.mod differ
diff --git a/source/sfcflx.mod b/source/sfcflx.mod
new file mode 100644
index 0000000000000000000000000000000000000000..e26a9eac705eafeb9398cebf687e28d8aedb7a9c
Binary files /dev/null and b/source/sfcflx.mod differ
diff --git a/source/model/time.f90 b/source/shared/time.f90
similarity index 94%
rename from source/model/time.f90
rename to source/shared/time.f90
index 1a83c425cb35e5716d655cbe7f779891ebc203d1..decac8274dfd38e9c656f91f4c0a39efce13ecea 100644
--- a/source/model/time.f90
+++ b/source/shared/time.f90
@@ -1,4 +1,6 @@
+MODULE TIME_LAKE_MOD
 
+contains
 SUBROUTINE JULIAN_DATE(year,month,day,hour,dt)
 
 ! The subroutine JULIAN_DATE updates the julian date
@@ -219,3 +221,5 @@ else
 endif
 
 END SUBROUTINE TIMEC
+
+END MODULE TIME_LAKE_MOD
diff --git a/source/snow_lake.mod b/source/snow_lake.mod
new file mode 100644
index 0000000000000000000000000000000000000000..e0b3a29bed44770d169b345c51fd9f45143d4a82
Binary files /dev/null and b/source/snow_lake.mod differ
diff --git a/source/snowsoil_mod.mod b/source/snowsoil_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..56fa9d460f84d9f560efb4de0f0a34fc7c97ec2a
Binary files /dev/null and b/source/snowsoil_mod.mod differ
diff --git a/source/snowtemp_lake.mod b/source/snowtemp_lake.mod
new file mode 100644
index 0000000000000000000000000000000000000000..6c8a3eea559e44c4f9059490284d225a3287ef37
Binary files /dev/null and b/source/snowtemp_lake.mod differ
diff --git a/source/soil_lake_mod.mod b/source/soil_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0742636608ae6c420a077bc653687cb9092ba557
Binary files /dev/null and b/source/soil_lake_mod.mod differ
diff --git a/source/surf_scheme1_mod.mod b/source/surf_scheme1_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..15f64ff4370f9974edbe84c9eccac708729be023
Binary files /dev/null and b/source/surf_scheme1_mod.mod differ
diff --git a/source/surf_scheme3_lake_mod.mod b/source/surf_scheme3_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..75e59404ed012d0c76d440e5cd4892f40a990479
Binary files /dev/null and b/source/surf_scheme3_lake_mod.mod differ
diff --git a/source/surf_scheme_inm.mod b/source/surf_scheme_inm.mod
new file mode 100644
index 0000000000000000000000000000000000000000..eb47a97b40961ec3e40b4e05ce14f3d7adbb859d
Binary files /dev/null and b/source/surf_scheme_inm.mod differ
diff --git a/source/t_solver_mod.mod b/source/t_solver_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..86857f48fe28ece2b2a5282803960e553b6d8580
Binary files /dev/null and b/source/t_solver_mod.mod differ
diff --git a/source/time_lake_mod.mod b/source/time_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..06d8ae59ad2e806d7e31733011a9018c66814bfa
Binary files /dev/null and b/source/time_lake_mod.mod differ
diff --git a/source/timevar.mod b/source/timevar.mod
new file mode 100644
index 0000000000000000000000000000000000000000..d30b1c9c0a16c2aa8fcc6a75174858455700b144
Binary files /dev/null and b/source/timevar.mod differ
diff --git a/source/tributaries.mod b/source/tributaries.mod
new file mode 100644
index 0000000000000000000000000000000000000000..0fbb7d6fa5d5279690981fec81ceedac4e4a81cf
Binary files /dev/null and b/source/tributaries.mod differ
diff --git a/source/turb_const.mod b/source/turb_const.mod
new file mode 100644
index 0000000000000000000000000000000000000000..f01fa35952cefb611a10c680249f3f81dc49c6cb
Binary files /dev/null and b/source/turb_const.mod differ
diff --git a/source/turb_lake_mod.mod b/source/turb_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..741bc35238a321c710669f727f20c26e1f1622b5
Binary files /dev/null and b/source/turb_lake_mod.mod differ
diff --git a/source/units.mod b/source/units.mod
new file mode 100644
index 0000000000000000000000000000000000000000..542f91a0df1ab9594d247e2a72e3274bec816230
Binary files /dev/null and b/source/units.mod differ
diff --git a/source/vertadv.mod b/source/vertadv.mod
new file mode 100644
index 0000000000000000000000000000000000000000..5f1f3f016dcba8fa875307aa9580088085d598ab
Binary files /dev/null and b/source/vertadv.mod differ
diff --git a/source/water_density.mod b/source/water_density.mod
new file mode 100644
index 0000000000000000000000000000000000000000..2146a243b2597558e4500fc4c31719a7b467d93c
Binary files /dev/null and b/source/water_density.mod differ
diff --git a/source/zerodim_model_lake_mod.mod b/source/zerodim_model_lake_mod.mod
new file mode 100644
index 0000000000000000000000000000000000000000..83864f0a14164b60e716dc01b244b919b24f527a
Binary files /dev/null and b/source/zerodim_model_lake_mod.mod differ