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
contains
SUBROUTINE MON_OUT(ix,iy,nx,ny,year,month,day,hour,time,z_out,nout)
use ARRAYS_WATERSTATE, only : Tw1, Sal1
use ARRAYS_TURB, only : E1, E2, eps1, S, Gen, Gen_seiches, row, TKE_turb_trans, TF, KT
use ARRAYS, only : u1, v1
use ARRAYS_GRID, only : dzeta_int, z_full, z_half
use ARRAYS_BATHYM, only : h1, l1, voldef, vol
use DRIVING_PARAMS, only : M
use PHYS_CONSTANTS, only : cw_m_row0
use INOUT_PARAMETERS, only : &
& lake_mon_out_unit_min, &
& lake_mon_out_unit_max
implicit none
! Input variables
! Reals
real(kind=ireals), intent(in) :: hour
real(kind=ireals), intent(in) :: time
real(kind=ireals), intent(in) :: z_out(1:nout) ! The output grid
! Integers
integer(kind=iintegers), intent(in) :: year, month, day
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: nout
! Local variables
! Reals
real(kind=ireals), allocatable :: accum_var (:,:,:,:)
real(kind=ireals), allocatable :: var (:,:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar (:,:,:)
real(kind=ireals), allocatable :: var_scalar (:,:,:)
real(kind=ireals), allocatable :: tsteps (:,:)
real(kind=ireals), allocatable :: Profile_out(:,:)
real(kind=ireals), allocatable :: work(:)
real(kind=ireals) :: voldef0, voldef0y
real(kind=ireals) :: work1, work2
real(kind=ireals), allocatable :: valmax(:)
real(kind=ireals), external:: DZETA
! Integers
integer(kind=iintegers), parameter :: n_var = 11, n_var_scalar = 5, n_var_max = 3
integer(kind=iintegers), parameter :: n_var1 = 4 ! Number of vector variables defined at cell interfaces
integer(kind=iintegers), allocatable :: month_old(:,:)
integer(kind=iintegers), allocatable :: out_unita(:,:)
integer(kind=iintegers) :: out_unit = lake_mon_out_unit_min
integer(kind=iintegers) :: i, j ! Loop indices
! Characters
character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: timestring*6
character :: coords_point*6
character :: formatline*50
! Logicals
logical :: firstcall, flag
logical, allocatable :: firstcallixiy(:,:)
data firstcall /.true./
SAVE
if (firstcall) then
allocate (firstcallixiy(1:nx,1:ny))
allocate (out_unita(1:nx,1:ny))
allocate (month_old(1:nx, 1:ny))
allocate (tsteps (1:nx, 1:ny))
allocate (var (1:n_var, 1:nx, 1:ny, 1:M+1) )
allocate (accum_var(1:n_var, 1:nx, 1:ny, 1:M+1) )
allocate (var_scalar (1:n_var_scalar, 1:nx, 1:ny) )
allocate (accum_var_scalar(1:n_var_scalar, 1:nx, 1:ny) )
allocate (valmax(1:n_var_max))
allocate (work(1:M+1))
firstcallixiy(:,:) = .true.
out_unita(:,:) = lake_mon_out_unit_min
month_old(:,:) = month
tsteps (:,:) = 0.d0
accum_var(:,:,:,:) = 0.d0
accum_var_scalar(:,:,:) = 0.d0
voldef0 = 0.
voldef0y = 0.
valmax(:) = 0.
endif
if (firstcallixiy(ix,iy)) then
call CHECK_UNIT(lake_mon_out_unit_min,lake_mon_out_unit_max,out_unita(ix,iy))
write (coords_point,'(2i3)') ix, iy
open(out_unita(ix,iy),file=outpath(1:len_trim(outpath))//'monthly/'// &
& 'depvol'//coords_point//'.dat', status='unknown')
write (out_unita(ix,iy),*) '1 - year'
write (out_unita(ix,iy),*) '2 - month'
write (out_unita(ix,iy),*) '3 - depth, m'
write (out_unita(ix,iy),*) '4 - ice thickness, m'
write (out_unita(ix,iy),*) '5 - volume, m**3'
write (out_unita(ix,iy),*) '6 - volume deficit, m**3'
write (out_unita(ix,iy),*) '7 - accumulated volume deficit, m**3'
write (out_unita(ix,iy),*) '8 - maximal monthly depth, m'
write (out_unita(ix,iy),*) '9 - maximal monthly ice thickness, m'
write (out_unita(ix,iy),*) '10 - maximal monthly volume, m**3'
endif
! Variables located at cell interfaces
var(1,ix,iy,1:M+1) = Tw1 (1:M+1)
var(2,ix,iy,1:M+1) = Sal1(1:M+1)
var(3,ix,iy,1:M+1) = row (1:M+1)
var(4,ix,iy,1:M+1) = sqrt(u1 (1:M+1)*u1 (1:M+1) + v1 (1:M+1)*v1 (1:M+1) )
! Variables located at cell centers
var(5,ix,iy,1:M) = E1 (1:M)
var(6,ix,iy,1:M) = TF(1:M)*E2(1:M) !eps1(1:M)
var(7,ix,iy,1:M) = S (1:M)
var(8,ix,iy,1:M) = Gen (1:M)
var(9,ix,iy,1:M) = Gen_seiches (1:M)
var(10,ix,iy,1:M) = TKE_turb_trans (1:M)
var(11,ix,iy,1:M) = KT (1:M)/cw_m_row0 !Kinematic heat conductance, m**2/s
var_scalar(1,ix,iy) = h1
var_scalar(2,ix,iy) = l1
var_scalar(3,ix,iy) = vol
var_scalar(4,ix,iy) = voldef - voldef0
var_scalar(5,ix,iy) = voldef - voldef0
voldef0 = voldef
valmax(1) = max(h1,valmax(1)) ! Calculating monthly maximum
valmax(2) = max(l1,valmax(2))
valmax(3) = max(vol,valmax(3))
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + &
& var_scalar(:,ix,iy)
tsteps(ix,iy) = tsteps(ix,iy) + 1._ireals
if (month_old(ix,iy) /= month) then
!print*, time/(365./12.*24.*3600.)
!read*
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:)/tsteps(ix,iy)
accum_var_scalar(1:3,ix,iy) = accum_var_scalar(1:3,ix,iy)/tsteps(ix,iy) ! excepting volume deficit
call DATEMINUS(1,year,month,day,hour,year1,month1,day1,hour1)
call TIMESTR(6,year1,month1,day1,hour1,timestring)
call CHECK_UNIT(lake_mon_out_unit_min,lake_mon_out_unit_max,out_unit)
write (coords_point,'(2i3)') ix, iy
open(out_unit,file=outpath(1:len_trim(outpath))//'monthly/'// &
& 'Profiles'//coords_point//timestring//'.dat', status='unknown')
write (out_unit,*) '1 - depth, m'
write (out_unit,*) '2 - temperature, C'
write (out_unit,*) '3 - salinity, kg/kg'
write (out_unit,*) '4 - density, kg/m**3'
write (out_unit,*) '5 - absolute value of velocity, m/s'
write (out_unit,*) '6 - turbulent kinetic energy, m**2/s**2'
write (out_unit,*) '7 - disspation rate, m**2/s**3'
write (out_unit,*) '8 - TKE production by buoyancy, m**2/s**3'
write (out_unit,*) '9 - TKE production by shear , m**2/s**3'
write (out_unit,*) '10 - TKE production by seiches (Goudsmit paramaterization) , m**2/s**3'
write (out_unit,*) '11 - TKE turbulent transport , m**2/s**3'
write (out_unit,*) '12 - kinematic heat conductance, m**2/s'
formatline = '(f12.6,f9.3,e12.4,f9.3,8e12.4)'
if (nout > 0) then
allocate (Profile_out(1:nout,1:n_var))
! Interpolating to given output levels
do j = 1, n_var1
work(1:M+1) = accum_var(j,ix,iy,1:M+1)
call LININTERPOL (z_full,work,M+1,z_out,Profile_out(1,j),nout,flag)
enddo
do j = n_var1+1, n_var
work(1:M) = accum_var(j,ix,iy,1:M)
call LININTERPOL (z_half,work,M,z_out,Profile_out(1,j),nout,flag)
enddo
j = nout
if (z_out(nout) > z_full(M+1)) then
do while (z_out(j) > z_full(M+1))
j = j - 1
enddo
endif
! Writing to file
do i = 1, j
write (unit = out_unit, fmt = formatline) &
& -z_out(i), Profile_out(i,1:n_var)
enddo
else
allocate (Profile_out(1:M+1,1:n_var))
forall (i = 1:M+1, j = 1:n_var1) Profile_out(i,j) = accum_var(j,ix,iy,i)
! Interpolating from cell centers to cell interfaces
do j = n_var1+1, n_var
work(1:M) = accum_var(j,ix,iy,1:M)
call LININTERPOL (z_half,work,M,z_full,Profile_out(1,j),M+1,flag)
enddo
! Writing to file
do i = 1, M+1
write (unit = out_unit, fmt = formatline) &
& -dzeta_int(i)*h1, Profile_out(i,1:n_var)
enddo
endif
deallocate(Profile_out)
close(out_unit)
read(year1,*) work1
read(month1,*) work2
write (out_unita(ix,iy),fmt=*) & !, fmt = formatline) &
& int(work1), int(work2), accum_var_scalar(1:n_var_scalar,ix,iy), &
& valmax(:)
tsteps(ix,iy) = 0.d0
accum_var(:,ix,iy,:) = 0.d0
accum_var_scalar(1:4,ix,iy) = 0.d0
if (month_old(ix,iy) == 12) accum_var_scalar(5,ix,iy) = 0.d0
valmax(:) = 0.d0
month_old(ix,iy) = month
endif
if (firstcall) firstcall = .false.
if (firstcallixiy(ix,iy)) firstcallixiy(ix,iy) = .false.
END SUBROUTINE MON_OUT
SUBROUTINE DAY_OUT(ix,iy,nx,ny,year,month,day,hour)
use ARRAYS_GRID, only : dzeta_int
use ARRAYS_WATERSTATE, only : Tw1, Sal1
use ARRAYS_TURB, only : E1, eps1, KT, Ri
use ARRAYS_BATHYM, only : h1
use ARRAYS, only : u1, v1
use ATMOS , only : discharge
use DRIVING_PARAMS, only : M
use PHYS_CONSTANTS, only : cw_m_row0
use INOUT_PARAMETERS, only : &
& lake_day_out_unit_min, &
& lake_day_out_unit_max
implicit none
! Input variables
! Reals
real(kind=ireals), intent(in) :: hour
! Integers
integer(kind=iintegers), intent(in) :: year, month, day
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
! Local variables
! Reals
real(kind=ireals), allocatable :: tsteps (:,:)
real(kind=ireals), allocatable :: accum_var(:,:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar (:,:,:)
real(kind=ireals), allocatable :: var_scalar (:,:,:)
real(kind=ireals), allocatable :: var (:,:,:,:)
real(kind=ireals), external :: DZETA
real(kind=ireals) :: work1, work2, work3
! Integers
integer(kind=iintegers), allocatable :: day_old(:,:)
integer(kind=iintegers), allocatable :: out_unita(:,:)
integer(kind=iintegers) :: i ! Loop index
integer(kind=iintegers) :: out_unit = lake_day_out_unit_min
integer(kind=iintegers) :: nvars = 7
integer(kind=iintegers) :: nvars_scalar = 2
character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: timestring*8
character :: coords_point*6
logical :: firstcall
logical, allocatable :: firstcallixiy(:,:)
data firstcall /.true./
SAVE
if (firstcall) then
allocate (firstcallixiy(1:nx,1:ny))
allocate (out_unita(1:nx,1:ny))
allocate (day_old (1:nx, 1:ny))
allocate (tsteps (1:nx, 1:ny))
allocate (var (1:nvars, 1:nx, 1:ny, 1:M+1) )
allocate (var_scalar(1:nvars_scalar, 1:nx, 1:ny ) )
allocate (accum_var (1:nvars, 1:nx, 1:ny, 1:M+1) )
allocate (accum_var_scalar (1:nvars_scalar, 1:nx, 1:ny) )
firstcallixiy(:,:) = .true.
out_unita(:,:) = lake_day_out_unit_min
day_old (:,:) = day
tsteps (:,:) = 0.d0
accum_var(:,:,:,:) = 0.d0
accum_var_scalar(:,:,:) = 0.d0
var (:,:,:,:) = 0.d0
var_scalar(:,:,:) = 0.d0
endif
if (firstcallixiy(ix,iy)) then
call CHECK_UNIT(lake_day_out_unit_min,lake_day_out_unit_max,out_unita(ix,iy))
write (coords_point,'(2i3)') ix, iy
open(out_unita(ix,iy),file=outpath(1:len_trim(outpath))//'daily/'// &
& 'fluxes'//coords_point//'.dat', status='unknown')
write (out_unita(ix,iy),*) '1 - year'
write (out_unita(ix,iy),*) '2 - month'
write (out_unita(ix,iy),*) '3 - discharge in x direction, m**3/s'
write (out_unita(ix,iy),*) '4 - discharge in y direction, m**3/s'
endif
!Vector variables
var(1,ix,iy,1:M+1) = Tw1 (1:M+1)
var(2,ix,iy,1:M+1) = Sal1(1:M+1)
var(3,ix,iy,1:M) = E1 (1:M)
var(4,ix,iy,1:M) = eps1(1:M)
var(5,ix,iy,1:M) = KT(1:M)/cw_m_row0
var(6,ix,iy,1:M+1) = sqrt(u1(1:M+1)**2 + v1(1:M+1)**2)
var(7,ix,iy,1:M) = Ri(1:M)
!Scalar variables
var_scalar(1,ix,iy) = discharge(1)
var_scalar(2,ix,iy) = discharge(2)
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + var_scalar(:,ix,iy)
tsteps(ix,iy) = tsteps(ix,iy) + 1.d0
if (day_old(ix,iy)/=day) then
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:)/tsteps(ix,iy)
call DATEMINUS (2,year,month,day,hour,year1,month1,day1,hour1)
call TIMESTR (8,year1,month1,day1,hour1,timestring)
call CHECK_UNIT(lake_day_out_unit_min,lake_day_out_unit_max,out_unit)
write (coords_point, '(2i3)') ix, iy
open(out_unit, file=outpath(1:len_trim(outpath))//'daily/'// &
& 'Profiles'//coords_point//timestring//'.dat',status='unknown')
write (out_unit,*) '1 - depth, m'
write (out_unit,*) '2 - temperature, C'
write (out_unit,*) '3 - salinity, kg/kg'
write (out_unit,*) '4 - turbulent kinetic energy, m**2/s**2'
write (out_unit,*) '5 - dissipation rate, m**2/s**3'
write (out_unit,*) '6 - kinematic heat conductance, m**2/s'
write (out_unit,*) '7 - modulus of horizontal velocity, m/s'
write (out_unit,*) '8 - Richardson number, n/d'
do i = 1, M+1
write (out_unit,'(f12.6,f11.5,6e12.4)') &
& -dzeta_int(i)*h1, accum_var(1:nvars,ix,iy,i)
enddo
close(out_unit)
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy)/tsteps(ix,iy)
read(year1,*) work1
read(month1,*) work2
read(day1,*) work3
write (out_unita(ix,iy),fmt=*) & !, fmt = formatline) &
& int(work1), int(work2), int(work3), accum_var_scalar(1:nvars_scalar,ix,iy)
day_old(ix,iy) = day
accum_var(:,ix,iy,:) = 0.d0
accum_var_scalar(:,ix,iy) = 0.d0
tsteps(ix,iy) = 0.d0
endif
if (firstcall) firstcall=.false.
if (firstcallixiy(ix,iy)) firstcallixiy(ix,iy) = .false.
END SUBROUTINE DAY_OUT
SUBROUTINE HOUR_OUT(ix,iy,nx,ny,year,month,day,hour,time)
use ARRAYS_WATERSTATE, only : Tw1, Sal1
use ARRAYS_TURB!, only : E1, eps1, row, H_mixed_layer, H_entrainment, Wedderburn
use ARRAYS_BATHYM, only : h1
use ARRAYS_GRID, only : dzeta_int, dzeta_05int
use ARRAYS, only : u1, v1
use TIMEVAR, only : hour_sec
use ARRAYS_OXYGEN, only: oxyg, DIC
use ARRAYS_METHANE, only : qwater
use PHYS_FUNC, only : HC_CORR_CARBEQUIL
use PHYS_CONSTANTS, only : Kelvin0
use METH_OXYG_CONSTANTS, only : &
& molm3tomgl_o2, molm3tomcM, molm3tomcM, pH
use DRIVING_PARAMS, only : M, scale_output
use INOUT_PARAMETERS, only : &
& lake_hour_out_unit_min, &
& lake_hour_out_unit_max
implicit none
! Input variables
! Reals
real(kind=ireals), intent(in) :: hour
real(kind=ireals), intent(in) :: time
! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: year, month, day
! Local variables
! Reals
real(kind=ireals), allocatable :: accum_var (:,:,:,:)
real(kind=ireals), allocatable :: var (:,:,:,:)
real(kind=ireals), allocatable :: var_scalar(:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar(:,:,:)
real(kind=ireals), external :: DZETA
! Integers
integer(kind=iintegers), parameter :: n_var = 30
integer(kind=iintegers), parameter :: n_var_scalar = 10
integer(kind=iintegers), allocatable :: hour_old(:,:)
integer(kind=iintegers), allocatable :: tsteps(:,:)
integer(kind=iintegers), allocatable :: n_unit(:,:)
integer(kind=iintegers) :: i
integer(kind=iintegers) :: out_unit = lake_hour_out_unit_min
character :: month1*2
character :: year1*4
character :: day1*2
character :: hour1*2
character :: timestring*10
character :: coords_point*6
character :: format_char*100
logical :: firstcall = .true.
logical, allocatable :: firstcall_ixiy(:,:)
SAVE
if (firstcall) then
allocate (accum_var (n_var,1:nx,1:ny,1:M+1) )
allocate (var (n_var,1:nx,1:ny,1:M+1) )
allocate (accum_var_scalar (n_var_scalar,1:nx,1:ny) )
allocate (var_scalar (n_var_scalar,1:nx,1:ny) )
allocate (tsteps (1:nx,1:ny) )
allocate (hour_old (1:nx,1:ny) )
allocate (n_unit(1:nx,1:ny))
allocate (firstcall_ixiy(1:nx,1:ny))
accum_var(:,:,:,:) = 0.d0
accum_var_scalar(:,:,:) = 0.d0
var(:,:,:,:) = 0.d0
var_scalar(:,:,:) = 0.d0
tsteps(:,:) = 0
hour_old(:,:) = int(hour)
firstcall_ixiy(:,:) = .true.
n_unit(:,:) = out_unit + 1
endif
var(1,ix,iy,1:M) = PEMF (1:M)
var(2,ix,iy,1:M+1) = PDENS_DOWN (1:M+1)
var(3,ix,iy,1:M+1) = PT_DOWN (1:M+1)
var(4,ix,iy,1:M+1) = PSAL_DOWN (1:M+1)
var(5,ix,iy,1:M) = k_turb_T_flux (1:M)
var(6,ix,iy,1:M) = T_massflux (1:M)
var(7,ix,iy,1:M+1) = row (1:M+1)
var(8,ix,iy,1:M+1) = Tw1 (1:M+1)
var(9,ix,iy,1:M+1) = Sal1 (1:M+1)
var(10,ix,iy,1:M) = E1 (1:M)
var(11,ix,iy,1:M) = eps1 (1:M)
if (scale_output%par == 1) then
var(12,ix,iy,1) = H_mixed_layer ! Scale
var(13,ix,iy,1) = Buoyancy0 ! Scale
var(14,ix,iy,1) = w_conv_scale ! Scale
var(15,ix,iy,1) = T_conv_scale ! Scale
elseif (scale_output%par == 0) then
var(12,ix,iy,1) = 1.d0 ! No scaling
var(13,ix,iy,1) = 1.d0 ! No scaling
var(14,ix,iy,1) = 1.d0 ! No scaling
var(15,ix,iy,1) = 1.d0 ! No scaling
else
print*, 'Scale_output is ', scale_output%par, &
& ', must be 0 or 1: STOP'
STOP
endif
var(16,ix,iy,1:M) = S (1:M)
var(17,ix,iy,1:M) = Gen (1:M)
var(18,ix,iy,1:M) = Gen_seiches (1:M)
var(19,ix,iy,1:M) = TKE_turb_trans (1:M)
var(20,ix,iy,2:M) = k3_mid (2:M)
var(21,ix,iy,1:M+1) = u1(1:M+1)
var(22,ix,iy,1:M+1) = v1(1:M+1)
var(23,ix,iy,1:M+1) = sqrt(u1(1:M+1)*u1(1:M+1) + v1(1:M+1)*v1(1:M+1))
!Biogeochemical substances
var(24,ix,iy,1:M+1) = oxyg(1:M+1,1)*molm3tomgl_o2
do i = 1, M+1
var(25,ix,iy,i) = DIC(i,1) / HC_CORR_CARBEQUIL(Tw1(i)+Kelvin0,pH) * molm3tomcM
enddo
var(26,ix,iy,1:M+1) = qwater(1:M+1,1)*molm3tomcM
!Coefficients of turbulent viscosity
var(27,ix,iy,1:M) = k2 (1:M )
var(28,ix,iy,1:M) = viscturb_ziriv (1:M )
var(29,ix,iy,1:M) = eps_ziriv (1:M )
var(30,ix,iy,1:M+1) = speed_ziriv (1:M+1)
var_scalar(1,ix,iy) = S_integr_positive
var_scalar(2,ix,iy) = S_integr_negative
var_scalar(3,ix,iy) = Gen_integr
var_scalar(4,ix,iy) = eps_integr
var_scalar(5,ix,iy) = TKE_balance
var_scalar(6,ix,iy) = H_entrainment
var_scalar(7,ix,iy) = TKE_turb_trans_integr
var_scalar(8,ix,iy) = Wedderburn
var_scalar(9,ix,iy) = LakeNumber
var_scalar(10,ix,iy) = Rossby_rad
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + var_scalar(:,ix,iy)
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) + var(:,ix,iy,:)
tsteps(ix,iy) = tsteps(ix,iy) + 1
if (int(hour)/=hour_old(ix,iy)) then
accum_var(:,ix,iy,:) = accum_var(:,ix,iy,:) / real(tsteps(ix,iy))
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) / real(tsteps(ix,iy))
call DATEMINUS(3,year,month,day,hour,year1,month1,day1,hour1)
call TIMESTR(10,year1,month1,day1,hour1,timestring)
call CHECK_UNIT(lake_hour_out_unit_min,lake_hour_out_unit_max,out_unit)
write (coords_point, '(2i3)') ix, iy
! Writing to the file Profiles<yyyymmddhh>.dat
open(out_unit, file=outpath(1:len_trim(outpath))// &
& 'hourly/'//'Profiles'//coords_point//timestring//'.dat', &
& status='unknown')
write (out_unit,*) '1 - depth, normalized'
write (out_unit,*) '2 - temperature, normalized'
write (out_unit,*) '3 - salinity, kg/kg'
write (out_unit,*) '4 - water density, kg/m**3'
write (out_unit,*) '5 - turbulent kinetic energy, normalized'
write (out_unit,*) '6 - dissipation rate, normalized'
write (out_unit,*) '7 - eddy viscosity for TKE, m**2/s'
write (out_unit,*) '8 - mass flux, m/s'
write (out_unit,*) '9 - downdraft temperature, C'
write (out_unit,*) '10 - downdraft salinity, kg/kg'
write (out_unit,*) '11 - downdraft density, kg/m**3'
write (out_unit,*) '12 - x-component of speed, m/s'
write (out_unit,*) '13 - y-component of speed, m/s'
write (out_unit,*) '14 - modulus of speed, m/s'
write (out_unit,*) '15 - oxygen, mg/l'
write (out_unit,*) '16 - co2, mcmol/l'
write (out_unit,*) '17 - ch4, mcmol/l'
write (out_unit,*) '18 - eddy viscosity from k-epsilon model, m**2/s'
write (out_unit,*) '19 - eddy viscosity from Ziriv model, m**2/s'
write (out_unit,*) '20 - dissipation rate from Ziriv model, m**2/s**3'
write (out_unit,*) '21 - longitudinal speed from Ziriv model, m/s'
do i=1,M
write (out_unit,'(f14.6,f12.5,e10.3,f10.3,3e10.3,4e16.7,3f11.5,3e16.7,3e10.3,f11.5)') &
& -dzeta_int(i)*h1 /accum_var(12,ix,iy,1), &
& (accum_var(8,ix,iy,i) - scale_output%par*maxval(accum_var(8,ix,iy,:)) ) &
& /accum_var(15,ix,iy,1), &
& accum_var(9,ix,iy,i), &
& accum_var(7,ix,iy,i), &
& accum_var(10,ix,iy,i) / (accum_var(13,ix,iy,1)*accum_var(12,ix,iy,1) / &
& accum_var(14,ix,iy,1) ), &
& accum_var(11,ix,iy,i) / accum_var(13,ix,iy,1), &
& accum_var(20,ix,iy,i), & !accum_var(10,i)**2 / accum_var(11,i),
& accum_var(1,ix,iy,i), &
& accum_var(3,ix,iy,i), &
& accum_var(4,ix,iy,i), &
& accum_var(2,ix,iy,i), &
& accum_var(21,ix,iy,i), &
& accum_var(22,ix,iy,i), &
& accum_var(23,ix,iy,i), &
& accum_var(24,ix,iy,i), &
& accum_var(25,ix,iy,i), &
& accum_var(26,ix,iy,i), &
& accum_var(27,ix,iy,i), &
& accum_var(28,ix,iy,i), &
& accum_var(29,ix,iy,i), &
& accum_var(30,ix,iy,i)
enddo
close (out_unit)
! Writing to the file EDMF_profiles<yyyymmddhh>.dat
open(out_unit, file=outpath(1:len_trim(outpath))// &
& 'hourly/'//'EDMF_profiles'//coords_point//timestring//'.dat', &
& status='unknown')
write (out_unit,*) '1 - depth, normalized by mixed layer depth'
write (out_unit,*) '2 - "k - flux" of temperature, normalized'
write (out_unit,*) '3 - mass flux of temperature, normalized'
write (out_unit,*) '4 - total flux of temperature, normalized'
write (out_unit,'( 4(i10,6x) )') 1,2,3,4
do i=1, M
write (out_unit,'(f18.6,3e16.7)') &
& -dzeta_05int(i)*h1/ accum_var(12,ix,iy,1), &
& accum_var(5,ix,iy,i)/(accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1)), &
& accum_var(6,ix,iy,i)/(accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1)), &
& (accum_var(5,ix,iy,i) + accum_var(6,ix,iy,i) ) / &
& (accum_var(14,ix,iy,1)*accum_var(15,ix,iy,1))
enddo
close (out_unit)
! Writing to the file TKE_budget<yyyymmddhh>.dat
open(out_unit, file=outpath(1:len_trim(outpath))// &
& 'hourly/'//'TKE_budget'//coords_point//timestring//'.dat', &
& status='unknown')
write (out_unit,*) '1 - depth, normalised with mixed layer depth'
write (out_unit,*) '2 - shear production, normalised with surface buoyancy flux'
write (out_unit,*) '3 - seiches production, normalised with surface buoyancy flux'
write (out_unit,*) '4 - buoyancy source, normalised with surface buoyancy flux'
write (out_unit,*) '5 - dissipation rate, normalised with surface buoyancy flux'
write (out_unit,*) '6 - turbulent transport, normalised with surface buoyancy flux'
write (out_unit,'( 6(i10,6x) )') 1,2,3,4,5,6
do i = 1, M
write (out_unit,'(f18.6,5e16.7)') &
& -dzeta_05int(i)*h1/ accum_var(12,ix,iy,1), &
& accum_var(17,ix,iy,i) / accum_var(13,ix,iy,1), &
& accum_var(18,ix,iy,i) / accum_var(13,ix,iy,1), &
& accum_var(16,ix,iy,i) / accum_var(13,ix,iy,1), &
& accum_var(11,ix,iy,i) / accum_var(13,ix,iy,1), &
& accum_var(19,ix,iy,i) / accum_var(13,ix,iy,1)
enddo
! write (out_unit, '(a)')
close (out_unit)
if (firstcall_ixiy(ix,iy)) then
call CHECK_UNIT(lake_hour_out_unit_min,lake_hour_out_unit_max,n_unit(ix,iy))
open(n_unit(ix,iy), file=outpath(1:len_trim(outpath))// &
& 'hourly/'//'TKE_integr_conv'//coords_point//timestring//'.dat', &
& status='unknown')
write (n_unit(ix,iy),*) 'Col. 1 - year'
write (n_unit(ix,iy),*) 'Col. 2 - month'
write (n_unit(ix,iy),*) 'Col. 3 - day'
write (n_unit(ix,iy),*) 'Col. 4 - hour'
write (n_unit(ix,iy),*) 'Col. 5 - the time from the start of integration, hours'
write (n_unit(ix,iy),*) 'Col. 6 - H_entrainment, m '
write (n_unit(ix,iy),*) 'Col. 7 - B_integr+, m**3/s**3'
write (n_unit(ix,iy),*) 'Col. 8 - B_integr-, m**3/s**3'
write (n_unit(ix,iy),*) 'Col. 9 - S_integr, m**3/s**3'
write (n_unit(ix,iy),*) 'Col. 10 - eps_integr, m**3/s**3'
write (n_unit(ix,iy),*) 'Col. 11 - TKE turb trans integr, m**3/s**3'
write (n_unit(ix,iy),*) 'Col. 12 - TKE_balance, m**3/s**3'
write (n_unit(ix,iy),*) 'Col. 13 - Wedderburn number, n/d'
write (n_unit(ix,iy),*) 'Col. 14 - Lake number, n/d'
write (n_unit(ix,iy),*) 'Col. 15 - Internal Rossby radius, m'
firstcall_ixiy(ix,iy) = .false.
endif
format_char = '(3i7,f8.2,f13.2,f10.2,9e15.6)'
write (n_unit(ix,iy), format_char) &
& year,month,day,hour,time/hour_sec, &
& accum_var_scalar(6,ix,iy), &
& accum_var_scalar(1,ix,iy), &
& accum_var_scalar(2,ix,iy), &
& accum_var_scalar(3,ix,iy), &
& accum_var_scalar(4,ix,iy), &
& accum_var_scalar(7,ix,iy), &
& accum_var_scalar(5,ix,iy), &
& accum_var_scalar(8,ix,iy), &
& accum_var_scalar(9,ix,iy), &
& accum_var_scalar(10,ix,iy)
hour_old(ix,iy) = int(hour)
accum_var(:,ix,iy,:) = 0.d0
accum_var_scalar(:,ix,iy) = 0.d0
tsteps(ix,iy) = 0
endif
if (firstcall) firstcall = .false.
!if (firstcall_ixiy(ix,iy)) firstcall_ixiy(ix,iy) = .false.
END SUBROUTINE HOUR_OUT
SUBROUTINE EVERYSTEP_OUT(ix,iy,nx,ny)
use ARRAYS!, only : time, nstep
use ARRAYS_WATERSTATE!, only : Tw1, Sal1
use ARRAYS_BATHYM!, only : h1
use ARRAYS_GRID!, only : dzeta_int, M
use ARRAYS_TURB!, only : &
! & S_integr_positive, &
! & S_integr_negative, &
! & Gen_integr, &
! & eps_integr, &
! & TKE_balance, &
! & E_integr, &
! & H_entrainment, &
! & E1, eps1, &
! & k_turb_T_flux
use DRIVING_PARAMS, only : everystep, M
use INOUT_PARAMETERS, only : &
& lake_everystep_out_unit_min, &
& lake_everystep_out_unit_max
implicit none
! Input variables
! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
! Local variables
! Reals
real(kind=ireals), external :: DZETA
real(kind=ireals), parameter :: ACC = 1.d-20
real(kind=ireals), parameter :: hour_sec = 60.*60.
! Integers
integer(kind=iintegers), allocatable :: n_unit(:,:,:)
integer(kind=iintegers) :: i ! Loop index
! Characters
character :: coords_point*6
character :: format_char*100
! Logicals
logical, allocatable :: firstcall(:,:)
SAVE
if (.not.allocated(firstcall)) then
allocate (firstcall(1:nx, 1:ny) )
allocate (n_unit (1:2, 1:nx, 1:ny) )
firstcall = .true.
n_unit = lake_everystep_out_unit_min
endif
if (firstcall(ix,iy)) then
write (coords_point, '(2i3)') ix, iy
if (everystep%par /= 2) then
call CHECK_UNIT(lake_everystep_out_unit_min,lake_everystep_out_unit_max, &
& n_unit(1,ix,iy))
open (n_unit(1,ix,iy), file=outpath(1:len_trim(outpath))//'everystep/'// &
& 'Profiles'//coords_point//'.dat', status='unknown')
write (n_unit(1,ix,iy),*) '1 - depth, m'
write (n_unit(1,ix,iy),*) '2 - temperature, C'
write (n_unit(1,ix,iy),*) '3 - salinity, kg/kg'
write (n_unit(1,ix,iy),*) '4 - turbulent kinetic energy, m**2/s**2'
write (n_unit(1,ix,iy),*) '5 - disspation rate, m**2/s**3'
write (n_unit(1,ix,iy),*) '6 - eddy diffusivity (TKE**2/dissipation), m**2/s'
write (n_unit(1,ix,iy),*) '7 - k-flux of temperature, K*m/s'
endif
call CHECK_UNIT(lake_everystep_out_unit_min,lake_everystep_out_unit_max, &
& n_unit(2,ix,iy))
open (n_unit(2,ix,iy), file=outpath(1:len_trim(outpath))//'everystep/'// &
& 'TKE_integr'//coords_point//'.dat', status='unknown')
write (n_unit(2,ix,iy),*) '1 - timestep '
write (n_unit(2,ix,iy),*) '2 - time, hours '
write (n_unit(2,ix,iy),*) '3 - B_integr+, m**3/s**3'
write (n_unit(2,ix,iy),*) '4 - B_integr-, m**3/s**3'
write (n_unit(2,ix,iy),*) '5 - S_integr, m**3/s**3'
write (n_unit(2,ix,iy),*) '6 - eps_integr, m**3/s**3'
write (n_unit(2,ix,iy),*) '7 - TKE_balance, m**3/s**3'
write (n_unit(2,ix,iy),*) '8 - E_integr, m**3/s**2'
endif
if (everystep%par /= 2) then
write (n_unit(1,ix,iy),*) 'nstep = ', nstep
format_char = '(f10.6,f9.3,5e12.4)'
do i=1,M
write (n_unit(1,ix,iy), format_char) &
& -dzeta_int(i)*h1,Tw1(i),Sal1(i),E1(i),eps1(i), &
& E1(i)**2/(eps1(i)+ACC), k_turb_T_flux(i)
enddo
endif
format_char = '(i10, 2f10.2, 6e15.6)'
write(n_unit(2,ix,iy), format_char) &
& nstep, time/hour_sec, H_entrainment, S_integr_positive, S_integr_negative, &
& Gen_integr, eps_integr, TKE_balance, E_integr
if (firstcall(ix,iy)) firstcall(ix,iy)=.false.
END SUBROUTINE EVERYSTEP_OUT
SUBROUTINE SERIES_OUT(ix,iy,nx,ny,year,month,day,hour,tsw)
use DRIVING_PARAMS, only : M, ns, Mice, dt_out
use NUMERIC_PARAMS , only : &
& ms, small_value
use ARRAYS_TURB, only : row, H_mixed_layer, H_entrainment, &
& signwaveheight, w_conv_scale, u_star, Ri_bulk, maxN, i_maxN, &
& ThermThick, ReTherm, RiTherm, TKE_budget_terms
use ARRAYS_WATERSTATE, only : Tw1, Ti1, lamw, salice
use ARRAYS_GRID, only : nsoilcols
use ARRAYS_BATHYM, only : &
& h1,hx1,hx2,hy1,hy2, &
& hx1t,hx2t,hy1t,hy2t, &
& l1,hs1,ls1, &
& voldef, vol, &
& area_int, bathymwater
use ARRAYS_METHANE, only : &
& fbbleflx_ch4_sum, &
& fbbleflx_ch4, &
& febul0, fplant, &
& fdiff_lake_surf, &
& fdiffbot, foutflow, &
& qwater, qsoil, &
& rprod_total_newC, &
& rprod_total_oldC, &
& h_talik, lammeth, &
& qwateroxidtot, qwateroxidML
use ARRAYS_OXYGEN, only : &
& fbbleflx_co2_sum, &
& fbbleflx_co2, &
& fbbleflx_o2_sum, &
& fbbleflx_o2, &
& oxyg, &
& fco2, fo2
use ARRAYS_SOIL, only : Tsoil1, lsm
use ARRAYS, only: &
& Tskin, &
& zinv, &
& time, &
& T_0dim, &
& roughness,emissivity,albedo,aM,bM,relhums
use ATMOS, only: &
& eflux0_kinem, &
& eflux0, &
& turb_density_flux, &
& hw, xlew, botflux, &
& tau, windwork, surfrad, &
& discharge, velfrict_bot, &
& shortwave
use INOUT_PARAMETERS, only : &
& lake_series_out_unit_min, &
& lake_series_out_unit_max
use METH_OXYG_CONSTANTS, only : &
& molmass_ch4, mfs
use TIMEVAR, only : &
hour_sec, day_sec, hour_min
use PHYS_CONSTANTS, only : &
Kelvin0, cw_m_row0, roa0, sabs, row0
use DZETA_MOD, only : VARMEAN
use ARRAYS_GRID, only : gsp
use SNOWSOIL_MOD, only : Tsn, itop
implicit none
! Input variables
! Reals
real(kind=ireals), intent(in) :: tsw
real(kind=ireals), intent(in) :: hour
! Integers
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: year, month, day
! Local variables
! Reals
!real(kind=ireals), parameter :: mfs = molmass_ch4*8.64d+7 ! the transform multiplier for methane flux
! from mol/(m**2*s) to mg/(m**2*day)
real(kind=ireals), parameter :: small_number = 1.d-5
real(kind=ireals) :: T_mean, voldef0 = 0.
real(kind=ireals) :: methfluxML, methfluxshalsed
real(kind=ireals), allocatable :: tsteps(:,:)
real(kind=ireals), allocatable :: var_scalar (:,:,:)
real(kind=ireals), allocatable :: accum_var_scalar(:,:,:)
! Integers
integer(kind=iintegers), parameter :: nfiles = 6
integer(kind=iintegers), parameter :: numaccum = 5 ! The number of scalar variables being accumulated (averaged)
integer(kind=iintegers), allocatable :: n_unit(:,:,:)
integer(kind=iintegers), allocatable :: count_out(:,:)
! Characters
character :: coords_point*6
character :: format_char*100, workchar*10
! Logicals
logical, allocatable :: firstcall(:,:)
!real(kind=ireals) :: dz(ms)
integer(kind=iintegers) :: i
!common /SOILDAT/ dz,itop
!real(kind=ireals) :: Tsn(1:ms), cs(1:ms)
!common /SNOW_CHAR/ Tsn,cs
SAVE
if (.not.allocated(firstcall)) then
allocate (firstcall(1:nx, 1:ny) )
allocate (count_out(1:nx, 1:ny) )
allocate (n_unit (1:nfiles, 1:nx, 1:ny) )
allocate (tsteps (1:nx, 1:ny))
allocate (var_scalar (1:numaccum, 1:nx, 1:ny) )
allocate (accum_var_scalar (1:numaccum, 1:nx, 1:ny) )
tsteps(:,:) = 0.d0
var_scalar(:,:,:) = 0.d0
accum_var_scalar(:,:,:) = 0.d0
firstcall(:,:) = .true.
n_unit = lake_series_out_unit_min
endif
if (firstcall(ix,iy)) then
write (coords_point, '(2i3)') ix, iy
i = 1
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(i,ix,iy))
open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
& 'layers'//coords_point//'.dat', status='unknown')
write (n_unit(i,ix,iy),*)'Col. 1 - year'
write (n_unit(i,ix,iy),*)'Col. 2 - month'
write (n_unit(i,ix,iy),*)'Col. 3 - day'
write (n_unit(i,ix,iy),*)'Col. 4 - hour'
write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
write (n_unit(i,ix,iy),*)'Col. 6 - water layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 7 - W mixed layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 8 - E mixed layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 9 - S mixed layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 10 - N mixed layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 11 - W lower layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 12 - E lower layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 13 - S lower layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 14 - N lower layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 15 - ice layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 16 - snow layer thickness, m'
write (n_unit(i,ix,iy),*)'Col. 17 - bottom ice thickness, m'
write (n_unit(i,ix,iy),*)'Col. 18 - reservoir volume, m**3'
write (n_unit(i,ix,iy),*)'Col. 19 - volume deficit (accumulated), m**3'
i = i + 1
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(i,ix,iy))
open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
& 'T_fluxes'//coords_point//'.dat',status='unknown')
write (n_unit(i,ix,iy),*)'Col. 1 - year'
write (n_unit(i,ix,iy),*)'Col. 2 - month'
write (n_unit(i,ix,iy),*)'Col. 3 - day'
write (n_unit(i,ix,iy),*)'Col. 4 - hour'
write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
write (n_unit(i,ix,iy),*)'Col. 6 - surface temperature, C'
write (n_unit(i,ix,iy),*)'Col. 7 - water skin temperature, C'
write (n_unit(i,ix,iy),*)'Col. 8 - water surface temperature, C'
write (n_unit(i,ix,iy),*)'Col. 9 - mean temperature of water coloumn, C'
write (n_unit(i,ix,iy),*)'Col. 10 - maximal temperature in the water coloumn, C'
write (n_unit(i,ix,iy),*)'Col. 11 - zero-dimensional model temperature, C'
write (n_unit(i,ix,iy),*)'Col. 12 - upper ice surface temperature, C'
write (n_unit(i,ix,iy),*)'Col. 13 - upper snow surface temperature, C'
write (n_unit(i,ix,iy),*)'Col. 14 - sensible heat flux, W/m**2'
write (n_unit(i,ix,iy),*)'Col. 15 - latent heat flux, W/m**2'
write (n_unit(i,ix,iy),*)'Col. 16 - downward heat flux at the upper lake surface, W/m**2'
write (n_unit(i,ix,iy),*)'Col. 17 - downward heat flux at the lake bottom, W/m**2'
write (n_unit(i,ix,iy),*)'Col. 18 - friction velocity at the surface (waterside), m/s'
write (n_unit(i,ix,iy),*)'Col. 19 - friction velocity at the bottom, m/s'
write (n_unit(i,ix,iy),*)'Col. 20 - wind work at the water surface, W/m**2'
write (n_unit(i,ix,iy),*)'Col. 21 - albedo of the lake-atmosphere interface, n/d'
write (n_unit(i,ix,iy),*)'Col. 22 - shortwave radiation penetrated below surface, W/m**2'
write (n_unit(i,ix,iy),*)'Col. 23 - significant wave height, m'
write (n_unit(i,ix,iy),*)'Col. 24 - bottom ice salinity, kg/kg'
write (n_unit(i,ix,iy),*)'Col. 25 - discharge in x direction, m**3/s'
write (n_unit(i,ix,iy),*)'Col. 26 - discharge in y direction, m**3/s'
write (n_unit(i,ix,iy),*)'Cols. 27... - top soil columns temperature, m'
i = i + 1
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(i,ix,iy))
open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
& 'conv_series'//coords_point//'.dat',status='unknown')
write (n_unit(i,ix,iy),*)'Col. 1 - year'
write (n_unit(i,ix,iy),*)'Col. 2 - month'
write (n_unit(i,ix,iy),*)'Col. 3 - day'
write (n_unit(i,ix,iy),*)'Col. 4 - hour'
write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
write (n_unit(i,ix,iy),*)'Col. 6 - surface temperature, C'
write (n_unit(i,ix,iy),*)'Col. 7 - surface density, kg/m**3'
write (n_unit(i,ix,iy),*)'Col. 8 - heat flux downwards at the surface, K*m/s'
write (n_unit(i,ix,iy),*)'Col. 9 - turbulent density flux at the surface, kg/m**2/s'
write (n_unit(i,ix,iy),*)'Col. 10 - inversion depth (mass flux diagnostics), m'
write (n_unit(i,ix,iy),*)'Col. 11 - mixed layer depth, m'
write (n_unit(i,ix,iy),*)'Col. 12 - entrainment depth, m'
write (n_unit(i,ix,iy),*)'Col. 13 - convective velocity scale, m/s'
write (n_unit(i,ix,iy),*)'Col. 14 - friction velocity, m/s'
write (n_unit(i,ix,iy),*)'Col. 15 - -w_star/u_star, m/s'
write (n_unit(i,ix,iy),*)'Col. 16 - bulk Richardson number, n/d'
write (n_unit(i,ix,iy),*)'Col. 17 - maximal Brunt-Vaisala frequency, s**(-1)'
write (n_unit(i,ix,iy),*)'Col. 18 - minimal thermal conductance, m**2/s'
write (n_unit(i,ix,iy),*)'Col. 19 - thermocline thickness, m'
write (n_unit(i,ix,iy),*)'Col. 20 - Reynolds number in thermocline , n/d'
write (n_unit(i,ix,iy),*)'Col. 21 - bulk Richardson number in thermocline, n/d'
i = i + 1
! The output of LakeMIP file 1 format
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(i,ix,iy))
open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
& 'LakeMIP_file1'//coords_point//'.dat',status='unknown')
! write (n_unit(i,ix,iy),*)'Col. 1 - time, days'
! write (n_unit(i,ix,iy),*)'Col. 1 - year'
! write (n_unit(i,ix,iy),*)'Col. 2 - month'
! write (n_unit(i,ix,iy),*)'Col. 3 - day in a month'
! write (n_unit(i,ix,iy),*)'Col. 4 - hour'
! write (n_unit(i,ix,iy),*)'Col. 5 - minute'
! write (n_unit(i,ix,iy),*)'Col. 6 - mixed-layer temperature &
! &(equals to the temperature of the upper water surface), degrees Celsius'
! write (n_unit(i,ix,iy),*)'Col. 7 - mean temperature of the water column, degrees Celsius'
! write (n_unit(i,ix,iy),*)'Col. 8 - bottom temperature, degrees Celsius'
! write (n_unit(i,ix,iy),*)'Col. 9 - mixed-layer depth, meters'
! write (n_unit(i,ix,iy),*)'Col. 10 - ice thickness, meters'
! write (n_unit(i,ix,iy),*)'Col. 11 - snow thickness, meters'
! write (n_unit(i,ix,iy),*)'Col. 12 - temperature at the ice upper surface, degrees Celsius'
! write (n_unit(i,ix,iy),*)'Col. 13 - temperature at the snow upper surface, degrees Celsius'
! write (n_unit(i,ix,iy),*)'Col. 14 - sensible heat flux at the lake-atmosphere &
! &interface, averaged over output interval, upwards, W/m**2'
! write (n_unit(i,ix,iy),*)'Col. 15 - latent heat flux at the lake-atmosphere &
! &interface, averaged over output interval, upwards, W/m**2'
! write (n_unit(i,ix,iy),*)'Col. 16 - momentum flux at the lake-atmosphere &
! &interface, averaged over output interval, positive, N/m**2'
! write (n_unit(i,ix,iy),*)'Col. 17 - upward long-wave radiation flux &
! &at the lake-atmosphere interface, averaged over output interval, W/m**2'
! write (n_unit(i,ix,iy),*)'Col. 18 - downward heat flux at the lake-atmosphere &
! &interface, averaged over output interval, W/m**2'
! write (n_unit(i,ix,iy),*)'Col. 19 - surface albedo, n/d'
i = i + 1
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(i,ix,iy))
open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
& 'methane_series'//coords_point//'.dat',status='unknown')
write (n_unit(i,ix,iy),*)'Col. 1 - year'
write (n_unit(i,ix,iy),*)'Col. 2 - month'
write (n_unit(i,ix,iy),*)'Col. 3 - day'
write (n_unit(i,ix,iy),*)'Col. 4 - hour'
write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
write (n_unit(i,ix,iy),*)'Col. 6 - the talik depth, m'
write (n_unit(i,ix,iy),*)'Col. 7 - lake surface methane concentration, mol/m**3'
write (n_unit(i,ix,iy),*)'Col. 8 - lake bottom methane concentration, mol/m**3'
write (n_unit(i,ix,iy),*)'Col. 9 - soil bottom methane concentration, mol/m**3'
write (n_unit(i,ix,iy),*)'Col. 10 - lake surface oxygen concentration, mol/m**3'
write (n_unit(i,ix,iy),*)'Col. 11 - lake bottom oxygen concentration, mol/m**3'
write (n_unit(i,ix,iy),*)'Col. 12 - total methane production due to young C decomposition, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 13 - total methane production due to old C decomposition, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 14 - methane ebullition flux at the surface, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 15 - methane plant-mediated flux at the lake bottom, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 16 - methane diffusion flux at the lake bottom, upwards, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 17 - methane turbulent flux at the lake surface, upwards, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 18 - methane ebullition flux at the surface, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 19 - methane plant-mediated flux at the lake bottom, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 20 - methane diffusion flux at the lake bottom, upwards, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 21 - methane turbulent flux at the lake surface, upwards, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 22 - methane turbulent flux at the bottom of mixed layer normalized by &
& surface area, upwards, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 23 - methane flux from sediments in the mixed layer normalized by &
& surface area, upwards, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 24 - methane bubble flux at the bottom of the mixed layer normalized by &
& surface area, upwards, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 25 - total methane oxidation in water normalized by &
& surface area, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 26 - methane oxidation in mixed layer normalized by &
& surface area, mg/(m**2*day)'
write (n_unit(i,ix,iy),*)'Col. 27 - co2 turbulent flux at the lake surface, upwards, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 28 - co2 ebullition flux at the surface, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 29 - oxygen turbulent flux at the lake surface, upwards, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 30 - oxygen ebullition flux at the surface, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 31 - methane flux through outlet, normalized by surface area, mol/(m**2*s)'
write (n_unit(i,ix,iy),*)'Col. 32..32+ns-1 - methane ebullition flux at the surface from different soil columns, mg/(m**2*day)'
i = i + 1
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(i,ix,iy))
open (n_unit(i,ix,iy),file=outpath(1:len_trim(outpath))//'time_series/'// &
& 'TKE_budget_terms_series'//coords_point//'.dat',status='unknown')
write (n_unit(i,ix,iy),*)'Col. 1 - year'
write (n_unit(i,ix,iy),*)'Col. 2 - month'
write (n_unit(i,ix,iy),*)'Col. 3 - day'
write (n_unit(i,ix,iy),*)'Col. 4 - hour'
write (n_unit(i,ix,iy),*)'Col. 5 - the time from the start of integration, hours'
write (n_unit(i,ix,iy),*)'Col. 6 - TKE tendency, m**2/s**3'
write (n_unit(i,ix,iy),*)'Col. 7 - TKE turbulent transport, m**2/s**3'
write (n_unit(i,ix,iy),*)'Col. 8 - TKE shear production, m**2/s**3'
write (n_unit(i,ix,iy),*)'Col. 9 - TKE buoyancy production/sink, m**2/s**3'
write (n_unit(i,ix,iy),*)'Col. 10 - TKE dissipation rate, m**2/s**3'
count_out(ix,iy) = int(time/(dt_out%par*hour_sec))
endif
! Accumulating variables
var_scalar(1,ix,iy) = hw
var_scalar(2,ix,iy) = xlew
var_scalar(3,ix,iy) = tau
var_scalar(4,ix,iy) = surfrad
var_scalar(5,ix,iy) = eflux0
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy) + var_scalar(:,ix,iy)
tsteps(ix,iy) = tsteps(ix,iy) + 1.d0
! Writing to files
if (int(time/(dt_out%par*hour_sec))>count_out(ix,iy) .or. &
& abs(time - (count_out(ix,iy) + 1)*dt_out%par*hour_sec) < small_number) then
T_mean = VARMEAN(Tw1,bathymwater,1_iintegers)
format_char = '(3i7,f10.4,f15.4,f10.4,8e12.4,3f10.4,1x,e12.4,1x,f10.1)'
write (n_unit(1,ix,iy),format_char) year,month,day ,hour, &
& time/hour_sec , &
& h1, hx1, hx2, hy1, hy2,&
& hx1t, hx2t, hy1t, hy2t,&
& l1, hs1, ls1, vol, &
& voldef - voldef0
voldef0 = voldef
format_char = '(3i7,f8.2,f13.2,15e14.5,2f11.2,f11.4,f6.2,2e14.5'
! Formatting for possible multiple soil columns in a lake
write(workchar,'(i2)') nsoilcols
format_char = format_char(1:len_trim(format_char))//','// &
& workchar (1:len_trim(workchar) )//'f11.4)'
write (n_unit(2,ix,iy),format_char) year,month,day,hour, &
& time/hour_sec , &
& tsw - Kelvin0,Tskin(1), &
& Tw1(1), &
& T_mean,maxval(Tw1(1:M+1)), &
& T_0dim,Ti1(1),Tsn(itop), &
& hw,xlew,eflux0,botflux,sqrt(tau/row0), &
& velfrict_bot, &
& windwork, &
& albedo,shortwave*(1.-sabs),signwaveheight, &
& salice(Mice+1), &
& discharge(1:2), &
& Tsoil1(1,1:nsoilcols)
format_char = '(3i7,f8.2,f13.2,2f10.4,2e15.5,5f10.4,7e15.5)'
write (n_unit(3,ix,iy),format_char) year,month,day,hour, &
& time/hour_sec, &
& Tw1(1),row(1),eflux0_kinem, &
& turb_density_flux, zinv, &
& H_mixed_layer, H_entrainment, &
& w_conv_scale, u_star, &
& min(max(- w_conv_scale/(u_star + small_value), &
& -1.e+3_ireals),1.e+3_ireals), &
& min(max(Ri_bulk, -1.e+3_ireals),1.e+3_ireals), &
& maxN, minval(lamw(max(i_maxN-3,1_iintegers):min(i_maxN+3,M)))/cw_m_row0, &
& ThermThick, ReTherm, RiTherm
format_char = '(15(e15.8,1x))'! '(19(e15.8,1x))' !'(f10.5,10f8.2,f10.4,3f8.2)'
accum_var_scalar(:,ix,iy) = accum_var_scalar(:,ix,iy)/tsteps(ix,iy)
write (n_unit(4,ix,iy),format_char) & !float(year), float(month), float(day), hour, &
& time/day_sec, & !max((hour - floor(hour+small_number))*hour_min,0.d0), &
& Tw1(1), T_mean, &
& Tw1(M+1), H_mixed_layer, &
& l1, hs1, Ti1(1), Tsn(itop), & !missing_value, missing_value, missing_value, missing_value,
& accum_var_scalar(1:numaccum,ix,iy), &
& albedo
accum_var_scalar(:,ix,iy) = 0.d0
tsteps(ix,iy) = 0.d0
! Methane flux at the bottom of mixed layer
if (h1 > 0.) then
methfluxML = lammeth(i_maxN)*(qwater(i_maxN+1,2) - qwater(i_maxN,2))/(gsp%ddz(i_maxN)*h1)*&
& bathymwater(i_maxN)%area_half/bathymwater(1)%area_int
! Methane flux from shallow sediments (i.e. sediments in the mixed layer)
methfluxshalsed = 0.
do i = i_maxN, 1, -1
methfluxshalsed = methfluxshalsed + bathymwater(i)%area_int * lsm%water(i) * gsp%ddz05(i-1)*h1
enddo
methfluxshalsed = methfluxshalsed/bathymwater(1)%area_int
else
methfluxML = 0.
methfluxshalsed = 0.
endif
format_char = '(3i7,f8.2,f13.2,f8.2,31e15.5)'
write (n_unit(5,ix,iy),format_char) year,month,day,hour, &
& time/hour_sec, h_talik, &
& qwater(1,1), qwater(M+1,1), &
& qsoil(ns,nsoilcols), oxyg(1,1), &
& oxyg(M+1,1), &
& rprod_total_newC(nsoilcols), &
& rprod_total_oldC(nsoilcols), &
& fbbleflx_ch4_sum(0)/area_int(1), fplant, &
& - fdiffbot(nsoilcols), fdiff_lake_surf, &
& fbbleflx_ch4_sum(0)/area_int(1)*mfs, fplant*mfs, &
& - fdiffbot(nsoilcols)*mfs, &
& fdiff_lake_surf*mfs, &
& methfluxML*mfs, methfluxshalsed*mfs, &
& fbbleflx_ch4_sum(i_maxN)/area_int(1)*mfs, &
& qwateroxidtot*mfs, qwateroxidML*mfs, &
& fco2, fbbleflx_co2_sum(0)/area_int(1), &
& fo2, fbbleflx_o2_sum(0) /area_int(1), &
& foutflow, &
& fbbleflx_ch4(0,1:nsoilcols)*febul0(1:nsoilcols)*mfs
! format_char = '(3i7,f8.2,f13.2,f8.2,24e15.5)' ! two-meth
! write (n_unit(5,ix,iy),format_char) year,month,day,hour, & ! two-meth
! & time/hour_sec, h_talik, & ! two-meth
! & qwater2(1,1,1:2), & ! two-meth
! & qwater2(M+1,1,1:2), & ! two-meth
! & qsoil2(ns,1:2), & ! two-meth
! & oxyg(1,1), & ! two-meth
! & oxyg(M+1,1), & ! two-meth
! & rprod_total_newC, & ! two-meth
! & rprod_total_oldC, & ! two-meth
! & febul2(1:2), fplant, & ! two-meth
! & - fdiff2(1:2), - fdiff_lake_surf2(1:2), & ! two-meth
! & febul2(1:2)*mfs, fplant*mfs,& ! two-meth
! & - fdiff2(1:2)*mfs, & ! two-meth
! & - fdiff_lake_surf2(1:2)*mfs ! two-meth
format_char = '(3i7,f8.2,f13.2,5E15.5)'
write (n_unit(6,ix,iy),format_char) year,month,day,hour, &
& time/hour_sec, &
& TKE_budget_terms(1:5)
count_out(ix,iy) = count_out(ix,iy) + 1
endif
if (firstcall(ix,iy)) firstcall(ix,iy)=.false.
END SUBROUTINE SERIES_OUT
SUBROUTINE PROFILE_OUTPUT &
& (ix, iy, nx, ny, &
& year, month, day, hour, &
& time, dt_out, &
& Profile, z_prof, nprof, &
& z_out, nout, &
& outpath1, filename, &
& unitindex, ndec, last, Profile2, Profile3)
! The subroutine WATER_TEMP_OUT writes the water temperature profiles
! with the timestep dt_out
use INOUT_PARAMETERS, only : &
& lake_series_out_unit_min, &
& lake_series_out_unit_max
use TIMEVAR, only : &
hour_sec, day_sec, hour_min
implicit none
! Input variables
! Reals
real(kind=ireals), intent(in) :: Profile(1:nprof) ! The profile
real(kind=ireals), intent(in) :: z_prof(1:nprof) ! The levels of the profile elements
real(kind=ireals), intent(in) :: hour ! The hour
real(kind=ireals), intent(in) :: z_out(1:nout) ! The output grid
real(kind=ireals), intent(in) :: time, dt_out ! The time from beginning of integration
! The timestep of output
! Integers
integer(kind=iintegers), intent(in) :: ix, iy ! The current coordinates of lake point
integer(kind=iintegers), intent(in) :: nx, ny ! The maximal coordinates of lake points
integer(kind=iintegers), intent(in) :: year, month, day ! Year, month and day
integer(kind=iintegers), intent(in) :: nprof ! The number of layers in water, (M+1) - number of levels
integer(kind=iintegers), intent(in) :: nout
integer(kind=iintegers), intent(in) :: unitindex
integer(kind=iintegers), intent(in) :: ndec ! Number of decimal digits
character(len=*), intent(in) :: outpath1
character(len=*), intent(in) :: filename
logical, intent(in) :: last
real(kind=ireals), intent(in), optional :: Profile2(1:nprof), Profile3(1:nprof) ! Optional profiles
! Local variables
! Reals
real(kind=ireals), parameter :: small_number = 1.d-5
real(kind=ireals), allocatable :: Profile_out(:,:)
! Integers
integer(kind=iintegers), parameter :: unitindex_max = 100
integer(kind=iintegers), allocatable, save :: n_unit(:,:,:)
integer(kind=iintegers), allocatable, save :: count_out(:,:)
integer(kind=iintegers) :: i, k !Loop index
integer(kind=iintegers) :: nvars
!integer(kind=iintegers) :: nform(1:2)
! Characters
character(len=6) :: coords_point
character(len=100) :: formats(1:2)
! Logicals
logical, allocatable, save :: firstcall(:,:)
logical :: flag
if (.not.allocated(firstcall)) then
allocate (firstcall(1:nx, 1:ny) )
allocate (count_out(1:nx, 1:ny) )
allocate (n_unit (1:unitindex_max, 1:nx, 1:ny) )
firstcall(:,:) = .true.
n_unit = lake_series_out_unit_min
endif
nvars = 1
if (present(Profile2)) nvars = 2
if (present(Profile2) .and. present(Profile3)) nvars = 3
if (firstcall(ix,iy)) then
write (coords_point, '(2i3)') ix, iy
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(unitindex,ix,iy))
open (n_unit(unitindex,ix,iy), &
& file=outpath1(1:len_trim(outpath1))//'time_series/'// &
& filename(1:len_trim(filename))//coords_point//'.dat', status='unknown')
write (n_unit(unitindex,ix,iy),*)'Col. 1 - year'
write (n_unit(unitindex,ix,iy),*)'Col. 2 - month'
write (n_unit(unitindex,ix,iy),*)'Col. 3 - day'
write (n_unit(unitindex,ix,iy),*)'Col. 4 - hour'
write (n_unit(unitindex,ix,iy),*) &
& 'Col. 5 - the time from the start of integration, days' !& 'Col. 5 - the time from the start of integration, days'
write (n_unit(unitindex,ix,iy),*) &
& 'Col. [6...] - Depth1, Var1, Depth2, Var2, ..., DepthN, VarN'
call CHECK_UNIT(lake_series_out_unit_min,lake_series_out_unit_max, &
& n_unit(unitindex+1,ix,iy))
open (n_unit(unitindex+1,ix,iy), &
& file=outpath1(1:len_trim(outpath1))//'time_series/'// &
& filename(1:len_trim(filename))//coords_point//'f2'//'.dat', status='unknown')
write (n_unit(unitindex+1,ix,iy),*)'Col. 1 - year'
write (n_unit(unitindex+1,ix,iy),*)'Col. 2 - month'
write (n_unit(unitindex+1,ix,iy),*)'Col. 3 - day'
write (n_unit(unitindex+1,ix,iy),*)'Col. 4 - hour'
write (n_unit(unitindex+1,ix,iy),*) &
& 'Col. 5 - the time from the start of integration, days' !& 'Col. 5 - the time from the start of integration, days'
write (n_unit(unitindex+1,ix,iy),*) &
& 'Col. 6 - Depth, m'
write (n_unit(unitindex+1,ix,iy),*) &
& 'Col. 7 - Var, m'
count_out(ix,iy) = int(time/(dt_out*hour_sec))
endif
if (int(time/(dt_out*hour_sec))>count_out(ix,iy) .or. &
& abs(time - (count_out(ix,iy) + 1)*dt_out*hour_sec) < small_number) then
if (nout > 0) then
allocate (Profile_out(1:nout,1:nvars))
call LININTERPOL (z_prof,Profile,nprof,z_out,Profile_out(1,1),nout,flag,.false.)
if (nvars >= 2) call LININTERPOL (z_prof,Profile2,nprof,z_out,Profile_out(1,2),nout,flag,.false.)
if (nvars == 3) call LININTERPOL (z_prof,Profile3,nprof,z_out,Profile_out(1,3),nout,flag,.false.)
k = nout
call DEFFORMAT
write (unit=n_unit(unitindex,ix,iy),fmt=formats(1)) year,month,day,hour, &
& time/day_sec, & !max((hour - floor(hour+small_number))*hour_min,0.d0), &
& (z_out(i),Profile_out(i,1:nvars), i = 1, nout)
do i = 1, nout
write (unit=n_unit(unitindex+1,ix,iy),fmt=formats(2)) year,month,day,hour, &
& time/day_sec, & !max((hour - floor(hour+small_number))*hour_min,0.d0), &
& -z_out(i),Profile_out(i,1:nvars)
enddo
deallocate (Profile_out)
else
allocate (Profile_out(1:nprof,1:nvars))
Profile_out(1:nprof,1) = Profile(1:nprof)
if (nvars >= 2) Profile_out(1:nprof,2) = Profile2(1:nprof)
if (nvars == 3) Profile_out(1:nprof,3) = Profile3(1:nprof)
k = nprof
call DEFFORMAT
write (unit=n_unit(unitindex,ix,iy),fmt=formats(1)) year,month,day,hour, &
& time/day_sec, & ! max((hour - floor(hour+small_number))*hour_min,0.d0), &
& (z_prof(i),Profile_out(i,1:nvars), i = 1, nprof)
do i = 1, nprof
write (unit=n_unit(unitindex+1,ix,iy),fmt=formats(2)) year,month,day,hour, &
& time/day_sec, & ! max((hour - floor(hour+small_number))*hour_min,0.d0), &
& -z_prof(i),Profile_out(i,1:nvars)
enddo
endif
if (last) count_out(ix,iy) = count_out(ix,iy) + 1
endif
if (firstcall(ix,iy) .and. last) firstcall(ix,iy) = .false.
contains
SUBROUTINE DEFFORMAT
implicit none
type charl
character(len=20) :: ch
integer(kind=iintegers) :: len
end type charl
type(charl) :: nc1, nc2, nc3, nv, ymdh
ymdh%ch = '3i7,f8.2,' !'3i7,f8.2,' ! year,month,day,hour
ymdh%len = len_trim(ymdh%ch)
write(nc1%ch,'(i3)') k ; nc1%len = len_trim(nc1%ch)
write(nv%ch, '(i3)') nvars ; nv%len = len_trim(nv%ch)
if (ndec >= 0) then
write(nc2%ch,'(i3)') 6 + ndec ; nc2%len = len_trim(nc2%ch)
write(nc3%ch,'(i3)') ndec ; nc3%len = len_trim(nc3%ch)
formats(1) = '('//ymdh%ch(1:ymdh%len)//'f13.4,'//nc1%ch(1:nc1%len)// &
& '(f10.2, '//nv%ch(1:nv%len)//'f'//nc2%ch(1:nc2%len)//'.'//nc3%ch(1:nc3%len)//'))'
formats(2) = '('//ymdh%ch(1:ymdh%len)//'f13.4,f10.2,'//nv%ch(1:nv%len)//'f'//nc2%ch(1:nc2%len)//'.'// &
& nc3%ch(1:nc3%len)//')'
else
formats(1) = '('//ymdh%ch(1:ymdh%len)//'f13.4,'//nc1%ch(1:nc1%len)//'(f10.2, '//nv%ch(1:nv%len)//'E23.10E3))'
formats(2) = '('//ymdh%ch(1:ymdh%len)//'f13.4,f10.2,'//nv%ch(1:nv%len)//'E23.10E3)'
endif
END SUBROUTINE DEFFORMAT
!100 format (3i7,f8.2,f13.4,<k>(f10.2, f<6+ndec>.<ndec>))
!101 format (3i7,f8.2,f13.4,f10.2,f<6+ndec>.<ndec>)
!102 format (3i7,f8.2,f13.4,<k>(f10.2, E15.7))
!103 format (3i7,f8.2,f13.4,f10.2,E15.7)
END SUBROUTINE PROFILE_OUTPUT
SUBROUTINE ACCUM_VAR &
& (dt, l1, febul, febul0, fdiff, fdiff_lake_surf, qwateroxidtot, &
& ice_meth_oxid_total, &
& rprod_total_newC, rprod_total_oldC, &
& febultot, febulbottot, fdifftot, fdiff_lake_surftot, &
& rprod_total_newC_integr, rprod_total_oldC_integr, &
& methoxidwat, ice_meth_oxid, add_to_winter)
use METH_OXYG_CONSTANTS, only : &
& molmass_ch4, mfs, mf
implicit none
! Input/output variables
real(kind=ireals), intent(in) :: dt
real(kind=ireals), intent(in) :: l1
real(kind=ireals), intent(in) :: febul, febul0
real(kind=ireals), intent(in) :: fdiff
real(kind=ireals), intent(in) :: fdiff_lake_surf
!real(kind=ireals), intent(in) :: febul(1:2) ! two-meth
!real(kind=ireals), intent(in) :: fdiff(1:2) ! two-meth
!real(kind=ireals), intent(in) :: fdiff_lake_surf(1:2) ! two-meth
real(kind=ireals), intent(in) :: rprod_total_newC
real(kind=ireals), intent(in) :: rprod_total_oldC
real(kind=ireals), intent(in) :: qwateroxidtot
real(kind=ireals), intent(in) :: ice_meth_oxid_total
logical, intent(inout) :: add_to_winter
real(kind=ireals), intent(inout) :: febultot(1:2), febulbottot(1:2)
real(kind=ireals), intent(inout) :: fdifftot(1:2)
real(kind=ireals), intent(inout) :: fdiff_lake_surftot(1:2)
real(kind=ireals), intent(inout) :: rprod_total_newC_integr(1:2)
real(kind=ireals), intent(inout) :: rprod_total_oldC_integr(1:2)
real(kind=ireals), intent(inout) :: methoxidwat(1:2)
real(kind=ireals), intent(inout) :: ice_meth_oxid
!real(kind=ireals), intent(inout) :: febultot(1:2,1:2) ! two-meth
!real(kind=ireals), intent(inout) :: fdifftot(1:2,1:2) ! two-meth
!real(kind=ireals), intent(inout) :: fdiff_lake_surftot(1:2,1:2) ! two-meth
!real(kind=ireals), intent(inout) :: rprod_total_newC_integr(1:2,1:2) ! two-meth
!real(kind=ireals), intent(inout) :: rprod_total_oldC_integr(1:2,1:2) ! two-meth
! Local variables
real(kind=ireals), parameter :: day_sec = 24.*60.*60.
integer(kind=iintegers) :: i ! loop index
if (l1 == 0. .and. (.not. add_to_winter)) then
i = 1
else
i = 2
if (add_to_winter) add_to_winter = .false.
endif
febultot(i) = febultot(i) + febul*mfs*dt/day_sec
febulbottot(i) = febulbottot(i) + febul0*mfs*dt/day_sec
fdifftot(i) = fdifftot(i) - fdiff*mfs*dt/day_sec
fdiff_lake_surftot(i) = fdiff_lake_surftot(i) - fdiff_lake_surf*mfs*dt/day_sec
methoxidwat(i) = methoxidwat(i) - qwateroxidtot*mfs*dt/day_sec
rprod_total_newC_integr(i) = rprod_total_newC_integr(i) + rprod_total_newC*mfs*dt/day_sec
rprod_total_oldC_integr(i) = rprod_total_oldC_integr(i) + rprod_total_oldC*mfs*dt/day_sec
ice_meth_oxid = ice_meth_oxid + ice_meth_oxid_total*mf
!do i = 1, 2 ! two-meth
! if (l1 == 0.) then ! two-meth
! febultot(1,i) = febultot(1,i) + febul(i)*mfs*dt/day_sec ! two-meth
! fdifftot(1,i) = fdifftot(1,i) - fdiff(i)*mfs*dt/day_sec ! two-meth
! fdiff_lake_surftot(1,i) = fdiff_lake_surftot(1,i) - fdiff_lake_surf(i)*mfs*dt/day_sec ! two-meth
! rprod_total_newC_integr(1,i) = rprod_total_newC_integr(1,i) + rprod_total_newC*mfs*dt/day_sec ! two-meth
! rprod_total_oldC_integr(1,i) = rprod_total_oldC_integr(1,i) + rprod_total_oldC*mfs*dt/day_sec ! two-meth
! else ! two-meth
! febultot(2,i) = febultot(2,i) + febul(i)*mfs*dt/day_sec ! two-meth
! fdifftot(2,i) = fdifftot(2,i) - fdiff(i)*mfs*dt/day_sec ! two-meth
! fdiff_lake_surftot(2,i) = fdiff_lake_surftot(2,i) - fdiff_lake_surf(i)*mfs*dt/day_sec ! two-meth
! rprod_total_newC_integr(2,i) = rprod_total_newC_integr(2,i) + rprod_total_newC*mfs*dt/day_sec ! two-meth
! rprod_total_oldC_integr(2,i) = rprod_total_oldC_integr(2,i) + rprod_total_oldC*mfs*dt/day_sec ! two-meth
! endif ! two-meth
!enddo ! two-meth
END SUBROUTINE ACCUM_VAR
!> Subroutine TEMPLOC calculates the temperature and other scalars in a specified location
!! of a lake given the vertical profile of horizontally-averaged field
!! and deviations of heights of constant-density layers
SUBROUTINE TEMPLOC(M,nvars,nvar,itherm,h,dt_out,time,hx1ml,hy1ml,Tw,gsp,gs,bathymwater,rtemp, &
& basename,firstcall_)
!Currently implemented for single-lake output only.
use NUMERIC_PARAMS, only : pi
use ARRAYS_GRID, only : gridspacing_type, gridsize_type
use ARRAYS_BATHYM, only : bathym, aki, bki
use INOUT_PARAMETERS, only : lake_misc_unit_min, lake_misc_unit_max
use TIMEVAR, only : hour_sec
use DRIVING_PARAMS, only : grarr1
implicit none
!Input variables
integer(kind=iintegers), intent(in) :: M !> Number of grid layers (number of levels - (M+1) )
integer(kind=iintegers), intent(in) :: nvars !> Number of scalar fields
integer(kind=iintegers), intent(in) :: nvar !> Number of scalar field
integer(kind=iintegers), intent(in) :: itherm(1:M+2) !> Grid levels between constant-density layers
real(kind=ireals), intent(in) :: h !> Lake depth, m
real(kind=ireals), intent(in) :: dt_out !> Output interval, s
real(kind=ireals), intent(in) :: time !> Time, s
real(kind=ireals), intent(in) :: hx1ml(1:M+1), hy1ml(1:M+1) !> Deviations of thickness of
!! constant-density layers
real(kind=ireals), intent(in) :: Tw(1:M+1) !> Profile of horizontally-averaged scalar
type(gridspacing_type), intent(in) :: gsp !> Grid spacing group
type(gridsize_type), intent(in) :: gs !> Grid size group
type(bathym), intent(in) :: bathymwater(1:M+1) !> Bathymetry group
type(grarr1), intent(in) :: rtemp !> Coordinates of a location in a lake
character(len=*), intent(in) :: basename
logical, intent(in) :: firstcall_ != .true.
!Output variables
!real(kind=ireals), intent(out) :: Temp !> Temperature in the point set by coordinates rtemp
!Local variables
real(kind=ireals), allocatable :: Temp(:) ! Scalar in the point set by coordinates rtemp
real(kind=ireals), allocatable :: dep(:,:), Hlars(:)
real(kind=ireals) :: amplx, amply, amplx0, amply0, sumh, sumh1, dsumh, &
& larthick, dlarthick, alpha, Lx, Ly, h_, x, y, z
real(kind=ireals), save :: tcount = 0.
integer(kind=iintegers) :: i, j, k !Loop indices
integer(kind=iintegers) :: nlevs
integer(kind=iintegers), allocatable, save :: out_unit(:)
character :: coords_point*6
if (.not. allocated(out_unit)) then
allocate(out_unit(1:nvars))
out_unit = lake_misc_unit_min
endif
if (firstcall_) then
call CHECK_UNIT(lake_misc_unit_min,lake_misc_unit_max,out_unit(nvar))
write (coords_point,'(2i3)') gs%ix, gs%iy
open(out_unit(nvar),file=outpath(1:len_trim(outpath))//'time_series/'// &
& basename(1:len_trim(basename))//coords_point//'.dat', status='unknown')
endif
if (int(time/(dt_out*hour_sec)) > tcount) then
allocate (Temp(1:rtemp%numarr))
Temp(:) = missing_value
! Counting number of constant-density layers
nlevs = 0
do i = 2, M+2
if (itherm(i) < 0) then
exit
else
nlevs = nlevs + 1
endif
enddo
allocate (dep(1:nlevs,1:2))
allocate (Hlars(1:nlevs))
forall (i=1:nlevs) Hlars(i) = h*sum(gsp%ddz05(itherm(i):itherm(i+1)-1))
! Cycle by locations of output
do k = 1, rtemp%numarr
do i = 1, nlevs
amplx0 = 0.; amply0 = 0.
if (i < nlevs) then
do j = nlevs, i+1, -1
amplx0 = amplx0 - aki(j,i)*0.5*pi*hx1ml(j)
amply0 = amply0 - bki(j,i)*0.5*pi*hy1ml(j)
enddo
endif
amplx = amplx0 - 0.5*hx1ml(i)*pi
amply = amply0 - 0.5*hy1ml(i)*pi
if (i == nlevs) then
x = 0.
else
x = sum(Hlars(i+1:nlevs))
endif
if (i == 1) then
Lx = bathymwater(1)%Lx
Ly = bathymwater(1)%Ly
else
Lx = bathymwater(itherm(i))%Lx_half
Ly = bathymwater(itherm(i))%Ly_half
endif
dep(i,2) = x + amplx0*sin(pi*rtemp%arr1(1,k)/Lx) + amply0*sin(pi*rtemp%arr1(2,k)/Ly)
dep(i,2) = h - dep(i,2)
dep(i,1) = x + Hlars(i) + &
& amplx *sin(pi*rtemp%arr1(1,k)/Lx) + amply *sin(pi*rtemp%arr1(2,k)/Ly)
dep(i,1) = h - dep(i,1)
!if (i == 1) print*, 'i=1:', Hlars(i), 0.5*hy1ml(i)*pi, dep(i,1:2)
enddo
!if (k == 7 .or. k == 3) then
!print*, 'hx1ml=', hx1ml(1:10)
!do i = 1, nlevs
! print*, 'k=', k, 'dep1,dep2', i, dep(i,1:2)
!enddo
!endif
!!read*
! Computing temperature at a given depth as linear interpolation between model levels
seekl:do i = 1, nlevs
if ( i /= nlevs) then
if ( rtemp%arr1(3,k) < dep(i,2) .and. rtemp%arr1(3,k) > dep(i+1,1) ) then
z = abs(rtemp%arr1(3,k) - dep(i,2)) + abs(rtemp%arr1(3,k) - dep(i+1,1))
x = abs(rtemp%arr1(3,k) - dep(i,2))/z
y = abs(rtemp%arr1(3,k) - dep(i+1,1))/z
Temp(k) = y*Tw(itherm(i+1)) + x*Tw(itherm(i+1)+1)
!if (k == 7 .and. Temp(k) > 16.) print*, 'Temploc1: ', dep(i,2), dep(i+1,1)
exit seekl
endif
if ( rtemp%arr1(3,k) > dep(i,2) .and. rtemp%arr1(3,k) < dep(i+1,1) ) then
z = abs(rtemp%arr1(3,k) - dep(i,2)) + abs(rtemp%arr1(3,k) - dep(i+1,1))
x = abs(rtemp%arr1(3,k) - dep(i,2))/z
y = abs(rtemp%arr1(3,k) - dep(i+1,1))/z
Temp(k) = y*Tw(itherm(i+1)) + x*Tw(itherm(i+1)+1)
!if (k == 7 .and. Temp(k) > 16.) print*, 'Temploc2: ', dep(i,2), dep(i+1,1)
exit seekl
endif
endif
if (rtemp%arr1(3,k) > dep(i,1) .and. rtemp%arr1(3,k) <= dep(i,2)) then
y = h*(dep(i,2) - dep(i,1))/Hlars(i)
! Near top of the layer
if (i /= 1) then
if ( rtemp%arr1(3,k) < dep(i,1) + y*gsp%ddz(itherm(i))*0.5 ) then
if (dep(i-1,2)-dep(i-1,1) > 0.) then
x = 0.5*gsp%ddz(itherm(i))*h*(dep(i-1,2)-dep(i-1,1))/Hlars(i-1)
Temp(k) = Tw(itherm(i)) + (Tw(itherm(i)+1) - Tw(itherm(i))) * &
& (rtemp%arr1(3,k) - dep(i,1) + x)/(x + y*gsp%ddz(itherm(i))*0.5)
else
Temp(k) = 0.5*(Tw(itherm(i)) + Tw(itherm(i)+1))
endif
!if (k == 7 .and. Temp(k) > 16.) print*, 'Temploc3: ', dep(i,2), dep(i,1)
exit seekl
endif
endif
! Near bottom of the layer
if (i /= nlevs) then
if ( rtemp%arr1(3,k) > dep(i,2) - y*gsp%ddz(itherm(i+1))*0.5 ) then
if (dep(i+1,2)-dep(i+1,1) > 0.) then
x = 0.5*gsp%ddz(itherm(i+1))*h*(dep(i+1,2)-dep(i+1,1))/Hlars(i+1)
Temp(k) = Tw(itherm(i+1)) + (Tw(itherm(i+1)+1) - Tw(itherm(i+1))) * &
& (rtemp%arr1(3,k) - dep(i,2) + y*gsp%ddz(itherm(i+1))*0.5) / &
& (x + y*gsp%ddz(itherm(i+1))*0.5)
else
Temp(k) = 0.5*(Tw(itherm(i+1)) + Tw(itherm(i+1)+1))
endif
!if (k == 7 .and. Temp(k) > 16.) print*, 'Temploc5: ', x, y, dep(i,2), i
exit seekl
endif
endif
x = dep(i,1)
if (i /= 1) x = x + y*gsp%ddz(itherm(i))*0.5
j = itherm(i)+1
! Linear interpolation between numerical levels
seekl_1 : do
if (x < rtemp%arr1(3,k) .and. &
& x + y*gsp%ddz(j) >= rtemp%arr1(3,k)) then
Temp(k) = Tw(j) + (Tw(j+1) - Tw(j)) * &
& (rtemp%arr1(3,k) - x)/(y*gsp%ddz(j))
!if (k == 7 .and. Temp(k) > 16.) print*, 'Temploc4: ', x, y
exit seekl_1
endif
x = x + y*gsp%ddz(j)
j = j + 1
enddo seekl_1
endif
enddo seekl
!if (k == 7 .and. Temp(k) > 16.) then
! print*, 'In Temploc: STOP'
! read* !STOP
!endif
!sumh = 0.
!cyc1 : do i = M+1, 1, -1
! if1 : if (itherm(i+1) > 0) then
! Lx = bathymwater(itherm(2))%Lx!bathymwater(itherm(i+1))%Lx ! Assuming cross-section
! ! to be constant with depth
! Ly = bathymwater(itherm(2))%Ly!bathymwater(itherm(i+1))%Ly
! !if (abs(rtemp%arr1(1,k)) > 0.5*Lx .or. abs(rtemp%arr1(2,k)) > 0.5*Ly) exit if1
! amplx = - 0.5*hx1ml(i)*pi !/(Lx*Ly)
! amply = - 0.5*hy1ml(i)*pi !/(Lx*Ly)
! larthick = h*sum(gsp%ddz05(itherm(i):itherm(i+1)-1)) ! Equilibrium thickness of constant-density layer
! dlarthick = amplx*sin(pi*rtemp%arr1(1,k)/Lx) + amply*sin(pi*rtemp%arr1(2,k)/Ly)
! sumh = sumh + larthick + dlarthick
! !print*, 'ay', hy1ml
! endif if1
! if (sumh >= h - rtemp%arr1(3,k)) then
! !print*, itherm(i)+1, amplx
! h_ = h - rtemp%arr1(3,k)
! sumh1 = sumh
! cyclrs : do j = itherm(i)+1, itherm(i+1)
! dsumh = h*gsp%ddz05(j-1)*(1. + dlarthick/larthick)
! if ( h_ <= sumh1 .and. h_ > sumh1 - 0.5*dsumh ) then
! alpha = (h_ - sumh1 + 0.5*dsumh)/dsumh !Assuming regular grid and the same thickness deviation
! !in adjacent layers
! Temp(k) = Tw(max(j-1,1))*alpha + Tw(j)*(1. - alpha)
! !print*, '1', j, nvar
! exit cyclrs
! elseif ( h_ <= sumh1 - 0.5*dsumh .and. h_ > sumh1 - dsumh ) then
! alpha = (sumh1 - 0.5*dsumh - h_)/dsumh !Assuming regular grid and the same thickness deviation
! !in adjacent layers
! Temp(k) = Tw(min(j+1,M+1))*alpha + Tw(j)*(1. - alpha)
! !print*, '2', j, nvar
! exit cyclrs
! endif
! sumh1 = sumh1 - dsumh
! enddo cyclrs
! exit cyc1
! endif
!enddo cyc1
enddo
write(out_unit(nvar),*) time/hour_sec, Temp(1:rtemp%numarr)
if (nvar == nvars) tcount = tcount + 1.
!print*, Temp
!read*
deallocate(Temp)
deallocate(dep)
deallocate(Hlars)
endif
!if (firstcall) firstcall = .false.
END SUBROUTINE TEMPLOC
END MODULE OUT_MOD