Newer
Older
use LAKE_DATATYPES, only : ireals, iintegers
use PHYS_CONSTANTS, only : &
& row0
!The coefficients of UNESCO formula
!for water density dependence on temperature

Victor Stepanenko
committed
real(kind=ireals), parameter:: a0 = 800.969d-7
real(kind=ireals), parameter:: a1 = 588.194d-7
real(kind=ireals), parameter:: a2 = 811.465d-8
real(kind=ireals), parameter:: a3 = 476.600d-10
!The coefficients of formula for
!water density dependence on temperature and salinity
!(McCutcheon et al.,Water Quality and Maidment.
! Handbook on hydrology, 1993)

Victor Stepanenko
committed
real(kind=ireals), parameter:: A11 = 288.9414d0
real(kind=ireals), parameter:: A12 = 508929.2d0
real(kind=ireals), parameter:: A13 = 68.12963d0
real(kind=ireals), parameter:: A14 = 3.9863d0
real(kind=ireals), parameter:: alpha1 = 8.2449d-1
real(kind=ireals), parameter:: alpha2 = 4.0899d-3
real(kind=ireals), parameter:: alpha3 = 7.6438d-5
real(kind=ireals), parameter:: alpha4 = 8.2467d-7
real(kind=ireals), parameter:: alpha5 = 5.3675d-9
real(kind=ireals), parameter:: beta1 = 5.724d-3
real(kind=ireals), parameter:: beta2 = 1.0227d-4
real(kind=ireals), parameter:: beta3 = 1.6546d-6
real(kind=ireals), parameter:: gamma1 = 4.8314d-4

Victor Stepanenko
committed
real(kind=ireals), parameter :: at1 = 999.843, at2 = 65.4891d-3, at3 = - 8.56272d-3, at4 = 0.059385d-3
real(kind=ireals), parameter :: beta = 0.75d-3 ! kg/g
real(kind=ireals), parameter :: beta_co2 = 0.284d-3 ! coefficient for co2, kg/g
real(kind=ireals), parameter :: beta_ch4 = -1.25d-3 ! coefficient for ch4, kg/g
real(kind=ireals), parameter :: co2_ref = 0.d0, ch4_ref = 0.d0 ! reference concentrations, kg/g

Victor Stepanenko
committed
real(kind=ireals), save :: DDENS_DTEMP0, DDENS_DSAL0, DENS_REF

Victor Stepanenko
committed
integer(kind=iintegers), parameter :: eos_Hostetler = 1
integer(kind=iintegers), parameter :: eos_TEOS = 2
integer(kind=iintegers), parameter :: eos_Kivu = 3
integer(kind=iintegers), parameter :: eos_UNESCO = 4
integer(kind=iintegers), parameter :: eos_McCatcheon = 5

Victor Stepanenko
committed
!Temperature and salinity reference values
real(kind=ireals), parameter :: temp_ref = 15. ! deg. C
real(kind=ireals), parameter :: sal_ref = 0. !kg/kg
contains
SUBROUTINE SET_DENS_CONSTS(eos)
implicit none

Victor Stepanenko
committed
integer(kind=iintegers), intent(in) :: eos
select case (eos)
case(eos_Hostetler)
call DDENS_HOSTETLER(temp_ref,DDENS_DTEMP0,DDENS_DSAL0)

Victor Stepanenko
committed
DENS_REF = DENS_HOSTETLER(temp_ref)
case(eos_TEOS)
print*, 'Not operational eos: STOP '
STOP
case(eos_Kivu)
call DDENS_KIVU(temp_ref,sal_ref,DDENS_DTEMP0,DDENS_DSAL0)

Victor Stepanenko
committed
DENS_REF = WATER_DENS_TS_KIVU(temp_ref,sal_ref)
case(eos_UNESCO)
call DDENS_UNESCO(temp_ref,DDENS_DTEMP0,DDENS_DSAL0)

Victor Stepanenko
committed
DENS_REF = WATER_DENS_T_UNESCO(temp_ref)
DDENS_DTEMP0 = DDENS_DTEMP(temp_ref,sal_ref)

Victor Stepanenko
committed
DENS_REF = WATER_DENS_TS(temp_ref,sal_ref)

Victor Stepanenko
committed

Victor Stepanenko
committed
FUNCTION SET_DDENS_DTEMP(eos,lindens,Temp,Sal)

Victor Stepanenko
committed
integer(kind=iintegers), intent(in) :: eos !> Equation of state identifier
integer(kind=iintegers), intent(in) :: lindens !> Switch for linear density function

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: Temp, Sal

Victor Stepanenko
committed
real(kind=ireals) :: SET_DDENS_DTEMP

Victor Stepanenko
committed
real(kind=ireals) :: x

Victor Stepanenko
committed
if (lindens == 0) then
select case (eos)
case(eos_Hostetler)
call DDENS_HOSTETLER(Temp,SET_DDENS_DTEMP,x)
case(eos_TEOS)
print*, 'Not operational eos: STOP '
STOP
case(eos_Kivu)
call DDENS_KIVU(Temp,Sal,SET_DDENS_DTEMP,x)
case(eos_UNESCO)
call DDENS_UNESCO(Temp,SET_DDENS_DTEMP,x)
case(eos_McCatcheon)
SET_DDENS_DTEMP = DDENS_DTEMP(Temp,Sal)
end select
else
SET_DDENS_DTEMP = DDENS_DTEMP0
endif

Victor Stepanenko
committed
FUNCTION SET_DDENS_DSAL(eos,lindens,Temp,Sal)

Victor Stepanenko
committed
integer(kind=iintegers), intent(in) :: eos !> Equation of state identifier
integer(kind=iintegers), intent(in) :: lindens !> Switch for linear density function

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: Temp, Sal

Victor Stepanenko
committed
real(kind=ireals) :: SET_DDENS_DSAL

Victor Stepanenko
committed
real(kind=ireals) :: x

Victor Stepanenko
committed
if (lindens == 0) then
select case (eos)
case(eos_Hostetler)
call DDENS_HOSTETLER(Temp,x,SET_DDENS_DSAL)
case(eos_TEOS)
print*, 'Not operational eos: STOP '
STOP
case(eos_Kivu)
call DDENS_KIVU(Temp,Sal,x,SET_DDENS_DSAL)
case(eos_UNESCO)
call DDENS_UNESCO(Temp,x,SET_DDENS_DSAL)
case(eos_McCatcheon)
SET_DDENS_DSAL = DDENS_DSAL(Temp,Sal)
end select
else
SET_DDENS_DSAL = DDENS_DSAL0
endif

Victor Stepanenko
committed
!>Subroutine DENSITY_W calculates density profile
SUBROUTINE DENSITY_W(M,eos,lindens,Tw,Sal,preswat,row)

Victor Stepanenko
committed
integer(kind=iintegers), intent(in) :: eos, M, lindens

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: Tw(1:M+1), Sal(1:M+1), preswat(1:M+1)

Victor Stepanenko
committed
real(kind=ireals), intent(out) :: row(1:M+1)

Victor Stepanenko
committed
integer(kind=iintegers) :: i
!Calculation of water density profile at the current timestep

Victor Stepanenko
committed
if (lindens == 0) then
select case (eos)
case(eos_Hostetler)
do i = 1, M+1
row(i) = DENS_HOSTETLER(Tw(i))
enddo
case(eos_TEOS)
do i = 1, M+1
row(i) = GSW_RHO(Sal(i),Tw(i),preswat(i))
enddo
case(eos_Kivu)
do i = 1, M+1
row(i) = WATER_DENS_TS_KIVU(Tw(i),Sal(i)) !0.d0
enddo
case(eos_UNESCO)
do i = 1, M+1
row(i) = WATER_DENS_T_UNESCO(Tw(i)) !0.d0
enddo
case(eos_McCatcheon)
do i = 1, M+1
row(i) = WATER_DENS_TS(Tw(i),Sal(i)) ! McCutcheon et al., 1993
enddo
end select
else

Victor Stepanenko
committed
row(i) = DENS_REF + DDENS_DTEMP0*(Tw(i) - temp_ref) + &
& DDENS_DSAL0 *(Sal(i) - sal_ref)

Victor Stepanenko
committed
endif
real(kind=ireals) FUNCTION WATER_DENS_T_UNESCO(Temp)
!The function WATER_DENS_T_UNESCO
!returns the density of water, kg/m**3,
!as a function of temperature
!according to UNESCO formula
implicit none

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
WATER_DENS_T_UNESCO = &
& row0*(1+a0+a1*Temp-a2*Temp**2+a3*Temp**3)
END FUNCTION WATER_DENS_T_UNESCO

Victor Stepanenko
committed
SUBROUTINE DDENS_UNESCO(Temp,ddens_dtemp,ddens_dsal)
implicit none

Victor Stepanenko
committed

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(out):: ddens_dtemp, ddens_dsal
ddens_dtemp = &
& row0*(a1 - 2.*a2*Temp + 3.*a3*Temp**2)
ddens_dsal = 0.

Victor Stepanenko
committed
FUNCTION WATER_DENS_TS_KIVU(Temp,Sal)
! The function WATER_DENS_TS_KIVU
! resturns the density of water, kg/m**3,
! as a function of temperature and salinity
! adjusted for Lake Kivu (Rwanda & DRC) according to
! (Schmid et al. 2003)
! Input variables:
! Temp --- water temperature, deg C
! Sal --- water salinity, kg/kg
implicit none

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(in):: Sal

Victor Stepanenko
committed
real(kind=ireals) :: WATER_DENS_TS_KIVU

Victor Stepanenko
committed
real(kind=ireals) :: Sal_g_kg,A,B,C
! Converting salinity units from kg/kg to g/kg
Sal_g_kg = Sal*1.d+3
WATER_DENS_TS_KIVU = (at1 + at2*Temp + at3*Temp**2 + at4*Temp**3)* &
& (1. + beta*Sal_g_kg + beta_co2*co2_ref + beta_ch4*ch4_ref)
SUBROUTINE DDENS_KIVU(Temp,Sal,ddens_dtemp,ddens_dsal)
implicit none

Victor Stepanenko
committed

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: Temp, Sal
real(kind=ireals), intent(out) :: ddens_dtemp, ddens_dsal

Victor Stepanenko
committed

Victor Stepanenko
committed
real(kind=ireals) :: Sal_g_kg

Victor Stepanenko
committed
! Converting salinity units from kg/kg to g/kg
Sal_g_kg = Sal*1.d+3

Victor Stepanenko
committed
ddens_dtemp = (at2 + 2.*at3*Temp + 3.*at4*Temp**2)* &
& (1. + beta*Sal_g_kg + beta_co2*co2_ref + beta_ch4*ch4_ref)
ddens_dsal = (at1 + at2*Temp + at3*Temp**2 + at4*Temp**3)* &
& beta*1.d+3

Victor Stepanenko
committed

Victor Stepanenko
committed
real(kind=ireals) FUNCTION WATER_DENS_TS(Temp,Sal)
! The function WATER_DENS_TS
! resturns the density of water, kg/m**3,
! as a function of temperature and salinity
! according to
! (McCutcheon et al.,Water Quality and Maidment.
! Handbook on Hydrology, 1993)
! Temp --- water temperature, deg C
! Sal --- water salinity, kg/kg
implicit none

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(in):: Sal

Victor Stepanenko
committed
real(kind=ireals) Sal_g_kg,A,B,C
! Converting salinity units from kg/kg to g/kg
Sal_g_kg = Sal*1.d+3
! The first term, dependent on temperature
WATER_DENS_TS = row0 * &
& ( 1.-(Temp+A11)*(Temp-A14)**2/(A12*(Temp+A13) ) )
A = alpha1 - alpha2*Temp + alpha3*Temp**2 &
& -alpha4*Temp**3 + alpha5*Temp**4
! Adding the second term, dependent on temperature and salinity
WATER_DENS_TS = WATER_DENS_TS + &
& A*Sal_g_kg + B*Sal_g_kg**1.5 + C*Sal_g_kg**2.
real(kind=ireals) FUNCTION DDENS_DTEMP(Temp,Sal)
! The function DDENS_DTEMP
! returns the derivative of water density
! on temperature, kg/(m**3*C)
! Input variables:
! Temp --- water temperature, deg C
! Sal --- water salinity, kg/kg

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(in):: Sal
real(kind=ireals) Sal_g_kg,A,B,C
& -(row0/A12)*(Temp-A14)/( (Temp+A13)**2 ) * &
& ( (Temp+A13)*( (Temp-A14)+2.*(Temp+A11) ) - &
& (Temp+A11)*(Temp-A14) )
! Converting salinity units from kg/kg to g/kg
Sal_g_kg = Sal*1.d+3
A = - alpha2 + 2.*alpha3*Temp - 3.*alpha4*Temp**2 + 4.*alpha5*Temp**3
B = + beta2 - 2.*beta3*Temp
DDENS_DTEMP = DDENS_DTEMP + A*Sal_g_kg + B*Sal_g_kg**1.5
real(kind=ireals) FUNCTION DDENS_DSAL(Temp,Sal)
! The function DDENS_DSAL
! returns the derivative of water density
! on salinity, , kg/(m**3*kg/kg)
! Input variables:
! Temp --- water temperature, deg C
! Sal --- water salinity, kg/kg

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(in):: Sal

Victor Stepanenko
committed
real(kind=ireals) Sal_g_kg,A,B,C
A = alpha1 - alpha2*Temp + alpha3*Temp**2 &
& - alpha4*Temp**3 + alpha5*Temp**4
B = - beta1 + beta2*Temp - beta3*Temp**2
C = gamma1
! Converting salinity units from kg/kg to g/kg
Sal_g_kg = Sal*1.d+3
END FUNCTION DDENS_DSAL
FUNCTION DENS_HOSTETLER(temp)
! Function DENS_HOSTETLER calculates the density of water
! dependent on temperature only, following Hostetler and Bartlein, 1990

Victor Stepanenko
committed
real(kind=ireals) :: DENS_HOSTETLER

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: temp

Victor Stepanenko
committed
real(kind=ireals), parameter :: row0 = 1.d+3
real(kind=ireals), parameter :: temp0 = 3.85
real(kind=ireals), parameter :: const = 1.9549d-5
real(kind=ireals), parameter :: const_power = 1.68
DENS_HOSTETLER = row0*(1.-const*abs(temp-temp0)**const_power)
SUBROUTINE DDENS_HOSTETLER(Temp,ddens_dtemp,ddens_dsal)
implicit none
! Parameters

Victor Stepanenko
committed
real(kind=ireals), parameter :: row0 = 1.d+3
real(kind=ireals), parameter :: temp0 = 3.85
real(kind=ireals), parameter :: const = 1.9549d-5
real(kind=ireals), parameter :: const_power = 1.68

Victor Stepanenko
committed
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(out):: ddens_dtemp, ddens_dsal
ddens_dtemp = row0*( - const_power*const*abs(temp-temp0)**(const_power-1)* &
& sign(1._ireals,temp-temp0))
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
ddens_dsal = 0.
END SUBROUTINE DDENS_HOSTETLER
!--------------------------------------------------------------------------
!--------------------------------------------------------------------------
! density and enthalpy, based on the 48-term expression for density from TEOS-2010
!--------------------------------------------------------------------------
!==========================================================================
FUNCTION GSW_RHO(sa_,ct_,p_)
!==========================================================================
! Calculates in-situ density from Absolute Salinity and Conservative
! Temperature, using the computationally-efficient 48-term expression for
! density in terms of SA, CT and p (McDougall et al., 2011).
!
! sa : Absolute Salinity [kg/kg]
! ct : Conservative Temperature [deg C]
! p : sea pressure [Pa]
!
! gsw_rho : in-situ density (48 term equation)
implicit none
integer, parameter :: r14 = selected_real_kind(14,30)
real (r14), parameter :: v01 = 9.998420897506056d2, v02 = 2.839940833161907d0
real (r14), parameter :: v03 = -3.147759265588511d-2, v04 = 1.181805545074306d-3
real (r14), parameter :: v05 = -6.698001071123802d0, v06 = -2.986498947203215d-2
real (r14), parameter :: v07 = 2.327859407479162d-4, v08 = -3.988822378968490d-2
real (r14), parameter :: v09 = 5.095422573880500d-4, v10 = -1.426984671633621d-5
real (r14), parameter :: v11 = 1.645039373682922d-7, v12 = -2.233269627352527d-2
real (r14), parameter :: v13 = -3.436090079851880d-4, v14 = 3.726050720345733d-6
real (r14), parameter :: v15 = -1.806789763745328d-4, v16 = 6.876837219536232d-7
real (r14), parameter :: v17 = -3.087032500374211d-7, v18 = -1.988366587925593d-8
real (r14), parameter :: v19 = -1.061519070296458d-11, v20 = 1.550932729220080d-10
real (r14), parameter :: v21 = 1.0d0, v22 = 2.775927747785646d-3, v23 = -2.349607444135925d-5
real (r14), parameter :: v24 = 1.119513357486743d-6, v25 = 6.743689325042773d-10
real (r14), parameter :: v26 = -7.521448093615448d-3, v27 = -2.764306979894411d-5
real (r14), parameter :: v28 = 1.262937315098546d-7, v29 = 9.527875081696435d-10
real (r14), parameter :: v30 = -1.811147201949891d-11, v31 = -3.303308871386421d-5
real (r14), parameter :: v32 = 3.801564588876298d-7, v33 = -7.672876869259043d-9
real (r14), parameter :: v34 = -4.634182341116144d-11, v35 = 2.681097235569143d-12
real (r14), parameter :: v36 = 5.419326551148740d-6, v37 = -2.742185394906099d-5
real (r14), parameter :: v38 = -3.212746477974189d-7, v39 = 3.191413910561627d-9
real (r14), parameter :: v40 = -1.931012931541776d-12, v41 = -1.105097577149576d-7
real (r14), parameter :: v42 = 6.211426728363857d-10, v43 = -1.119011592875110d-10
real (r14), parameter :: v44 = -1.941660213148725d-11, v45 = -1.864826425365600d-14
real (r14), parameter :: v46 = 1.119522344879478d-14, v47 = -1.200507748551599d-15
real (r14), parameter :: v48 = 6.057902487546866d-17
real (kind=ireals), intent(in) :: sa_, ct_, p_
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
real (r14) :: sa, ct, p, sqrtsa, v_hat_denominator, &
& v_hat_numerator, gsw_rho
sa = sa_*1.d+3 ! Converting kg/kg to g/kg
p = p_*1.d-4 ! Converting pressure from Pa to dbar
ct = ct_
sqrtsa = sqrt(sa)
v_hat_denominator = v01 + ct*(v02 + ct*(v03 + v04*ct)) &
+ sa*(v05 + ct*(v06 + v07*ct) &
+ sqrtsa*(v08 + ct*(v09 + ct*(v10 + v11*ct)))) &
+ p*(v12 + ct*(v13 + v14*ct) + sa*(v15 + v16*ct) &
+ p*(v17 + ct*(v18 + v19*ct) + v20*sa))
v_hat_numerator = v21 + ct*(v22 + ct*(v23 + ct*(v24 + v25*ct))) &
+ sa*(v26 + ct*(v27 + ct*(v28 + ct*(v29 + v30*ct))) + v36*sa &
+ sqrtsa*(v31 + ct*(v32 + ct*(v33 + ct*(v34 + v35*ct))))) &
+ p*(v37 + ct*(v38 + ct*(v39 + v40*ct)) &
+ sa*(v41 + v42*ct) &
+ p*(v43 + ct*(v44 + v45*ct + v46*sa) &
+ p*(v47 + v48*ct)))
GSW_RHO = v_hat_denominator/v_hat_numerator
return
END FUNCTION GSW_RHO

Victor Stepanenko
committed
use LAKE_DATATYPES, only : ireals, iintegers
real(kind=ireals) FUNCTION TURB_DENS_FLUX(tempflux,salflux,Temp,Sal)
! Function TURB_DENS_FLUX _____
! returns the turbulent density flux w'ro' in water
! at given temperature,
! salinity, ____
! temperature flux w'T'
! ____
! and salinity flux w's'
! Input variables:
! tempflux --- kinematic heat flux, m*C/s
! salflux --- salinity flux, m*(kg/kg)/s
! Temp --- temperature , C
! Sal --- salinity , kg/kg
use water_density, only : &
& ddens_dtemp, &
& ddens_dsal
implicit none

Victor Stepanenko
committed
real(kind=ireals), intent(in):: tempflux
real(kind=ireals), intent(in):: salflux
real(kind=ireals), intent(in):: Temp
real(kind=ireals), intent(in):: Sal
TURB_DENS_FLUX = &
& ddens_dtemp(Temp,Sal)*tempflux + &
& ddens_dsal (Temp,Sal)*salflux
END FUNCTION TURB_DENS_FLUX
! real(kind=ireals) FUNCTION dirdif()

Victor Stepanenko
committed
! real(kind=ireals) cloud
! real(kind=ireals) dirdif0,b,S0,sinh0
! SAVE
! b=1./3.
! S0=1367.
! cloud = 0.
! dirdif0 = (shortwave-b*S0*sinh0())/(b*(S0*sinh0()-shortwave))
! END FUNCTION dirdif
FUNCTION SINH0(year,month,day,hour,phi)
! SINH0 is sine of solar angle ( = cosine of zenith angle)

Victor Stepanenko
committed
use TIME_LAKE_MOD, only : JULIAN_DAY

Victor Stepanenko
committed
real(kind=ireals) :: sinh0

Victor Stepanenko
committed
integer(kind=iintegers), intent(in) :: year
integer(kind=iintegers), intent(in) :: month
integer(kind=iintegers), intent(in) :: day

Victor Stepanenko
committed
real(kind=ireals) , intent(in) :: hour
real(kind=ireals) , intent(in) :: phi

Victor Stepanenko
committed
real(kind=ireals) delta
real(kind=ireals) theta
real(kind=ireals) pi
real(kind=ireals) phi_rad

Victor Stepanenko
committed
integer(kind=iintegers) nday
pi=4.*datan(1.d0)
nday = JULIAN_DAY(year,month,day)
delta = 23.5d0*pi/180.d0*cos(2*pi*(float(nday)-173.d0)/365.d0)
sinh0 = sin(phi_rad)*sin(delta) + &
& cos(phi_rad)*cos(delta)*cos(theta)

Victor Stepanenko
committed
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
!> 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
implicit none

Victor Stepanenko
committed
real(kind=ireals) t,p,aw,bw,ai,bi,a,b,es
integer(kind=iintegers) phase
aw = 7.6326
bw = 241.9
ai = 9.5
bi = 265.5
a = phase*aw+(1-phase)*ai
b = phase*bw+(1-phase)*bi
es=610.7*10.**(a*t/(b+t))
QS=0.622*es/p
END FUNCTION

Victor Stepanenko
committed
!> Derivative of saturated specific humidity on temperature
FUNCTION DQSDT(phase,t,p)
! QS - specific humidity, kg/kg, for saturated water vapour
implicit none
real(kind=ireals) :: t,p,aw,bw,ai,bi,a,b,qs,es
real(kind=ireals) :: DQSDT
integer(kind=iintegers) :: phase
!phase = 1 - liquid, 0 - ice
aw = 7.6326
bw = 241.9
ai = 9.5
bi = 265.5
a = phase*aw+(1-phase)*ai
b = phase*bw+(1-phase)*bi
es = 610.7*10.**(a*t/(b+t))
qs = 0.622*es/p
DQSDT = qs*a*b*log(10.)/(b+t)**2
END FUNCTION DQSDT
real(kind=ireals) FUNCTION ES(phase,t,p)

Victor Stepanenko
committed
real(kind=ireals) t,p,aw,bw,ai,bi,a,b
integer(kind=iintegers) phase
aw = 7.6326
bw = 241.9
ai = 9.5
bi = 265.5
a = phase*aw+(1-phase)*ai
b = phase*bw+(1-phase)*bi
ES = 610.7*10.**(a*t/(b+t))
END FUNCTION
real(kind=ireals) FUNCTION MELTPNT(C,p,npar)

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: C ! salinity, kg/kg
real(kind=ireals), intent(in) :: p ! pressure, Pa
integer(kind=iintegers), intent(in) :: npar

Victor Stepanenko
committed
real(kind=ireals), parameter :: dtdc = 66.7, pnorm = 101325
real(kind=ireals), parameter :: Pa2dbar = 1.d-4
select case(npar)
case(0)
MELTPNT = 0.
case(1)
MELTPNT = 0. - C*dtdc
case(2)
MELTPNT = FP_THETA(C*1.E+3_ireals, (p-pnorm)*Pa2dbar,'air-sat') ! convertion to 0/00
end select
END FUNCTION MELTPNT
FUNCTION FP_THETA(s,p,sat)
! potential temperature freezing point of seawater, as in
! Jackett, McDougall, Feistel, Wright and Griffies (2004), submitted JAOT
!
! s : salinity (psu)
! p : gauge pressure (dbar)
! (absolute pressure - 10.1325 dbar)
! sat : string variable
! 'air-sat' - air saturated water
! 'air-free' - air free water
!
! fp_theta : potential freezing temperature (deg C, ITS-90)
!
! check value : fp_theta(35,200,'air-sat') = -2.076426227617581 deg C
! fp_theta(35,200,'air-free') = -2.074408175943127 deg C
!
! DRJ on 2/6/04
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
character*(*) sat
sqrts = sqrt(s)
tf_num = 2.5180516744541290d-03 + &
s*(-5.8545863698926184d-02 + &
sqrts*( 2.2979985780124325d-03 - &
sqrts * 3.0086338218235500d-04)) + &
p*(-7.0023530029351803d-04 + &
p*( 8.4149607219833806d-09 + &
s * 1.1845857563107403d-11));
tf_den = 1.0000000000000000d+00 + &
p*(-3.8493266309172074d-05 + &
p * 9.1686537446749641d-10) + &
s*s*sqrts* 1.3632481944285909d-06
fp_theta = tf_num/tf_den;
if(sat.eq.'air-sat') then
fp_theta = fp_theta -2.5180516744541290d-03 + &
s * 1.428571428571429d-05
elseif(sat.eq.'air-free') then
continue
else
print *, '*** Error in fp_theta.f90: invalid third argument ***'
print *
stop
endif
return
END FUNCTION FP_THETA
FUNCTION WATER_FREEZE_MELT(temperature, grid_size, &
& melting_point, switch)
! The function WATER_FREEZE_MELT checks, if the
! heat storage of a grid cell is larger or equal
! the latent heat, necessary for phase transition
use PHYS_CONSTANTS, only : &
& cw_m_row0, &
& ci_m_roi, &
& row0_m_Lwi, &
& roi_m_Lwi
use NUMERIC_PARAMS, only : &
& min_ice_thick, &
& min_water_thick, &
& T_phase_threshold
implicit none
! Input variables

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: temperature ! Temperature, Celsius
real(kind=ireals), intent(in) :: grid_size
real(kind=ireals), intent(in) :: melting_point

Victor Stepanenko
committed
integer(kind=iintegers), intent(in) :: switch ! +1 for water -> ice transition
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
! -1 for ice -> water transition
logical :: WATER_FREEZE_MELT
if (switch == +1) then
if ( ( melting_point - T_phase_threshold - temperature) * &
& cw_m_row0*grid_size > min_ice_thick*roi_m_Lwi .or. &
& ( melting_point - T_phase_threshold - temperature) > 0.2d0) &
& then
WATER_FREEZE_MELT = .true.
else
WATER_FREEZE_MELT = .false.
endif
elseif (switch == -1) then
if ( (temperature - T_phase_threshold - melting_point) * &
& ci_m_roi*grid_size > min_water_thick*row0_m_Lwi .or. &
& (temperature - T_phase_threshold - melting_point) > 0.2d0) &
& then
WATER_FREEZE_MELT = .true.
else
WATER_FREEZE_MELT = .false.
endif
else
write(*,*) 'Wrong switch in WATER_FREEZE_MELT: STOP'
STOP
endif
END FUNCTION WATER_FREEZE_MELT
!> Subroutine TURB_SCALES calculates turbulent scales, both bulk and local

Victor Stepanenko
committed
SUBROUTINE TURB_SCALES(gsp, ls, bathymwater, wst, RadWater, &
& k_turb_T_flux, T_massflux, row, eflux0_kinem, &
& turb_density_flux, Buoyancy0, tau, kor, i_maxN, H_mixed_layer, maxN, w_conv_scale, &
& T_conv_scale, Wedderburn, LakeNumber, Rossby_rad, ThermThick, ReTherm, RiTherm, trb)

Victor Stepanenko
committed
use DRIVING_PARAMS, only : M, eos, lindens

Victor Stepanenko
committed
use PHYS_CONSTANTS, only : g, row0, cw_m_row0, g_d_Kelvin0, niu_wat
use ARRAYS_BATHYM, only : bathym, layers_type
use ARRAYS_GRID, only : gridspacing_type

Victor Stepanenko
committed
use WATER_DENSITY, only : DDENS_DTEMP0, DDENS_DSAL0
use NUMERIC_PARAMS, only : small_value
use ARRAYS_WATERSTATE, only : waterstate_type

Victor Stepanenko
committed
use RADIATION, only : rad_type
use ARRAYS_TURB, only : turb_type
use WATER_DENSITY, only : SET_DDENS_DTEMP, SET_DDENS_DSAL
type(gridspacing_type), intent(in) :: gsp
type(layers_type), intent(in) :: ls
type(bathym), intent(in) :: bathymwater(1:M+1)
type(waterstate_type), intent(in) :: wst

Victor Stepanenko
committed
type(rad_type), intent(in) :: RadWater

Victor Stepanenko
committed
real(kind=ireals), intent(in) :: k_turb_T_flux(1:M)
real(kind=ireals), intent(in) :: T_massflux(1:M)
real(kind=ireals), intent(in) :: row(1:M+1)
real(kind=ireals), intent(in) :: eflux0_kinem
real(kind=ireals), intent(in) :: turb_density_flux
real(kind=ireals), intent(in) :: Buoyancy0
real(kind=ireals), intent(in) :: tau
real(kind=ireals), intent(in) :: kor
integer(kind=iintegers), intent(out) :: i_maxN

Victor Stepanenko
committed
real(kind=ireals), intent(out) :: H_mixed_layer
real(kind=ireals), intent(out) :: maxN

Victor Stepanenko
committed
real(kind=ireals), intent(out) :: w_conv_scale
real(kind=ireals), intent(out) :: T_conv_scale
real(kind=ireals), intent(out) :: Wedderburn !Wedderburn number
real(kind=ireals), intent(out) :: LakeNumber
real(kind=ireals), intent(out) :: Rossby_rad !Internal Rossby radius

Victor Stepanenko
committed
real(kind=ireals), intent(out) :: ThermThick !Thermocline thickness
real(kind=ireals), intent(out) :: ReTherm !Reynolds number in thermocline
real(kind=ireals), intent(out) :: RiTherm !Richardson number in thermocline
type(turb_type), intent(inout) :: trb

Victor Stepanenko
committed
real(kind=ireals), external:: DZETA

Victor Stepanenko
committed
real(kind=ireals), parameter :: Ttherm0 = 14., Ttherm1 = 8. !temperatures of top and bottom of thermocline, C
! Auxilliary variables
integer(kind=iintegers) :: maxlocat, i

Victor Stepanenko
committed
real(kind=ireals), allocatable :: bvf(:)
real(kind=ireals) :: zv, zm, x, y, x1, x2

Victor Stepanenko
committed
real(kind=ireals) :: u1, u2, v1, v2, Sal1, Sal2
! Mixed-layer depth - depth of maximal heat flux
! maxlocat = 1
! do i = 2, M
! if (k_turb_T_flux(i) + T_massflux(i) > &
! & k_turb_T_flux(maxlocat) + T_massflux(maxlocat)) maxlocat = i
! enddo
! H_mixed_layer = dzeta_05int(maxlocat)*h1
! Mixed-layer depth - depth of maximal Bruent-Vaisala frequency
! allocate (bvf(2:M-1))
! do i = 2, M-1
! bvf(i) = g/row0*(row(i+1) - row(i-1))/(h1*(ddz(i-1) + ddz(i)))
! enddo
! maxlocat = maxloc(bvf,1) + lbound(bvf,1) - 1
! deallocate (bvf)
! H_mixed_layer = dzeta_int(maxlocat)*h1
call MIXED_LAYER_CALC(row,gsp%ddz,gsp%ddz2,gsp%dzeta_int,gsp%dzeta_05int, &
& ls%h1,M,i_maxN,H_mixed_layer,maxN)

Victor Stepanenko
committed
w_conv_scale = max( (Buoyancy0 - &
& (RadWater%integr(1) - RadWater%integr(i_maxN))*g_d_Kelvin0/cw_m_row0) * &
& H_mixed_layer, 0._ireals)**(1._ireals/3._ireals)
Wedderburn = max(min(g*(row(M+1) - row(1))*H_mixed_layer*H_mixed_layer/&
& ((tau+small_value)*max(bathymwater(1)%Lx,bathymwater(1)%Ly)),&
! Depth of volume center (zv) and of center of mass (zm)
zm = 0.; x1 = 0.
zv = 0.; x2 = 0.
do i = 1, M+1