Skip to content
Snippets Groups Projects
lake_modules.f90 103 KiB
Newer Older
  • Learn to ignore specific revisions
  • Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    real(kind=ireals), parameter :: hour_sec = 60.*60.
    real(kind=ireals), parameter :: day_sec = 24.*60.*60.
    real(kind=ireals), parameter :: year_sec  = 60.*60.*24.*365.
    
    real(kind=ireals), parameter :: month_sec  = 60.*60.*24.*365.25/12.
    
    real(kind=ireals), parameter :: hour_min = 60.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    END MODULE TIMEVAR
    
    MODULE UNITS
    
    use LAKE_DATATYPES, only : ireals
    
    real(kind=ireals), parameter :: kg_to_g = 1.E+3
    
    real(kind=ireals), parameter :: m3_to_l = 1.E+3
    
    real(kind=ireals), parameter :: m_to_nm = 1.E+9
    
    MODULE NUMERIC_PARAMS
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    integer(kind=iintegers), parameter :: KL = 1, NT = 52592, &
    & Num_Soil = 11, Num_Veget = 13, ML = 131, MS = 100 
    integer(kind=iintegers), target, save :: ML_ = ML, MS_ = MS
    real(kind=ireals), parameter :: dznorm = 0.04, dzmin = 0.015, &
    & UpperLayer = 0.02 ! UpperLayer = 0.1, UpperLayer = 0.008 
    
    integer(kind=iintegers), parameter :: vector_length = 250
    real(kind=ireals),    parameter :: pi = 4.*datan(1.d0)
    
    real(kind=ireals), parameter :: small_value = 1.d-20
    
    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-2
    
    
     !   ML --- total number of layers
     !   MS --- maximal number of layers in snow
     !   dznorm --- standard layer depth in snow (m)
     !   dzmin --- min layer depth in snow (m)
     !   depth --- depth of the lowest layer in soil (m)
     !   UpperLayer --- upper soil layer depth (m)
    
    END MODULE NUMERIC_PARAMS
    
    
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    real(kind=ireals) :: whc, hcond, sncr
    
    
    SAVE
    
    contains
    SUBROUTINE DEFINE_PHYS_PARAMETERS()
    implicit none
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    !cccccccccccccccccccccc SNOW cccccccccccccccccccccccc
    
    whc = 0.04   ! water holding capacity of snow
    hcond = 0.01 ! hydraulic conductivity in snow (m/s)
    SNCR = 0.01  ! M OF WATER, FROM WHICH SNOW COVERS ALL THE CELL
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
    
    END SUBROUTINE DEFINE_PHYS_PARAMETERS
    END MODULE PHYS_PARAMETERS
    
    MODULE PHYS_CONSTANTS
    
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    use NUMERIC_PARAMS, only : pi
    
    
    real(kind=ireals) :: g
    real(kind=ireals) :: omega
    
    real(kind=ireals) :: roa0, row0, rofrsn
    
    real(kind=ireals), target :: roi
    real(kind=ireals) :: cw, csnow
    real(kind=ireals), target :: ci
    real(kind=ireals), target :: Lwi
    real(kind=ireals) :: Lwv, Liv
    real(kind=ireals) :: kappa
    real(kind=ireals) :: sigma
    real(kind=ireals) :: cp
    real(kind=ireals) :: aMag, bMag
    real(kind=ireals) :: Rd, Rwv, R_univ
    real(kind=ireals) :: Avogadro_number
    
    real(kind=ireals) :: Planck_const, Planck_const_
    
    real(kind=ireals) :: clight
    real(kind=ireals) :: lambda_PAR0
    real(kind=ireals) :: short2PAR
    
    real(kind=ireals) :: niu_atm, niu_wat
    real(kind=ireals) :: surf_tension_wat
    real(kind=ireals) :: roughness0
    
    real(kind=ireals) :: charnock_const
    
    real(kind=ireals) :: albedoofice, albedoofsnow, albedoofwater, albedoofsoil
    real(kind=ireals) :: albedoofice_lw, albedoofsnow_lw, albedoofwater_lw, albedoofsoil_lw
    
    real(kind=ireals), target :: emissivityofice, emissivityofwater 
    real(kind=ireals), target :: emissivityofsnow, emissivityofsoil
    
    real(kind=ireals) :: aMagw, bMagw
    real(kind=ireals) :: aMagi, bMagi
    real(kind=ireals) :: alsal, almeth, aloxyg, alcarbdi
    real(kind=ireals) :: sabs
    
    real(kind=ireals) :: sabs0
    
    real(kind=ireals) :: Kelvin0
    real(kind=ireals) :: pref
    real(kind=ireals) :: esatsurf0
    
    real(kind=ireals) :: salice0, poricebot, asalice
    
    real(kind=ireals) :: snmeltwat, icetopfr
    
    real(kind=ireals) :: z0_bot, clin_bot, cquad_bot
    
    !Derived constants
    
    real(kind=ireals) :: ci_m_roi, cw_m_row0
    real(kind=ireals) :: row0_m_Lwi, roi_m_Lwi
    real(kind=ireals) :: row0_d_roi, roi_d_row0
    real(kind=ireals) :: roa0_d_row0, roa0_d_roi
    real(kind=ireals) :: row0_d_roa0
    real(kind=ireals) :: Rd_d_cp, Rd_d_Rwv
    
    real(kind=ireals) :: Wm22Em2, Planck_const_m_clight_m_Avogadro_number
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    real(kind=ireals) :: g_d_Kelvin0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    
    SAVE
    
    contains
    SUBROUTINE DEFINE_PHYS_CONSTANTS 
    implicit none
    
    g     = 9.814    ! acceleration due to gravity,   m/s**2
    omega = 7.29d-5  ! angular velocity of the Earth's rotation,  1/s
    sigma = 5.67d-8  ! Stefan-Boltzman constant,W/(m**2*K**4)
    
    roa0  = 1.273    ! reference density of air,kg/m**3
    row0  = 1000.    ! reference density of water,    kg/m**3
    
    rofrsn  = 100.   ! reference density of fresh snow,    kg/m**3
    
    roi   = 917.     ! density of ice,    kg/m**3
    
    cw    = 3990.    ! heat capacity of water,  J/(kg*K)
    ci    = 2150.    ! heat capacity of ice,    J/(kg*K)
    cp    = 1005.    ! heat capacity of air at constant pressure, J/(kg*K)
    csnow = 2100.    ! the snow specific heat content,J/(kg K)
    
    
    lamair = 0.024   ! molecular conductivity of air, J/(m*s*K)
    lami  = 2.2      ! molecular conductivity of ice, J/(m*s*K)
    
    lamw0 = 0.561    ! molecular conductivity of water,     J/(m*s*K)
    
    lamsal0 = lamw0/(cw*row0)*1.E-2 !molecular diffusivity coefficient for salinity, m**2/s
    
    niu_atm = 1.8d-5 ! molecular viscosity of air,    m**2/s
    niu_wat = 1.307d-6 ! molecular viscosity of water at T = 10 C,m**2/s
    surf_tension_wat = 72.d-3 ! surface tension of water at 25 C, N/m
    
    
    Lwi   = 333500. !267900.  ! latent heat of freezing, J/kg
    
    Lwv   = 2501000. ! latent heat of evaporation,    J/kg
    Liv   = 2834600. ! latent heat of sublimation,    J/kg
    
    
    Rd    = 287.05   ! gas constant for dry air,J/(kg*K)
    Rwv   = 461.91   ! gas constant for water vapor,  J/(kg*K)
    R_univ = 8.31441 ! universal gas constant,  J/(mol*K)
    Avogadro_number = 6.022141d+23 !mol**(-1)
    Planck_const = 6.62607d-34     !J*s
    
    Planck_const_ = Planck_const/(2.*pi) !J*s
    
    clight = 3.d+8   ! light speed in vacuum,   m / s
    lambda_PAR0 = 550.d-9 ! ~the wavelength of the PAR spectrum center, m
    short2PAR = 0.48 ! the global to PAR radiation ratio,   n/d
    
    aMagw = 7.6326   ! coefficient for Magnus formula for water surface
    bMagw = 241.9    ! coefficient for Magnus formula for water surface
    aMagi = 9.5! coefficient for Magnus formula for ice   surface
    bMagi = 265.5    ! coefficient for Magnus formula for ice   surface
    
    albedoofwater = 0.06 ! 0.07 ! albedo of water surface for visible radiation, n/d ! 0.07
    
    albedoofice   = 0.40 ! 0.5! albedo of ice   surface,    n/d ! 0.35
    
    albedoofsnow  = 0.85 ! 0.8! albedo of snow  surface,    n/d ! 0.62
    
    albedoofsoil  = 0.3  ! albedo of soil  surface,   n/d
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    !     albedo_lw ~ 1 - emissivity
    
    albedoofwater_lw = 0.02 ! albedo of water surface for atmospheric radiation, n/d
    
    albedoofice_lw   = 0.03  ! albedo of ice surface for atmospheric radiation,   n/d
    
    albedoofsnow_lw  = 0.01  ! albedo of snow surface for atmospheric radiation,  n/d
    
    albedoofsoil_lw  = 0.0  ! albedo of soil surface for atmospheric radiation,  n/d
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    emissivityofwater = 0.98 !emissivity of water surface, n/d ! 0.95
    
    emissivityofice   = 0.97 ! emissivity of ice   surface, n/d ! 0.95
    
    emissivityofsnow  = 0.99 ! emissivity of snow  surface, n/d
    emissivityofsoil  = 0.9  ! emissivity of soil  surface, n/d
    
    roughness0 = 1.d-4 ! the reference roughness of water, m
    
    charnock_const = 0.0144 ! Charnock constant, n/d, 0.0144 - classical value
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
          
    
    Kelvin0 = 273.15
    pref    = 1.d+5 ! Reference atmospheric pressure, Pa
    esatsurf0 = 610.7 ! Saturation partial water vapor pressure at 0 Celsius, Pa
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    sabs   = 0. !0.35 !0.4 ! the part of solar radiation, absorbed at the water surface = fraction of NIR 
    
    alsal  = 1.  ! the ratio of eddy diffusivity for salinity to that of heat
    almeth = 1.  ! the ratio of eddy diffusivity for dissolved methane to that of heat
    aloxyg = 1.  ! the ratio of eddy diffusivity for dissolved oxygen to that of heat
    alcarbdi = 1.  ! the ratio of eddy diffusivity for dissolved carbon dioxode to that of heat
    
    salice0 = 0.! 4.d-3 ! reference salinity of ice
    
    poricebot = 0.05 ! porosity of ice at the bottom of ice layer
    
    asalice = 1./(86400.*365.) !Characteristic time of salt removal from ice
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    snmeltwat = 0. ! The fraction of snowmelt transferred to water through cracks in ice, 
    
                   ! the rest is partly frozen on a top of the ice layer
    
    
    icetopfr = 1. ! The switch for freezing of water at the ice top
    
    z0_bot = 1.3d-4 ! Roughness parameter of bottom surface, m
    
    clin_bot  = 4.E-3 !7.d-5 !7.d-2 !5.d-3 !  Linear friction coefficient at the sloping bottom
    
    cquad_bot = 8.E-2 !7.E-2 !7.d-5 !7.d-2 !5.d-3 ! Quadratic friction coefficient at the sloping bottom
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    ci_m_roi  = ci * roi
    cw_m_row0 = cw * row0
    row0_m_Lwi = row0 * Lwi
    roi_m_Lwi = roi * Lwi
    
    row0_d_roi = row0/roi
    roi_d_row0 = roi/row0
    roa0_d_row0 = roa0/row0
    roa0_d_roi = roa0/roi
    row0_d_roa0 = row0/roa0
    
    Rd_d_cp = Rd/cp
    Rd_d_Rwv = Rd/Rwv
    
    Wm22Em2 = lambda_PAR0/(Avogadro_number*Planck_const*clight) ! converting multiplyer from W/m2 to Einstein/(m*m*s)
                                                                ! for PAR region
    
    Planck_const_m_clight_m_Avogadro_number = Planck_const_ * clight * Avogadro_number
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    g_d_Kelvin0 = g/Kelvin0
    
    
    END SUBROUTINE DEFINE_PHYS_CONSTANTS
    END MODULE PHYS_CONSTANTS
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    real(kind=ireals) :: tempair
    real(kind=ireals) :: relhum
    real(kind=ireals) :: humair
    real(kind=ireals) :: pressure
    real(kind=ireals) :: hw
    real(kind=ireals) :: xlew
    real(kind=ireals) :: cdmw
    real(kind=ireals) :: hw_input, xlew_input, cdmw_input, surfrad_input
    real(kind=ireals) :: uwind, vwind
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    real(kind=ireals) :: wind, wind10, windwork
    
    real(kind=ireals) :: netrad
    
    real(kind=ireals) :: Radbal
    real(kind=ireals) :: hflux
    real(kind=ireals) :: surfrad
    real(kind=ireals) :: hbal
    real(kind=ireals) :: zref
    real(kind=ireals) :: botflux
    
    real(kind=ireals) :: shortwave, sabspen
    
    real(kind=ireals) :: precip
    real(kind=ireals) :: Sflux0
    real(kind=ireals) :: Elatent
    
    real(kind=ireals) :: Radbal_surf, ShortBal_toplayer
    
    real(kind=ireals) :: eflux, eflux0, eflux0_kinem
    real(kind=ireals) :: hskin
    real(kind=ireals) :: tau, tau_wav
    
    real(kind=ireals) :: velfrict, velfrict_prev, velfrict_bot
    
    real(kind=ireals) :: cdm1
    real(kind=ireals) :: turb_density_flux
    real(kind=ireals) :: cloud
    
    real(kind=ireals) :: discharge(1:2)
    
    
    ! Former common-block wind
    real(kind=ireals) :: urel, vrel, u, v
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
         
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    
    !> Radiation data structures and operations on them
    MODULE RADIATION
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    implicit none
    
    type, public :: rad_type
      sequence
      integer(kind=iintegers) :: nbands, nlevels
    
      real(kind=ireals), allocatable       :: extinct(:,:), extinct_CDOC(:), flux(:,:), integr(:)
    
      procedure(RAD_UPDATE), pointer, pass :: RAD_UPDATE
      procedure(RAD_INIT)  , pointer, pass :: RAD_INIT
    
    endtype rad_type
    type(rad_type) :: RadWater, RadIce, RadDeepIce
    
    
    !>@todo: check wavelength bound between "IR-A" and "IR-B", ext coefs and energy fractions
    integer(kind=iintegers), parameter :: nbands = 4 !wavebands: UV, PAR, NIR(IR-A), IR-B
    real(kind=ireals), parameter :: bandbounds(1:nbands+1) = (/280.,400.,700.,1000.,3000./) !Bounds between bands, nm
    
    real(kind=ireals), parameter :: extwatbands(1:nbands) = (/3.13,0.81,1.73,1.087E+3/) !Extinction coefficients for bands
    real(kind=ireals), parameter :: fracbands(1:nbands) = (/0.046,0.429,0.208,0.317/) !The energy fraction of bands in incident shortwave radiation
    
    contains
    
    !>Subroutine integrates Beer-Lambert equation in all spectral bands
    SUBROUTINE RAD_UPDATE(this,thickness,flux0)
    implicit none
    type(rad_type), intent(inout) :: this
    !Input variables
    real(kind=ireals), intent(in) :: thickness(0:this%nlevels-1)
    real(kind=ireals), intent(in) :: flux0(this%nbands)
    !Local variables
    integer(kind=iintegers) :: i, j
    
    this%flux(0,1:this%nbands) = flux0(1:this%nbands)
    this%integr(0) = sum(this%flux(0,1:this%nbands))
    do j = 1, this%nbands
      do i = 1, this%nlevels
        this%flux(i,j) = this%flux(i-1,j) * exp(-this%extinct(i,j)*thickness(i-1))
      enddo
    enddo
    forall(i=1:this%nlevels) this%integr(i) = sum(this%flux(i,1:this%nbands))
    END SUBROUTINE RAD_UPDATE
    
    !> Initializes the radiation variable structure
    
    SUBROUTINE RAD_INIT(this,nbands,nlevels,extCDOC)
    
    type(rad_type), intent(inout) :: this
    integer(kind=iintegers), intent(in) :: nbands, nlevels
    
    logical, optional, intent(in) :: extCDOC 
    
    this%nbands = nbands
    this%nlevels = nlevels
    allocate(this%flux   (0:nlevels,this%nbands))
    allocate(this%extinct(1:nlevels,this%nbands))
    
    if (present(extCDOC)) then
      if (extCDOC) allocate(this%extinct_CDOC(1:nlevels))
    endif
    
    allocate(this%integr (0:nlevels            ))
    END SUBROUTINE RAD_INIT
    
    
    SUBROUTINE RAD_POINTERS
    implicit none
    RadWater   % RAD_UPDATE => RAD_UPDATE
    RadWater   % RAD_INIT   => RAD_INIT
    RadIce     % RAD_UPDATE => RAD_UPDATE
    RadIce     % RAD_INIT   => RAD_INIT
    RadDeepIce % RAD_UPDATE => RAD_UPDATE
    RadDeepIce % RAD_INIT   => RAD_INIT
    END SUBROUTINE RAD_POINTERS
    
    use LAKE_DATATYPES, only : ireals, iintegers
    
    !use DRIVING_PARAMS
    
    implicit none
    
     ! Water state variables
    
     real(kind=ireals), allocatable, target :: Tw1(:),Tw2(:),Ti1(:),Ti2(:),Tis1(:), &
    
     & Tis2(:),RS(:),lamw(:),lamw_back(:),lamsal(:),Sal1(:),Sal2(:), &
    
     & preswat(:), salice(:), porice(:), ci_m_roi_v(:), lami_v(:)
    
     !real(kind=ireals), pointer :: SR(:),SRi(:),SRdi(:)
    
     type, public :: waterstate_type
       sequence
       real(kind=ireals), pointer :: Tw1(:),Tw2(:),Ti1(:),Ti2(:),Tis1(:), &
    
       & Tis2(:),RS(:),SR(:),lamw(:),lamsal(:),Sal1(:),Sal2(:),row(:), &
    
       & SRi(:),SRdi(:),preswat(:),u1(:),v1(:),u2(:),v2(:),wArea(:),w(:),extwatarr(:), &
    
       & salice(:), porice(:), ci_m_roi_v(:), lami_v(:), lamw_back(:)
    
     endtype waterstate_type
     type(waterstate_type) :: wst
    
    
     !type, public :: intbal_type
     !  sequence
     !  real(kind=ireals), allocatable :: terms(:), termsdt(:)
     !  real(kind=ireals) :: storage, storage0, dstorage, resid
     !  integer(kind=iintegers) :: nterms
     !  procedure(intbal_update), pointer, nopass :: intbal_update => INTBAL_UPDATE
     !endtype intbal_type
     !
     !contains
     !
     !SUBROUTINE INTBAL_UPDATE(this,terms,storage1,storage2,dt)
     !implicit none
     !!Input/output variables
     !type(intbal_type) :: this
     !real(kind=ireals), intent(in), dimension(:) :: terms
     !real(kind=ireals), intent(in) :: storage1, storage2, dt
     !!Local variables
     !logical, save :: firstcall = .true.
     !integer(kind=iintegers) :: i
    
     !if (firstcall) then
     !  this%storage0 = storage1
     !endif
    
     !!do i = 1, this%nterms
     !this%terms(:) = terms(:)
     !this%termsdt(:) = this%termsdt(:) + terms(:)*dt
     !this%storage = storage2
     !this%dstorage = storage2 - this%storage0
     !this%resid = this%dstorage - sum(this%termsdt(:))
     !!enddo
    
     !if (firstcall) firstcall = .false.
     !END SUBROUTINE INTBAL_UPDATE
    
    MODULE ARRAYS_GRID
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    implicit none
    
    ! Grid size group 
    integer(kind=iintegers), save, target :: nsoilcols ! Number of soil columns per one lake
    type, public :: gridsize_type
    
      integer(kind=iintegers), pointer :: M, Mice, ns, ms, ml, nsoilcols, nx, ny
      integer(kind=iintegers) :: ix, iy, isoilcol
    
    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(:)
    real(kind=ireals), allocatable, target :: dzetai_int(:), dzetai_05int(:)      
    
    integer(kind=iintegers), allocatable, target :: isoilcolc(:), ksoilcol(:)
    
    type, public :: gridspacing_type
    
      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(:)
      real(kind=ireals), pointer :: dzetai_int(:), dzetai_05int(:) 
    
      integer(kind=iintegers), pointer :: isoilcolc(:), ksoilcol(:)
    
    type(gridspacing_type) :: gsp
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    MODULE ARRAYS_METHANE
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !Methane group
     real(kind=ireals) :: veg
    
     real(kind=ireals) :: fplant, fdiff_lake_surf, foutflow
    
     real(kind=ireals), allocatable :: fdiffbot(:)
     real(kind=ireals), allocatable :: febul0(:)
    ! real(kind=ireals) :: febul2(1:2) ! two-meth
    ! real(kind=ireals) :: fdiff2(1:2), fdiff_lake_surf2(1:2) ! two-meth
     real(kind=ireals), target :: febultot(1:2), febulbottot(1:2)
     real(kind=ireals) :: fdifftot(1:2), fdiff_lake_surftot(1:2)
     real(kind=ireals) :: methoxidwat(1:2), ice_meth_oxid
    
     real(kind=ireals) :: qwateroxidtot, qwateroxidML
    
     real(kind=ireals) :: ice_meth_oxid_total
     real(kind=ireals) :: meth_cont_wat_total0
    ! real(kind=ireals) :: febultot2(1:2,1:2), fdifftot2(1:2,1:2), fdiff_lake_surftot2(1:2,1:2) ! two-meth
     real(kind=ireals) :: plant_sum, bull_sum, oxid_sum, rprod_sum
     real(kind=ireals), allocatable :: rprod_total_oldC(:), rprod_total_newC(:)
     real(kind=ireals) :: rprod_total_newC_integr(1:2), rprod_total_oldC_integr(1:2)
    ! real(kind=ireals) :: rprod_total_newC_integr2(1:2,1:2), rprod_total_oldC_integr2(1:2,1:2) ! two-meth
     real(kind=ireals) :: h_talik
     real(kind=ireals) :: anox
     real(kind=ireals), allocatable :: rootss(:), TgrAnn(:)
     real(kind=ireals), allocatable, target :: qsoil(:,:), qwater(:,:), qmethane(:)
     real(kind=ireals), allocatable :: fbbleflx_ch4(:,:), fbbleflx_ch4_sum(:)
     real(kind=ireals), allocatable :: fracb0(:)
     real(kind=ireals), allocatable :: methgenmh(:) 
    ! real(kind=ireals), allocatable :: qsoil2(:,:), qwater2(:,:,:)  ! two-meth
     real(kind=ireals), allocatable :: lammeth(:)
     real(kind=ireals) :: tot_ice_meth_bubbles
    ! real(kind=ireals) :: tot_ice_meth_bubbles2(1:2)  ! two-meth
    
    END MODULE ARRAYS_METHANE
    
    
    MODULE ARRAYS_OXYGEN
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !Oxygen group      
    real(kind=ireals), allocatable, target :: oxyg(:,:)
    real(kind=ireals), allocatable :: lamoxyg(:)
    real(kind=ireals), allocatable :: fbbleflx_o2(:,:), fbbleflx_o2_sum(:)
    real(kind=ireals), allocatable :: Chl_a(:) ! Chlorophyll-a concentration, mg/l
    real(kind=ireals), allocatable :: prodox(:) ! O_2 production due to photosynthesis
    real(kind=ireals), allocatable :: resp(:) ! O_2 sink by respiration
    real(kind=ireals), allocatable :: bod(:) ! Biochemical oxygen demand 
    real(kind=ireals), allocatable :: sod(:) ! Sedimentary oxygen demand (at the lake margins)
    real(kind=ireals) :: sodbot ! Sedimentary oxygen demand at the lake bottom
    
    real(kind=ireals), target :: H_photic_zone ! The photic zone depth, m
    
    real(kind=ireals) :: fo2 ! Diffusive flux of oxygen to the atmosphere, mol/(m**2*s)
    integer(kind=iintegers) :: itroph ! Trophic status of a lake
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    real(kind=ireals), allocatable :: oxygsoil(:) ! Mean oxygen concentration in the uppermost (aerobic)
    
                                                  ! soil layer
    
    real(kind=ireals), allocatable, target :: DIC(:,:)
    
    real(kind=ireals), allocatable, target :: DOC(:,:) !DOC(:,1) stands for autochtonous DOC, DOC(:,2) -- for allochtonous one
    real(kind=ireals), allocatable, target :: POCL(:), POCD(:)
    
    real(kind=ireals), allocatable, target :: lamcarbdi(:), lampart(:)
    
    real(kind=ireals) :: fco2 ! Diffusive flux of carbon dioxide to the atmosphere, mol/(m**2*s)
    
    real(kind=ireals), allocatable :: fbbleflx_co2(:,:), fbbleflx_co2_sum(:)
    
    
    !Phosphorus group
    real(kind=ireals), allocatable, target :: phosph(:)
    
    
    END MODULE ARRAYS_OXYGEN
    
    
    
    MODULE ARRAYS_TURB
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    real(kind=ireals), allocatable :: S(:),Gen(:),F(:),TKE_turb_trans(:),Gen_seiches(:)
    
    real(kind=ireals), allocatable :: KT(:),k2(:),E1(:),eps1(:)
    
    real(kind=ireals), allocatable :: KC(:),KLengT(:),RSR(:) 
    
    real(kind=ireals), allocatable :: E2(:),eps2(:), &
    
    & k1(:),k3(:),k4(:),k5(:),u3(:), &
    
    & v3(:),C1aup(:),Re(:),Ri(:),E_it1(:), &
    
    & E_it2(:),C_num(:),E_it3(:),Feps(:),E_it21(:), &
    & eps_it21(:),eps_it1(:),eps_it2(:), &
    & l(:),k2_mid(:),k3_mid(:),k4_mid(:),k5_mid(:), &
    & E12(:),eps12(:),knum(:),k2t(:), &
    & Eeps(:),Ecent(:),Epscent(:),res_E(:), &
    & res_eps(:),dresd(:,:,:,:),dres2dE(:),dres2deps(:), &
    & veg_sw(:) 
    
    real(kind=ireals), allocatable :: eps_ziriv(:), viscturb_ziriv(:), speed_ziriv(:)
    
    real(kind=ireals), allocatable, target :: row(:), row2(:), rowc(:)
    
    real(kind=ireals), allocatable :: WU_(:),WV_(:),GAMT(:),GAMU(:),GAMV(:), TF(:),KLengu(:)
    real(kind=ireals), allocatable :: PEMF    (:) , PDENS_DOWN (:), &
    &                       PT_DOWN (:) , PSAL_DOWN  (:), &
    &                       pt_down_f(:)
    real(kind=ireals), allocatable :: k_turb_T_flux(:), T_massflux(:)
    
    
    integer(kind=iintegers) :: i_maxN
    
    integer(kind=iintegers), allocatable :: itherm(:)
    
    
    real(kind=ireals) :: Buoyancy0,H_mixed_layer,maxN,w_conv_scale,u_star,T_conv_scale,H_entrainment, &
    
    & signwaveheight = 0.d0, Ri_bulk, Wedderburn, Rossby_rad, LakeNumber, ThermThick, ReTherm, RiTherm
    
    real(kind=ireals) :: S_integr_positive, S_integr_negative, Gen_integr, &
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    & eps_integr, E_integr, TKE_turb_trans_integr 
    
    real(kind=ireals) :: Seps_integr_positive, Seps_integr_negative, &
    & Geneps_integr, epseps_integr
    real(kind=ireals) :: TKE_balance, eps_balance
    
    real(kind=ireals) :: Eseiches
    
    real(kind=ireals) :: TKE_budget_terms(1:5)
    
    ! Data structure for turbulent characteristics
    type, public :: turb_type
      sequence
      real(kind=ireals), allocatable :: Rp(:), Rpdens(:)
      real(kind=ireals), pointer :: Ri(:)
      real(kind=ireals), pointer :: Buoyancy0, H_mixed_layer, maxN, w_conv_scale, &
      & u_star, T_conv_scale, H_entrainment, signwaveheight, Ri_bulk, Wedderburn, &
      & Rossby_rad, LakeNumber, ThermThick, ReTherm, RiTherm
    end type
    type(turb_type) :: trb
    
    
    
    END MODULE ARRAYS_TURB
    
    
    
    MODULE ARRAYS_BATHYM
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    implicit none
    
    !Bathymetry group
    type, public :: bathym
      sequence
    
      real(kind=ireals) :: area_int, area_half, Lx, Ly, Lx_half, Ly_half, dzSL, dzSLc, rad_int
    
      real(kind=ireals), allocatable :: sarea_int(:), sarea_half(:)
    
      integer(kind=iintegers) :: itop, ibot, icent
    
    end type
    type(bathym), allocatable :: bathymwater(:), bathymice(:) 
    type(bathym), allocatable :: bathymdice(:), bathymsoil(:,:,:)
    real(kind=ireals), allocatable :: area_int(:), area_half(:), Lx(:), Ly(:)
    
    real(kind=ireals), allocatable :: aki(:,:), bki(:,:)
    
    real(kind=ireals), target :: h1, l1, hs1, ls1, vol, botar
    real(kind=ireals), target :: hx1, hx2, hy1, hy2, hx1t, hx2t, hy1t, hy2t
    
    real(kind=ireals), allocatable, target :: hx1ml(:), hx2ml(:), hy1ml(:), hy2ml(:)
    
    type, public :: layers_type
      sequence
      real(kind=ireals), pointer :: h1, hx1, hx2, hy1, hy2, l1, hs1, ls1, vol, &
    
      & dhw, dhw0, dhi, dhi0, dls, dls0, dhiimp, snmeltice, botar, H_photic_zone, &
    
      & hx1t, hx2t, hy1t, hy2t
    
      real(kind=ireals), pointer :: dhwhigh, dhwlow, dhwp, dhwe, dhwtrib, dhwls, dhwsnmelt
    
      real(kind=ireals), pointer :: dhihigh, dhilow, dhip, dhis, dhif
    
      real(kind=ireals), pointer :: hx1ml(:), hx2ml(:), hy1ml(:), hy2ml(:)
    
    end type
    type(layers_type) :: ls
    
    real(kind=ireals) :: voldef = 0.d0
    
    real(kind=ireals), target :: dhw,dhw0,dhi,dhi0,dls,dls0,dhiimp,snmeltice
    
    real(kind=ireals), target :: dhwhigh, dhwlow, dhwp, dhwe, dhwtrib, dhwls, dhwsnmelt
    
    real(kind=ireals), target :: dhihigh, dhilow, dhip, dhis, dhif
    
    
    real(kind=ireals) :: dhwfsoil = 0.
    
    
    integer(kind=iintegers), parameter :: form_ellipse = 1, form_rectangle = 2
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    use LAKE_DATATYPES, only : ireals, iintegers
    !use DRIVING_PARAMS
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     ! Soil group
     real(kind=ireals), allocatable:: rosoil(:), csoil(:), &
     & lamsoil(:), dzs(:), dzss(:), zsoil(:), wsoil(:)
     real(kind=ireals), allocatable :: lammoist(:), rosdry(:), filtr(:), wa(:)
     real(kind=ireals), allocatable :: Tsoil1(:,:), Tsoil2(:), Tsoil3(:,:)
     real(kind=ireals), allocatable :: wl1(:,:), wl2(:), wl3(:), wl4(:,:)
     real(kind=ireals), allocatable :: wi1(:,:), wi2(:,:)
     real(kind=ireals), allocatable :: Sals1(:,:), Sals2(:,:)
     real(kind=ireals), allocatable :: soilflux(:,:,:)
     real(kind=ireals), allocatable :: methsoilflux(:), co2soilflux(:), o2soilflux(:)
     type, public :: lsh_type
       sequence
       real(kind=ireals), allocatable :: water(:), ice(:), dice(:)
     end type lsh_type
    
     type(lsh_type) :: lsh, lsm, lsc, lso, lsp, lsu, lsv
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     real(kind=ireals), allocatable :: WLM0(:),WLM7(:),bH(:),PSIMAX(:),POR(:),FLWMAX(:),DLMAX(:)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     MODULE ARRAYS
     
    !MODULE ARRAYS contains allocatable arrays of the model  
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     use LAKE_DATATYPES, only : ireals, iintegers
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
     integer(kind=iintegers), parameter :: soil_indic = 1
     integer(kind=iintegers), parameter :: ice_indic = 2
    
     integer(kind=iintegers), parameter :: ice_sal_indic = 12
     integer(kind=iintegers), parameter :: snow_indic = 3
     integer(kind=iintegers), parameter :: water_indic = 4
     integer(kind=iintegers), parameter :: deepice_indic = 5
     integer(kind=iintegers), parameter :: water_salinity_indic = 6
     integer(kind=iintegers), parameter :: soil_salinity_indic = 7
     integer(kind=iintegers), parameter :: water_methane_indic = 8
     integer(kind=iintegers), parameter :: water_oxygen_indic = 9
     integer(kind=iintegers), parameter :: water_carbdi_indic = 10
     integer(kind=iintegers), parameter :: water_carbon_indic = 11
    
     integer(kind=iintegers), parameter :: water_particles_indic = 13
    
     integer(kind=iintegers) :: snow,ice,water,deepice,nstep = 0
    
     real(kind=ireals) :: time, Erad, dep_av
     real(kind=ireals) :: zinv
     integer(kind=iintegers) :: time_int,time_int_old,man,par,flag_snow_init
    
    
    !Computational time group
    
     integer(kind=iintegers), parameter :: ncomp = 10
     real(kind=ireals) :: comptime(1:ncomp) = 0.d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
    !Parameter calibration group      
    
     real(kind=ireals), pointer, save :: calibr_par1(:), calibr_par2, cost_function(:)
    
     logical, save :: calibrate_parameter = .false.
    
       real(kind=ireals), pointer :: qmethane(:),qsoil(:,:),qwater(:,:), &
    
       & oxyg(:,:),DIC(:,:),DOC(:,:),POCL(:),POCD(:),phosph(:)
    
    !Surface characteristics group
    
     real(kind=ireals) :: roughness,emissivity,albedo,albedo_lw,aM,bM,relhums
    
     real(kind=ireals) totalevaps,totalmelts,totalprecips
     real(kind=ireals), allocatable :: Tskin(:)
     real(kind=ireals) :: T_0dim
    
     real(kind=ireals) :: SR_botsnow
     real(kind=ireals) :: dt_keps
     real(kind=ireals) :: dt05, dt_inv, dt_inv05
     real(kind=ireals) :: deltaskin
    
    
     real(kind=ireals), allocatable :: z_watersoil(:)
    
     real(kind=ireals), allocatable, target :: u1(:),v1(:),u2(:),v2(:),um(:),vm(:),uv(:)
    
     real(kind=ireals), allocatable, target :: wArea(:),w(:)
    
     real(kind=ireals), allocatable :: dep_2d(:,:)
    
     integer(kind=iintegers), allocatable :: init(:,:), num(:)
    
     real(kind=ireals), allocatable :: utend(:), vtend(:), Twtend(:), Saltend(:)
    
     real(kind=ireals), allocatable :: gradpx_tend(:), gradpy_tend(:)
    
     real(kind=ireals), allocatable :: ueffl(:), veffl(:), Tweffl(:), Saleffl(:)
    
     real(kind=ireals), pointer :: workc(:)
    
    
     MODULE TURB_CONST
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !PHYSICAL CONSTANTS
    
      !real(kind=ireals), parameter:: niu       = 1.007d-6
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !DIMENSIONLESS CONSTANTS IN E-EPS PARAMETERIZATION, according to Goudsmit et al., 2002
    
      real(kind=ireals), parameter:: CE0       = 0.09d0 
      real(kind=ireals), parameter:: CEt0      = 0.072d0
      real(kind=ireals), parameter:: ceps1     = 1.44d0
      real(kind=ireals), parameter:: ceps2     = 1.92d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
      real(kind=ireals), parameter:: ceps3_unstable = 1.14d0
      real(kind=ireals), parameter:: ceps3_stable   = -0.4d0
    !        ceps3 = 1.14d0
    !         ceps3 < -0.4 should be used for stable stratification (Burchard, 2002)
    !         ceps3 = 1. ! following (Kochergin and Sklyar, 1992)
    !         ceps3 = 1.14d0 ! Baum and Capony (1992)
    !         Other values for ceps3:
    !         0.<ceps3<0.29 for stable and ceps3 = 1.44 for unstable conditions (Rodi, 1987);
    !         ceps3 >= 1 (Kochergin and Sklyar, 1992)
    
      real(kind=ireals), parameter:: sigmaE    = 1.d0
    
    !  real(kind=ireals), parameter:: sigmaeps  = 1.3d0
    
      real(kind=ireals), save :: CL 
      real(kind=ireals), save :: alpham
    
      real(kind=ireals), parameter :: C_wstar = 0.4d0   ! Deardorff, 1970
    
      real(kind=ireals), parameter:: lam_T     = 1.d0 
      real(kind=ireals), parameter:: lam_gen   = 1.d0 
      real(kind=ireals), parameter:: lamTM     = 0.14d0
    
    ! real(kind=ireals), parameter:: cmu_0     = 0.094d0
    
    ! Other values for cmu_0, proposed in literature, are
    ! 0.121, 0.077 - see Burchard, 2002
    
    
      real(kind=ireals), parameter:: lam_E0     = CE0/sigmaE   
      real(kind=ireals), save :: lam_eps0
    
    
    ! The coefficients in formulas for stability functions (Canuto et al., 2001)
    
      real(kind=ireals), parameter:: CEcoef01   = 0.127d0
      real(kind=ireals), parameter:: CEcoef02   = 0.1526d-1
      real(kind=ireals), parameter:: CEcoef03   = 0.16d-3
    
      real(kind=ireals), parameter:: CEtcoef01  = 0.1190d0
      real(kind=ireals), parameter:: CEtcoef02  = 0.4294d-2
      real(kind=ireals), parameter:: CEtcoef03  = 0.66d-3
    
      real(kind=ireals), parameter:: CEcoef11   = 1.d0
      real(kind=ireals), parameter:: CEcoef12   = 0.2d0
      real(kind=ireals), parameter:: CEcoef13   = 0.315d-1
      real(kind=ireals), parameter:: CEcoef14   = 0.58d-2
      real(kind=ireals), parameter:: CEcoef15   = 0.4d-2
      real(kind=ireals), parameter:: CEcoef16   = 0.4d-4
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
     
    
    ! The parameters of Mellor and Yamada model (1982)      
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), parameter:: A1  =  0.92d0
      real(kind=ireals), parameter:: B1  =  1.66d+1
      real(kind=ireals), parameter:: C1  =  0.8d-1
      real(kind=ireals), parameter:: A2  =  0.74d0
      real(kind=ireals), parameter:: B2  =  1.01d+1
      real(kind=ireals), save :: CL_K_KL_MODEL
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !Co is according to Satyanarayana et al., 1999
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !CONSTANTS IN CONVECTION PARAMETERIZATION
    
      real(kind=ireals), parameter:: KC0       = 0.d0 !1.d+5 !500.
      real(kind=ireals), parameter:: dtdz0     = 2.d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    !CONSTANTS OF SIMOES'S PARAMETERIZATION
    
      real(kind=ireals), parameter:: Cs        = 1.5d0  !0.15
      real(kind=ireals), parameter:: Cs1       = 100.d0    
    
     
    !lo=0.4*ddz*h1*3.                                                         VS,06.2007
    
      real(kind=ireals), parameter:: L0        = 0.1d0 !0.1d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), parameter:: CON0      = 0.05d0
      real(kind=ireals), parameter:: CON1      = 0.11d0
      real(kind=ireals), parameter:: CON2      = 1.56d0
      real(kind=ireals), parameter:: CONUV     = 1.d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), parameter :: kar_tilde = 2.d-1 ! Burchard (2002), p.80
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      real(kind=ireals), parameter :: grav_wave_Gill_const = 0.7 !Burchard(2002), Gill(1982)
    
    
      !real(kind=ireals), parameter :: min_turb = 1.d-10 !Minimal turbulent viscosity/diffusivity coefficient, m**2/s
    
      real(kind=ireals), parameter :: min_diff = 1.d-10 !2.d-7 !Minimal turbulent diffusivity coefficient, m**2/s
    
      real(kind=ireals), parameter :: Pr_stable = 1. !1.E+1 !Turbulent Prandtl number at very stable stratification
    
      real(kind=ireals), parameter :: min_visc = min_diff*Pr_stable
    
      real(kind=ireals), parameter :: CDN10m = 1.3E-3 !momentum exchange coefficient in the surface air at 10 m
      real(kind=ireals), save :: CDN !momentum exchange coefficient in the surface air
      real(kind=ireals), parameter :: z_CDN = 2.5 !Height of wind measurements, for CDN calculation only
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      SUBROUTINE TURB_CONST_DERIVED_SET
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      use DRIVING_PARAMS, only : kepsbc
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      CL = CE0**(0.75d0) ! This value results from assumption
                         ! that shear production is balanced by dissipation
                         ! (Burchard, 2002)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      CL_K_KL_MODEL = 2.d0**(1.5d0)/B1
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      alpham = (1.5*CE0**0.5*sigmaE*kar_tilde**(-2))**0.5
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      select case(kepsbc%par)
        case(1)
          sigmaeps  = kar*kar/(CE0**0.5 * (ceps2 - ceps1) ) ! according to Burchard, 2002,
                                                            ! fulfills the law of the wall
                                                            ! in k-eps model
        case(2)
          sigmaeps = (4./3.*alpham + 1.)*(alpham + 1.) * &
          & kar_tilde*kar_tilde/(ceps2*CE0**0.5)  ! according to Burchard, 2002,
                                                  ! fulfills the analytical solution for 
                                                  ! unstratified shear-free flow with 
                                                  ! wave breaking TKE input
        case default
          sigmaeps  = 1.3d0
      end select
    !  sigmaeps = 1.3d0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      lam_eps0   = CE0/sigmaeps 
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
      CDN = CDN10m/(1-CDN10m**0.5/kar*log(10.0/z_CDN))**2.0
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     END SUBROUTINE TURB_CONST_DERIVED_SET
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
     END MODULE TURB_CONST
    
    !> Module EVOLUTION_VARIABLES contains variables
    !! that have to be saved for the next time step (~prognostic variables)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    MODULE EVOLUTION_VARIABLES
    
    
    real(kind=ireals), allocatable, public :: l1_2d(:,:)
    real(kind=ireals), allocatable, public :: h1_2d(:,:)
    real(kind=ireals), allocatable, public :: hx1_2d(:,:), hx2_2d(:,:)
    real(kind=ireals), allocatable, public :: hy1_2d(:,:), hy2_2d(:,:)
    real(kind=ireals), allocatable, public :: hx1t_2d(:,:), hx2t_2d(:,:)
    real(kind=ireals), allocatable, public :: hy1t_2d(:,:), hy2t_2d(:,:)
    real(kind=ireals), allocatable, public :: hx1ml_2d(:,:,:), hx2ml_2d(:,:,:)
    real(kind=ireals), allocatable, public :: hy1ml_2d(:,:,:), hy2ml_2d(:,:,:)
    
    real(kind=ireals), allocatable, public :: gradpx_tend_2d(:,:,:), gradpy_tend_2d(:,:,:)
    
    real(kind=ireals), allocatable, public :: ls1_2d(:,:)
    real(kind=ireals), allocatable, public :: hs1_2d(:,:)
    real(kind=ireals), allocatable, public :: time_2d(:,:) 
    real(kind=ireals), allocatable, public :: cdm_2d(:,:)
    real(kind=ireals), allocatable, public :: u_2d(:,:,:)
    real(kind=ireals), allocatable, public :: v_2d(:,:,:)
    real(kind=ireals), allocatable, public :: Tsoil1_2d(:,:,:,:)
    real(kind=ireals), allocatable, public :: wi1_2d(:,:,:,:)
    real(kind=ireals), allocatable, public :: wl1_2d(:,:,:,:)
    real(kind=ireals), allocatable, public :: Tw1_2d(:,:,:)
    real(kind=ireals), allocatable, public :: Ti1_2d(:,:,:)
    real(kind=ireals), allocatable, public :: Tis1_2d(:,:,:)
    real(kind=ireals), allocatable, public :: dz_2d(:,:,:)
    real(kind=ireals), allocatable, public :: T_2d(:,:,:)
    real(kind=ireals), allocatable, public :: wl_2d(:,:,:)
    real(kind=ireals), allocatable, public :: dens_2d(:,:,:)
    real(kind=ireals), allocatable, public :: E_2d(:,:,:)
    real(kind=ireals), allocatable, public :: eps_2d(:,:,:)
    real(kind=ireals), allocatable, public :: dhwfsoil_2d(:,:) 
    real(kind=ireals), allocatable, public :: Sal1_2d(:,:,:)
    real(kind=ireals), allocatable, public :: Sals1_2d(:,:,:,:)
    
    real(kind=ireals), allocatable, public :: ueffl_2d(:,:,:) !>@deprecated
    real(kind=ireals), allocatable, public :: veffl_2d(:,:,:)!>@deprecated
    real(kind=ireals), allocatable, public :: Tweffl_2d(:,:,:)!>@deprecated
    real(kind=ireals), allocatable, public :: Saleffl_2d(:,:,:)!>@deprecated
    
    real(kind=ireals), allocatable, public :: Elatent_2d(:,:)
    real(kind=ireals), allocatable, public :: dhw_2d(:,:)
    real(kind=ireals), allocatable, public :: dhw0_2d(:,:)
    real(kind=ireals), allocatable, public :: dhi_2d(:,:)
    real(kind=ireals), allocatable, public :: dhi0_2d(:,:)
    real(kind=ireals), allocatable, public :: dls0_2d(:,:)
    real(kind=ireals), allocatable, public :: velfrict_2d(:,:)
    real(kind=ireals), allocatable, public :: roughness_2d(:,:)
    real(kind=ireals), allocatable, public :: eflux0_kinem_2d(:,:)
    real(kind=ireals), allocatable, public :: lamw_2d(:,:,:)
    real(kind=ireals), allocatable, public :: snmelt_2d(:,:)
    
    real(kind=ireals), allocatable, public :: snowmass_2d(:,:)
    
    real(kind=ireals), allocatable, public :: Tskin_2d(:,:)
    real(kind=ireals), allocatable, public :: k_turb_T_flux_2d(:,:,:)
    real(kind=ireals), allocatable, public :: qwater_2d(:,:,:)
    real(kind=ireals), allocatable, public :: qsoil_2d(:,:,:,:)
    real(kind=ireals), allocatable, public :: oxyg_2d(:,:,:)
    real(kind=ireals), allocatable, public :: oxygsoil_2d(:,:,:)
    real(kind=ireals), allocatable, public :: carbdi_2d(:,:,:)
    
    real(kind=ireals), allocatable, public :: phosph_2d(:,:,:)
    
    real(kind=ireals), allocatable, public :: DOC_2d(:,:,:,:)
    
    real(kind=ireals), allocatable, public :: POCL_2d(:,:,:)
    real(kind=ireals), allocatable, public :: POCD_2d(:,:,:)
    
    real(kind=ireals), allocatable, public :: tot_ice_meth_bubbles_2d(:,:)
    real(kind=ireals), allocatable, public :: febul0_2d(:,:,:)
    real(kind=ireals), allocatable, public :: Eseiches_2d(:,:)
    
    real(kind=ireals), allocatable, public :: salice_2d(:,:,:)
    
    real(kind=ireals), allocatable, public :: porice_2d(:,:,:)
    
    
    integer(kind=iintegers), allocatable, public :: fl_sn_2d(:,:)
    integer(kind=iintegers), allocatable, public :: fl_sn_init_2d(:,:) 
    integer(kind=iintegers), allocatable, public :: itop_2d(:,:)
    integer(kind=iintegers), allocatable, public :: nstep_2d(:,:)
    integer(kind=iintegers), allocatable, public :: i_maxN_2d(:,:)
    
    integer(kind=iintegers), allocatable, public :: itherm_2d(:,:,:)
    
    
    logical, public :: arrays_allocated = .false.
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    SAVE
    
    contains
    
    !> Subroutine ALLOCATE_ARRAYS allocates arrays of EVOLUTION_VARIABLES module
    
    SUBROUTINE ALLOCATE_ARRAYS(nx, ny, M, Mice, ns, ms, ml, nsoilcols)
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    implicit none
    
    ! Input variables
    
     integer(kind=iintegers), intent(in) :: nx, ny
     integer(kind=iintegers), intent(in) :: M, Mice, ns, ms, ml, nsoilcols
    
    Victor Stepanenko's avatar
    Victor Stepanenko committed
    
     allocate ( l1_2d(nx,ny) ) ; l1_2d = 0
     allocate ( h1_2d(nx,ny) ) ; h1_2d = 0