diff --git a/source/model/init.f90 b/source/model/init.f90
index c0047797e4489dcda3e0701ab5513e418789d2c0..9caf9b8768caa3c73ed71e1b0abc6aaf0048e7f9 100644
--- a/source/model/init.f90
+++ b/source/model/init.f90
@@ -222,6 +222,10 @@ END MODULE GRID_OPERATIONS_LAKE
allocate (ci_m_roi_v(1:Mice+1), lami_v(1:Mice))
allocate ( Tskin(1:2) )
+ !Time group
+ time_group%time => time
+ time_group%nstep => nstep
+
! Radiation group
!allocate (fracbands(1:nbands))
call RAD_POINTERS()
diff --git a/source/model/lake.f90 b/source/model/lake.f90
index ba60c978da654990a2fb5deeb6316d735e6e4c1e..2e5e9d8aa0ab81b03588e4fc1a78d55d7d22a1da 100644
--- a/source/model/lake.f90
+++ b/source/model/lake.f90
@@ -222,7 +222,9 @@ use PHYS_FUNC, only: &
& WL_MAX, &
& MIXED_LAYER_CALC, &
& HC_CORR_CARBEQUIL, &
-& DIFFMIN_HS
+& DIFFMIN_HS, &
+& LONGWAVE_RADIATION, &
+& SHORTWAVE_RADIATION
use EVOLUTION_VARIABLES, only : &
& UPDATE_CURRENT_TIMESTEP, &
@@ -286,35 +288,23 @@ implicit none
!Input variables
!Reals
-!> Air temperature, K
-real(kind=ireals), intent(in) :: tempair1
-!> Specific humidity, kg/kg
-real(kind=ireals), intent(in) :: humair1
-!> Air pressure, Pa
-real(kind=ireals), intent(in) :: pressure1
-!> x- and y-components of wind, m/s
-real(kind=ireals), intent(in) :: uwind1, vwind1
-!> Longwave downward radiation, W/m**2
-real(kind=ireals), intent(in) :: longwave1
-!> Net radiation, W/m**2
-real(kind=ireals), intent(in) :: netrad1
-!> Height of level above the lake surface, where the forcing is given, m
-real(kind=ireals), intent(in) :: zref1
-!> Global shortwave radiation, W/m**2
-real(kind=ireals), intent(in) :: shortwave1
-!> Precipitation intensity, m/s
-real(kind=ireals), intent(in) :: precip1
-real(kind=ireals), intent(in) :: Sflux01
-!> Cloudiness, fraction
-real(kind=ireals), intent(in) :: cloud1
-!> Surface temperature, K
-real(kind=ireals), intent(in) :: SurfTemp1
+real(kind=ireals), intent(in) :: tempair1 !> Air temperature, K
+real(kind=ireals), intent(in) :: humair1 !> Specific humidity, kg/kg
+real(kind=ireals), intent(in) :: pressure1 !> Air pressure, Pa
+real(kind=ireals), intent(in) :: uwind1, vwind1 !> x- and y-components of wind speed, m/s
+real(kind=ireals), intent(in) :: longwave1 !> Longwave downward radiation, W/m**2
+real(kind=ireals), intent(in) :: netrad1 !> Net radiation, W/m**2
+real(kind=ireals), intent(in) :: zref1 !> Height above the lake surface, where the forcing is given, m
+real(kind=ireals), intent(in) :: shortwave1 !> Global shortwave radiation, W/m**2
+real(kind=ireals), intent(in) :: precip1 !> Precipitation intensity, m/s
+real(kind=ireals), intent(in) :: Sflux01 !> Salinity flux at the water surface by rain
+real(kind=ireals), intent(in) :: cloud1 !> Cloud sky fraction, 0-1, n/d
+real(kind=ireals), intent(in) :: SurfTemp1 !> Surface temperature, K
real(kind=ireals), intent(in) :: ch4_pres_atm, co2_pres_atm, o2_pres_atm
-!> Time step for the LAKE model, s
-real(kind=ireals), intent(in) :: dtl
+real(kind=ireals), intent(in) :: dtl !> Time step for the LAKE model, s
-real(kind=ireals), intent(in) :: hour
+real(kind=ireals), intent(in), target :: hour
real(kind=ireals), intent(in) :: h10, l10, ls10, hs10
real(kind=ireals), intent(in) :: Ts0, Tb0, Tbb0
@@ -322,8 +312,8 @@ real(kind=ireals), intent(in) :: h_ML0
real(kind=ireals), intent(in) :: extwat, extice
real(kind=ireals), intent(in) :: kor_in
-!> Tributary inflow discharge, m**3/s
-real(kind=ireals), intent(in) :: trib_inflow
+
+real(kind=ireals), intent(in) :: trib_inflow !> Tributary inflow discharge, m**3/s
real(kind=ireals), intent(in) :: Sals0, Salb0
real(kind=ireals), intent(in) :: fetch
real(kind=ireals), intent(in) :: phi, lam
@@ -348,7 +338,7 @@ integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: nx0, ny0, nx_max, ny_max
integer(kind=iintegers), intent(in) :: ndatamax
-integer(kind=iintegers), intent(in) :: year, month, day
+integer(kind=iintegers), intent(in), target :: year, month, day
integer(kind=iintegers), intent(in) :: init_T
!> MPI parameters
integer(kind=iintegers), intent(in) :: comm3d, rank_comm3d, coords(1:3)
@@ -487,6 +477,11 @@ parparams%parallel = parallel
!call TIMEC(0)
+time_group%year => year
+time_group%month => month
+time_group%day => day
+time_group%hour => hour
+
gs%ix = ix
gs%iy = iy
@@ -585,15 +580,28 @@ else
netrad = netrad1
endif
+!Incident global shortwave radiation is computed by empirical formulae, if missing value is provided
+!in the LAKE arguments
+if (shortwave1 == missing_value) then
+ if (cloud1 /= missing_value) then
+ shortwave = SHORTWAVE_RADIATION(time_group,phi,cloud/10.)
+ else
+ print*, 'Global solar radiaion and cloud fraction are not provided: STOP'
+ STOP
+ endif
+endif
-if (longwave1 == missing_value .and. firstcall) then
- write(*,*) 'Atmospheric radiation is missing in input file and will be calculated &
- &from net longwave radiation or net radiation'
+if (longwave1 == missing_value) then
+ if (firstcall) write(*,*) 'Atmospheric radiation is missing in input file and will be calculated &
+ &from net longwave radiation, net radiation or empirical formulae'
if (cloud1 == missing_value) then
- write(*,*) 'Cloudiness data is missing: atmospheric radation to be computed from net radiation'
+ if (firstcall) write(*,*) 'Cloudiness data is missing: atmospheric radation to be computed &
+ &from net radiation'
if (netrad1 == missing_value) then
STOP 'Net radiation data is missing!: STOP'
endif
+ else
+ longwave = LONGWAVE_RADIATION(tempair,humair,pressure,cloud/10.)
endif
endif
@@ -2910,7 +2918,7 @@ endif
deallocate(radflux)
deallocate(a,b,c,d)
deallocate(Temp)
-
+nullify(time_group%year,time_group%month,time_group%day,time_group%hour)
if (firstcall) firstcall = .false.
diff --git a/source/model/lake_modules.f90 b/source/model/lake_modules.f90
index a63a5265b1eabe08c3ff5c1f316a955a2d436863..b87268d0cb841be24a406ab24f1f2d33a07a883c 100644
--- a/source/model/lake_modules.f90
+++ b/source/model/lake_modules.f90
@@ -20,6 +20,7 @@ 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
+real(kind=ireals), parameter :: Pa_to_hPa = 1.E-2
END MODULE UNITS
@@ -150,6 +151,7 @@ real(kind=ireals) :: Planck_const, Planck_const_
real(kind=ireals) :: clight
real(kind=ireals) :: lambda_PAR0
real(kind=ireals) :: short2PAR
+real(kind=ireals) :: solar_constant
real(kind=ireals) :: lami, lamw0, lamair
real(kind=ireals) :: lamsal0
real(kind=ireals) :: niu_atm, niu_wat
@@ -227,6 +229,7 @@ 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
+solar_constant = 1367. ! the solar constant (mean extraterrestrial solar radiation), W/m**2
aMagw = 7.6326 ! coefficient for Magnus formula for water surface
bMagw = 241.9 ! coefficient for Magnus formula for water surface
@@ -348,7 +351,6 @@ SAVE
END MODULE ATMOS
-
!> Radiation data structures and operations on them
MODULE RADIATION
@@ -422,6 +424,7 @@ END SUBROUTINE RAD_POINTERS
END MODULE RADIATION
+
MODULE ARRAYS_WATERSTATE
use LAKE_DATATYPES, only : ireals, iintegers
@@ -716,9 +719,17 @@ END MODULE ARRAYS_SOIL
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
+ type, public :: time_type
+ sequence
+ real (kind=ireals) , pointer :: time, hour
+ integer(kind=iintegers), pointer :: year, month, day, nstep
+ end type time_type
+ type(time_type) :: time_group
+
+ integer(kind=iintegers) :: snow,ice,water,deepice
+ integer(kind=iintegers), target :: nstep = 0
+ real(kind=ireals) :: Erad, dep_av
+ real(kind=ireals), target :: time
real(kind=ireals) :: zinv
integer(kind=iintegers) :: time_int,time_int_old,man,par,flag_snow_init
diff --git a/source/model/phys_func.f90 b/source/model/phys_func.f90
index f8035e73c3b54a1dc3bb0274d3a8c68c14c75819..eab6a6e5a0c3ca02f22cf7da95395adf69f74c54 100644
--- a/source/model/phys_func.f90
+++ b/source/model/phys_func.f90
@@ -610,6 +610,107 @@ END FUNCTION GSW_RHO
END FUNCTION SINH0
+ !> Function LONGWAVE_RADIATION returns incident
+ !! longwave (atmospheric) radiation on the earth surface
+ FUNCTION LONGWAVE_RADIATION(tempair,humair,pressure,cloud_sky_fraction) result(LR)
+ use UNITS, only : Pa_to_hPa
+ use PHYS_CONSTANTS, only : Kelvin0, Rd_d_Rwv, sigma
+ implicit none
+ real(kind=ireals), intent(in) :: tempair !> Surface air temperature, Celsius
+ real(kind=ireals), intent(in) :: humair !> Surface air specific humidity, kg/kg
+ real(kind=ireals), intent(in) :: pressure!> Surface atmospheric pressure, Pa
+ real(kind=ireals), intent(in) :: cloud_sky_fraction !> Cloud-sky fraction, 0-1, n/d
+ !> Output variables
+ real(kind=ireals) :: LR
+ !> Local variables
+ real(kind=ireals) :: emis, a1, a2, a3, parpres, Rwv_d_Rd
+
+ integer(kind=iintegers), parameter :: CAN = 1 !(Anstrom, 1915; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: CBR = 2 !(Brunt, 1932; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: CSW = 3 !(Swinbank, 1963; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: CIJ = 4 !(Idso and Jackson, 1969; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: CBT = 5 !(Brutsaert, 1975, 1982; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: CID = 6 !(Idso, 1981; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: CKZ = 7 !(Konzelmann et al., 1994; Marthews et al., Theor Appl Climatol, 2012)
+ integer(kind=iintegers), parameter :: nCSE = CAN !switch for clear-sky emissivity parameterization
+
+ Rwv_d_Rd = 1./Rd_d_Rwv
+ parpres = humair*pressure*Rwv_d_Rd/(1. + humair*(Rwv_d_Rd - 1.))*Pa_to_hPa
+
+ select case(nCSE)
+ case(CAN)
+ a1 = 0.83
+ a2 = 0.18
+ a3 = 0.067
+ emis = a1 - a2*10**( - a3*parpres)
+ case(CBR)
+ a1 = 0.51
+ a2 = 0.066
+ emis = a1 + a2*sqrt(parpres)
+ case(CSW)
+ a1 = 0.0000092
+ emis = a1*(tempair + Kelvin0)**2
+ case(CIJ)
+ a1 = 0.261
+ a2 = 0.000777
+ emis = 1. - a1*exp( - a2*tempair**2)
+ case(CBT)
+ a1 = 1.24
+ emis = a1*(parpres/(tempair + Kelvin0))**(1./7.)
+ case(CID)
+ a1 = 0.7
+ a2 = 0.0000595
+ a3 = 1500.
+ emis = a1 + a2*parpres*exp(a3/(tempair + Kelvin0))
+ case(CKZ)
+ a1 = 0.23
+ a2 = 0.484
+ emis = a1 + a2*(parpres/(tempair + Kelvin0))**(1./8.)
+ end select
+
+ !Correction of emissivity for cloudiness
+ a1 = 0.22
+ emis = emis*(1. + a1*cloud_sky_fraction)
+ emis = min(emis,1._ireals)
+ LR = emis*sigma*(tempair + Kelvin0)**4
+
+ END FUNCTION LONGWAVE_RADIATION
+
+ !> Function SHORTWAVE_RADIATION returns incident
+ !! shortwave (solar) radiation (direct+diffuse) on the earth surface
+ FUNCTION SHORTWAVE_RADIATION(time_group,phi,cloud_sky_fraction) result(SR)
+ use PHYS_CONSTANTS, only : solar_constant
+ use ARRAYS, only : time_type
+ implicit none
+ !> Input variables
+ type(time_type), intent(in) :: time_group !>The group of time/date variables, time assumed local, not GMT
+ real(kind=ireals), intent(in) :: phi !>The latitude, deg.
+ real(kind=ireals), intent(in) :: cloud_sky_fraction !> Cloud-sky fraction, 0-1, n/d
+ !> Output variables
+ real(kind=ireals) :: SR
+ !> Local variables
+ real(kind=ireals) :: sinsolh, trans_atm, SR_dir, SR_dif, a1, a2, SR_atten
+
+ sinsolh = SINH0(time_group%year,time_group%month,time_group%day,time_group%hour,phi)
+ trans_atm = 0.6 ! atmospheric transmissivity for direct solar radiation
+ SR_atten = trans_atm**(1./sinsolh)
+ !Direct solar radiation at the horizontal surface (Burger's law)
+ SR_dir = solar_constant*SR_atten*sinsolh
+ !Diffuse solar radiation at the horizontal surface (Liu and Jordan, Solar Energy, 1960;
+ !Jemaa et al., Energy Procedia, 2013)
+ a1 = 0.271
+ a2 = 0.294
+ SR_dif = solar_constant*(a1 - a2*SR_atten)*sinsolh
+ SR = SR_dir + SR_dif
+ !Correction of global radiation by cloudiness (Kasten and Czeplak, Solar Energy, 1980;
+ !Jemaa et al., Energy Procedia, 2013)
+ a1 = 0.75
+ a2 = 3.4
+ SR = SR*(1.-a1*cloud_sky_fraction**a2)
+
+ END FUNCTION SHORTWAVE_RADIATION
+
+
real(kind=ireals) FUNCTION QS(phase,t,p)
! QS - specific humidity, kg/kg, for saturated water vapour
diff --git a/source/shared/time.f90 b/source/shared/time.f90
index decac8274dfd38e9c656f91f4c0a39efce13ecea..4308cb9040a0d2d9697b51c47a43eccece262dc3 100644
--- a/source/shared/time.f90
+++ b/source/shared/time.f90
@@ -28,8 +28,9 @@ data ndaym /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
SAVE ndaym
-if (mod(year,4_iintegers)==0) then
- ndaym(2) = 29 !Checking for leap-year
+if (mod(year,4 ) == 0 .and. .not. &
+& (mod(year,100) == 0 .and. mod(year,400) /= 0)) then ! Checking for leap-year
+ ndaym(2) = 29
else
ndaym(2) = 28
endif
@@ -68,7 +69,8 @@ data ndaym /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
SAVE ndaym
-if (mod(year,4)==0) then ! Checking for leap-year
+if (mod(year,4 ) == 0 .and. .not. &
+& (mod(year,100) == 0 .and. mod(year,400) /= 0)) then ! Checking for leap-year
ndaym(2) = 29
else
ndaym(2) = 28