Newer
Older
program CALCULATE_OBL
use EQUATION
use RICHARDSON
use PACANOWSKI_PHILANDER
use PACANOWSKI_PHILANDER_PLUS
use K_EPSILON
use io
use io_metadata
!use vertical_mixing
use vermix
implicit none
real :: t !Time
integer :: nt !Number of steps in t (time)
real :: z !Depth
integer :: nz !Number of steps in z (depth)
real :: dt !Timestep
real, allocatable :: dz(:) !z-step
integer :: i, k
integer :: status, num
real :: time_begin, time_end, time_current !time for output
real, allocatable :: Theta(:)
real, allocatable :: Salin(:)
real, allocatable :: U(:)
real, allocatable :: V(:)
real, allocatable :: Km(:)
real, allocatable :: Kh(:)
real, allocatable :: Flux_heat_surf(:), Flux_heat_bot(:)
real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:)
real, allocatable :: Tau_x_surf(:), Tau_y_surf(:), Tau_x_bot(:), Tau_y_bot(:)
real, allocatable :: flux_u_surf(:), flux_u_bot(:)
real, allocatable :: flux_v_surf(:), flux_v_bot(:)
real, allocatable :: f_Heat_right(:)
real, allocatable :: f_Sal_right(:)
real, allocatable :: f_u_right(:)
real, allocatable :: f_v_right(:)
real, allocatable :: Rho(:), N2(:), S2(:), Ri_grad(:)
real :: hml, lab_hml
integer :: maxcellnumber
real, allocatable :: TKE(:), eps(:), P_TKE(:), B_TKE(:)
real, allocatable :: Theta_write(:,:)
real, allocatable :: Salin_write(:,:)
real, allocatable :: U_write(:,:)
real, allocatable :: V_write(:,:)
real, allocatable :: Rho_write(:,:)
real, allocatable :: N2_write(:,:)
real, allocatable :: S2_write(:,:)
real, allocatable :: Ri_grad_write(:,:)
real, allocatable :: Kh_write(:,:)
real, allocatable :: Km_write(:,:)
real, allocatable :: hml_write(:)
real, allocatable :: maxcellnumber_write(:)
real, allocatable :: Theta_C_write(:,:)
real, allocatable :: lab_hml_write(:)
closure_mode = 1 ! 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
allocate(Theta(1:nz))
allocate(Salin(1:nz))
allocate(U(1:nz))
allocate(V(1:nz))
allocate(flux_Heat_surf(1:nt))
allocate(flux_Heat_bot(1:nt))
allocate(flux_Sal_surf(1:nt))
allocate(flux_Sal_bot(1:nt))
allocate(Tau_x_surf(1:nt))
allocate(Tau_y_surf(1:nt))
allocate(Tau_x_bot(1:nt))
allocate(Tau_y_bot(1:nt))
allocate(flux_u_surf(1:nt))
allocate(flux_u_bot(1:nt))
allocate(flux_v_surf(1:nt))
allocate(flux_v_bot(1:nt))
allocate(f_Heat_right(1:nz))
allocate(f_Sal_right(1:nz))
allocate(f_u_right(1:nz))
allocate(f_v_right(1:nz))
allocate(Km(1:nz))
allocate(Kh(1:nz))
allocate(Rho(1:nz))
allocate(N2(1:nz))
allocate(S2(1:nz))
allocate(Ri_grad(1:nz)) !
allocate(TKE(1:nz)) !
allocate(eps(1:nz)) !
allocate(P_TKE(1:nz)) !
allocate(B_TKE(1:nz)) !
allocate(Km_TKE(1:nz)) !
allocate(Km_eps(1:nz)) !
allocate(Theta_write(nz, nt))
allocate(Salin_write(nz, nt))
allocate(U_write(nz, nt))
allocate(V_write(nz, nt))
allocate(Rho_write(nz, nt))
allocate(N2_write(nz, nt))
allocate(S2_write(nz, nt))
allocate(Ri_grad_write(nz, nt))
allocate(Kh_write(nz, nt))
allocate(Km_write(nz, nt))
allocate(hml_write(nt))
allocate(maxcellnumber_write(nt))
allocate(Theta_C_write(nz, nt))
call Theta_init (Theta, nz, dz)
call Salin_init(Salin, nz, dz)
call set_Flux_heat_surf (Flux_heat_surf, nt)
call set_Flux_heat_bot (Flux_heat_bot, nt)
call set_Flux_sal_surf (Flux_sal_surf, nt)
call set_Tau_x_surf(Tau_x_surf, flux_u_surf, nt)
call set_Tau_y_surf(Tau_y_surf, flux_v_surf, nt)
call set_Tau_x_bot(Tau_x_bot, flux_u_bot, nt)
call set_Tau_y_bot(Tau_y_bot, flux_v_bot, nt)
call set_f_Heat_right (f_Heat_right, nz)
call set_f_Sal_right (f_Sal_right, nz)
call set_f_u_right (f_v_right, V, nz)
open (199, file= 'output_Daria/hml.txt', status ='replace')
open (20, file= 'output_Daria/surf_temp.txt', status ='replace')
status = 0
num = 0
time_begin = 0.0
time_end = t
time_current = time_begin
if (closure_mode.eq.3 .or. closure_mode.eq.4) then
call TKE_init(TKE, nz)
call eps_init(eps, nz, z)
endif
do k = 1, nz
Km(k) = 10.0**(-1.0)
Kh(k) = 10.0**(-1.0)
end do
if (closure_mode.eq.1) then
call solve_state_eq (Rho, Theta, Salin, nz)
call get_N2 (N2, Rho, dz, nz)
call get_S2 (S2, U, V, dz, nz)
call get_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
call get_lab_hml(lab_hml, flux_u_surf(i), flux_v_surf(i), dt*i, z)
!call get_Km(Km, Ri_grad, nz)
!call get_Kh(Kh, Ri_grad, nz)
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
call solve_scalar_eq (Theta, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq (Salin, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
else if (closure_mode.eq.2) then
call solve_state_eq (Rho, Theta, Salin, nz)
call get_N2 (N2, Rho, dz, nz)
call get_S2 (S2, U, V, dz, nz)
call get_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
call get_lab_hml(lab_hml, flux_u_surf(i), flux_v_surf(i), dt*i, z)
call get_Km_plus(Km, Ri_grad, flux_u_surf(i), flux_v_surf(i), hml, nz, dz, z)
call get_Kh_plus(Kh, Km, nz)
call get_TKE_generation(P_TKE, Km, S2, nz)
call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
call solve_scalar_eq (Theta, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq (Salin, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
else if (closure_mode.eq.3) then
call solve_state_eq (Rho, Theta, Salin, nz)
call get_N2 (N2, Rho, dz, nz)
call get_S2 (S2, U, V, dz, nz)
call get_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
call get_lab_hml(lab_hml, flux_u_surf(i), flux_v_surf(i), dt*i, z)
call TKE_bc(TKE, flux_u_surf(i), flux_v_surf(i), nz)
call eps_bc(eps, flux_u_surf(i), flux_v_surf(i), nz, dz)
call get_Km_k_eps (Km, TKE, eps, nz)
call get_Kh_k_eps (Kh, Km, nz)
call get_TKE_generation(P_TKE, Km, S2, nz)
call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
call TKE_bc(TKE, flux_u_surf(i), flux_v_surf(i), nz)
call eps_bc(eps, flux_u_surf(i), flux_v_surf(i), nz, dz)
call solve_TKE_eq (TKE, Km, nz, dz, dt, P_TKE, B_TKE, eps)
call solve_eps_eq (eps, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
call solve_scalar_eq (Theta, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq (Salin, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
else if (closure_mode.eq.4) then
call solve_state_eq (Rho, Theta, Salin, nz)
call get_N2 (N2, Rho, dz, nz)
call get_S2 (S2, U, V, dz, nz)
call get_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
call get_lab_hml(lab_hml, flux_u_surf(i), flux_v_surf(i), dt*i, z)
call TKE_bc(TKE, flux_u_surf(i), flux_v_surf(i), nz)
call eps_bc(eps, flux_u_surf(i), flux_v_surf(i), nz, dz)
call get_Km_k_eps (Km, TKE, eps, nz)
call get_Kh_k_eps (Kh, Km, nz)
call get_TKE_generation(P_TKE, Km, S2, nz)
call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
call TKE_bc(TKE, flux_u_surf(i), flux_v_surf(i), nz)
call eps_bc(eps, flux_u_surf(i), flux_v_surf(i), nz, dz)
call get_Km_TKE(Km_TKE, Km, nz)
call get_Km_eps(Km_eps, Km, nz)
call solve_TKE_eq_semiimplicit (TKE, Km_TKE, nz, dz, dt, P_TKE, B_TKE, eps)
call solve_eps_eq_semiimplicit (eps, Km_eps, nz, dz, dt, P_TKE, B_TKE, TKE)
call limit_min_array(eps, eps_min, nz)
call solve_scalar_eq (Theta, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq (Salin, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
else if (closure_mode.eq.5) then
call solve_state_eq(Rho, Theta, Salin, nz)
call vermix_init(rho, u, v, Tau_x_surf(i), Tau_y_surf(i), Ri_grad, Kh, Km)
call solve_scalar_eq(Theta, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq(Salin, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
call solve_vector_eq(U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
call solve_vector_eq(V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
!do k = 1, nz
! Theta (k) = Theta(k) + Theta_0
!end do
Theta_write(:, i) = Theta
Salin_write(:, i) = Salin
U_write(:, i) = U
V_write(:, i) = V
Rho_write(:, i) = Rho
N2_write(:, i) = N2
S2_write(:, i) = S2
Ri_grad_write(:, i) = Ri_grad
Kh_write(:, i) = Kh
Km_write(:, i) = Km
hml_write(i) = hml
maxcellnumber_write(i) = maxcellnumber
write(199,*) time_current, hml, lab_hml
write(20,*) time_current, Theta(nz)
i = i+1
time_current = time_current+dt
enddo
Theta_C_write = Theta_write - 273.15
! Writing 2D arrays (variables with depth and time)
call write_2d_real_nc(Theta_write, 'output.nc', meta_theta)
call write_2d_real_nc(Salin_write, 'output.nc', meta_salin)
call write_2d_real_nc(U_write, 'output.nc', meta_u)
call write_2d_real_nc(V_write, 'output.nc', meta_v)
call write_2d_real_nc(Rho_write, 'output.nc', meta_rho)
call write_2d_real_nc(N2_write, 'output.nc', meta_n2)
call write_2d_real_nc(S2_write, 'output.nc', meta_s2)
call write_2d_real_nc(Ri_grad_write, 'output.nc', meta_ri_grad)
call write_2d_real_nc(Kh_write, 'output.nc', meta_kh)
call write_2d_real_nc(Km_write, 'output.nc', meta_km)
call write_2d_real_nc(Theta_C_write, 'output.nc', meta_theta_C)
! Writing 1D arrays (scalar variables over time)
call write_1d_real_nc(hml_write, 'output.nc', meta_hml)
call write_1d_real_nc(maxcellnumber_write, 'output.nc', meta_maxcellnumber)
call write_1d_real_nc(Tau_x_surf, 'output.nc', meta_tau_u)
call write_1d_real_nc(Tau_y_surf, 'output.nc', meta_tau_v)
call write_1d_real_nc(lab_hml_write, 'output.nc', meta_lab_hml)
! Writing 2D arrays (variables with depth and time) to ASCII files
call write_2d_real_ascii(Theta_write, 'output_Theta.txt', meta_theta)
call write_2d_real_ascii(Salin_write, 'output_Salin.txt', meta_salin)
call write_2d_real_ascii(U_write, 'output_U.txt', meta_u)
call write_2d_real_ascii(V_write, 'output_V.txt', meta_v)
call write_2d_real_ascii(Rho_write, 'output_Rho.txt', meta_rho)
call write_2d_real_ascii(N2_write, 'output_N2.txt', meta_n2)
call write_2d_real_ascii(S2_write, 'output_S2.txt', meta_s2)
call write_2d_real_ascii(Ri_grad_write, 'output_Ri_grad.txt', meta_ri_grad)
call write_2d_real_ascii(Kh_write, 'output_Kh.txt', meta_kh)
call write_2d_real_ascii(Km_write, 'output_Km.txt', meta_km)
! Writing 1D arrays (scalar variables over time) to ASCII files
call write_1d_real_ascii(hml_write, 'output_hml.txt', meta_hml)
call write_1d_real_ascii(lab_hml_write, 'output_lab_hml.txt', meta_lab_hml)
call write_1d_real_ascii(maxcellnumber_write, 'output_maxcellnumber.txt', meta_maxcellnumber)
call write_1d_real_ascii(Tau_x_surf, 'output_tau_u.txt', meta_tau_u)
call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)
close(199)
close(20)
end program