Skip to content
Snippets Groups Projects
Commit b19e6e91 authored by Daria Gladskikh's avatar Daria Gladskikh :ocean:
Browse files

Changes in code organization

parent 5955562e
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
calculate_obl.f90
+ 165
− 310
View file @ b19e6e91
program CALCULATE_OBL program CALCULATE_OBL
use TIME_AND_SPACE
use INITIAL_CONDITIONS use INITIAL_CONDITIONS
use FORCING_AND_BOUNDARY use FORCING_AND_BOUNDARY
use RHS
use EQUATION use EQUATION
use RICHARDSON use RICHARDSON
use PACANOWSKI_PHILANDER use PACANOWSKI_PHILANDER
use PACANOWSKI_PHILANDER_PLUS use PACANOWSKI_PHILANDER_PLUS
use K_EPSILON use K_EPSILON
!use INITIAL_CONDITIONS
use io use io
use io_metadata use io_metadata
implicit none implicit none
integer :: nz !Number of steps in z (depth) real :: t !Time
integer :: nt !Number of steps in 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 :: dt !Timestep
real, allocatable :: dz(:) !z-step real, allocatable :: dz(:) !z-step
real :: t !Time
real :: z !Depth
integer :: i, k integer :: i, k
integer :: status, num integer :: status, num
real :: time_begin, time_end, time_current !time for output real :: time_begin, time_end, time_current !time for output
real :: dz_start, dz_last
!real :: Theta_start
!real :: Salin_start
!real :: Kh_start
!real :: U_start
!real :: V_start
!real :: Km_start
integer :: init_dz
!integer :: init_theta, init_sal, init_kh
!integer :: bound_flux_heat_bot, bound_flux_heat_surf, bound_f_heat_right
real :: flux_heat_surf_tmp
real :: flux_heat_bot_tmp
!integer :: init_u, init_v, init_km
integer :: bound_flux_u_bot
integer :: bound_flux_v_bot
integer :: init_tau_x, init_tau_y
integer :: bound_f_u_right, bound_f_v_right
real :: flux_u_bot_tmp
real :: flux_v_bot_tmp
real :: tau_x_tmp, tau_y_tmp
real, allocatable :: Theta(:) real, allocatable :: Theta(:)
real, allocatable :: f_heat_right(:)
real, allocatable :: Flux_heat_surf(:), flux_heat_bot(:)
real, allocatable :: Salin(:) real, allocatable :: Salin(:)
real, allocatable :: f_sal_right(:)
real, allocatable :: flux_sal_surf(:), flux_sal_bot(:)
real, allocatable :: Kh(:)
real, allocatable :: U(:) real, allocatable :: U(:)
real, allocatable :: V(:) real, allocatable :: V(:)
real, allocatable :: Km(:)
real, allocatable :: Kh(:)
real, allocatable :: f_u_right(:)
real, allocatable :: f_v_right(:) 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_u_surf(:), flux_u_bot(:)
real, allocatable :: flux_v_surf(:), flux_v_bot(:) real, allocatable :: flux_v_surf(:), flux_v_bot(:)
real, allocatable :: Km(:) real, allocatable :: f_Heat_right(:)
real, allocatable :: f_Sal_right(:)
real, allocatable :: f_u_right(:)
real, allocatable :: f_v_right(:)
real :: f_cor, Ug, Vg integer :: closure_mode
real, allocatable :: Rho(:), N2(:), S2(:), Ri_grad(:) real, allocatable :: Rho(:), N2(:), S2(:), Ri_grad(:)
real, allocatable :: tau(:), tau_x(:), tau_y(:)
!real, allocatable :: Rho_write(:,:), N2_write(:,:), S2_write(:,:), Ri_grad_write(:,:), Kh_write(:,:), Km_write(:,:)
real :: hml, lab_hml real :: hml, lab_hml
integer :: maxcellnumber integer :: maxcellnumber
real :: ustar
real, allocatable :: TKE(:), eps(:), P_TKE(:), B_TKE(:) real, allocatable :: TKE(:), eps(:), P_TKE(:), B_TKE(:)
real, allocatable :: Km_TKE(:), Km_eps(:) real, allocatable :: Km_TKE(:), Km_eps(:)
...@@ -109,49 +73,52 @@ program CALCULATE_OBL ...@@ -109,49 +73,52 @@ program CALCULATE_OBL
real, allocatable :: Theta_C_write(:,:) real, allocatable :: Theta_C_write(:,:)
real, allocatable :: lab_hml_write(:) real, allocatable :: lab_hml_write(:)
closure_mode = 3 ! 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 3 - k-epsilon semiimplicit
call set_Time_Space(t, nt, z, nz)
allocate(dz(1:nz))
!call set_time_space (t, nt, z, nz, dt, dz) call set_dz_dt(dz, dt, nz, nt, t, z)
!write (*,*) t, nt, z, nz, dt, dz
!Time and space
open (1, file= 'time_space_for_calc.txt', status ='old')
status = 0
do while (status.eq.0)
read (1, *, iostat = status) t, nt, z, nz
end do
close (1)
allocate(dz(1:nz))
allocate(Theta(1:nz)) allocate(Theta(1:nz))
allocate(Kh(1:nz))
allocate(f_heat_right(1:nz))
allocate(flux_heat_surf(1:nt))
allocate(flux_heat_bot(1:nt))
allocate(Salin(1:nz)) allocate(Salin(1:nz))
allocate(f_sal_right(1:nz))
allocate(flux_sal_surf(1:nt))
allocate(flux_sal_bot(1:nt))
allocate(U(1:nz)) allocate(U(1:nz))
allocate(V(1:nz)) allocate(V(1:nz))
allocate(Km(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_surf(1:nt))
allocate(flux_u_bot(1:nt)) allocate(flux_u_bot(1:nt))
allocate(flux_v_surf(1:nt)) allocate(flux_v_surf(1:nt))
allocate(flux_v_bot(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_u_right(1:nz))
allocate(f_v_right(1:nz)) allocate(f_v_right(1:nz))
allocate(Km(1:nz))
allocate(Kh(1:nz))
allocate(Rho(1:nz)) allocate(Rho(1:nz))
allocate(N2(1:nz)) allocate(N2(1:nz))
allocate(S2(1:nz)) allocate(S2(1:nz))
allocate(Ri_grad(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(tau(1:nt)) allocate(Km_TKE(1:nz)) !
allocate(tau_x(1:nt)) allocate(Km_eps(1:nz)) !
allocate(tau_y(1:nt))
allocate(Theta_write(nz, nt)) allocate(Theta_write(nz, nt))
allocate(Salin_write(nz, nt)) allocate(Salin_write(nz, nt))
...@@ -168,67 +135,6 @@ program CALCULATE_OBL ...@@ -168,67 +135,6 @@ program CALCULATE_OBL
allocate(Theta_C_write(nz, nt)) allocate(Theta_C_write(nz, nt))
allocate(lab_hml_write(nt)) allocate(lab_hml_write(nt))
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))
dt = t / nt
!Как будем задавать начальные распределения
init_dz = 1 ! 1 - const, 2 - function(z), 3 - file
!init_theta = 2 ! 1 - const, 2 - function(z), 3 - file
!init_kh = 1 ! 1 - const, 2 - function(z), 3 - file
!bound_flux_heat_bot = 1 ! 1 - const, 2 - function(z), 3 - file
!bound_flux_heat_surf = 3 ! 1 - const, 2 - function(z), 3 - file
!bound_f_heat_right = 1 ! 1 - const, 2 - function(z), 3 - file
!init_sal = 1 ! 1 - const, 2 - function(z), 3 - file
!init_u = 1 ! 1 - const, 2 - function(z), 3 - file
!init_v = 1 ! 1 - const, 2 - function(z), 3 - file
!init_km = 1 ! 1 - const, 2 - function(z), 3 - file
bound_flux_u_bot = 1! 1 - const, 2 - function(z), 3 - file
bound_flux_v_bot = 1 ! 1 - const, 2 - function(z), 3 - file
bound_f_u_right = 1 ! 1 - const, 2 - function(z), 3 - file
bound_f_v_right = 1 ! 1 - const, 2 - function(z), 3 - file
init_tau_x = 1 ! 1 - const, 2 - function(z), 3 - file
init_tau_y = 1 ! 1 - const, 2 - function(z), 3 - file
open (33, file= 'dz.txt', status ='old') !'dz_PAPA.txt'
status = 0
num = 0
if (init_dz.eq.1) then
do i = 1, nz
dz(i) = z / nz
end do
else if (init_dz.eq.2) then
read (33, *, iostat = status) dz_start, dz_last !just example, not for use
dz(1) = dz_start
dz(nz) = dz_last
do i = 2, nz-1
dz(i) = dz(1) + 0.5 * (dz(nz)-dz(1)) * i
end do
else if (init_dz.eq.3) then
do while (status.eq.0)
read (33, *, iostat = status) dz_start
num = num + 1
dz(num) = dz_start
if (num >= nz) then
exit
endif
end do
endif
close (33)
call Theta_init (Theta, nz, dz) call Theta_init (Theta, nz, dz)
...@@ -238,120 +144,29 @@ program CALCULATE_OBL ...@@ -238,120 +144,29 @@ program CALCULATE_OBL
call V_init(V, nz) call V_init(V, nz)
!do i=1, nz
! write (*,*) Theta(i), Salin(i), U(i), V(i)
!enddo
call set_Flux_heat_bot (Flux_heat_bot, nt)
call set_Flux_heat_surf (Flux_heat_surf, nt) call set_Flux_heat_surf (Flux_heat_surf, nt)
call set_Tau_x_surf(Tau_x, nt) call set_Flux_heat_bot (Flux_heat_bot, nt)
call set_Tau_y_surf(Tau_y, nt)
do i=1, nt
write (*,*) Flux_heat_bot(i), Flux_heat_surf(i), Tau_x(i), Tau_y(i)
enddo
open (11, file= 'Flux_u_bottom.txt', status ='old')
status = 0
num = 0
if (bound_flux_u_bot.eq.1) then
read (11, *, iostat = status) flux_u_bot_tmp
do i = 1, nt
flux_u_bot(i) = 0. !flux_u_bot_tmp;
end do
else if (bound_flux_u_bot.eq.2) then
read (11, *, iostat = status) flux_u_bot_tmp
do i = 1, nt
flux_u_bot(i) = flux_u_bot_tmp + 0 * z
end do
else if (bound_flux_u_bot.eq.3) then
do while (status.eq.0)
read (11, *, iostat = status) flux_u_bot_tmp
num = num + 1
flux_u_bot(num) = flux_u_bot_tmp
if (num >= nt) then
exit
endif
end do
endif
close (11)
open (13, file= 'Flux_v_bottom.txt', status ='old')
status = 0
num = 0
if (bound_flux_v_bot.eq.1) then
read (13, *, iostat = status) flux_v_bot_tmp
do i = 1, nt
flux_v_bot(i) = 0 !flux_v_bot_tmp;
end do
else if (bound_flux_v_bot.eq.2) then
read (13, *, iostat = status) flux_v_bot_tmp
do i = 1, nt
flux_v_bot(i) = flux_v_bot_tmp + 0 * z
end do
else if (bound_flux_v_bot.eq.3) then
do while (status.eq.0)
read (13, *, iostat = status) flux_v_bot_tmp
num = num + 1
flux_v_bot(num) = flux_v_bot_tmp
if (num >= nt) then
exit
endif
end do
endif
close (13)
call set_Flux_sal_surf (Flux_sal_surf, nt)
do i = 1, nt call set_Flux_sal_bot (Flux_sal_bot, nt)
flux_u_surf(i) = tau_x(i) / Rho_0
flux_v_surf(i) = tau_y(i) / Rho_0
end do
call set_Tau_x_surf(Tau_x_surf, flux_u_surf, nt)
f_cor = 0 call set_Tau_y_surf(Tau_y_surf, flux_v_surf, nt)
Ug = 0
Vg = 0
do i = 1, nz call set_Tau_x_bot(Tau_x_bot, flux_u_bot, nt)
f_u_right(i) = - f_cor * (V(i) - Vg)
end do
do i = 1, nz call set_Tau_y_bot(Tau_y_bot, flux_v_bot, nt)
f_v_right(i) = + f_cor * (U(i) - Ug)
end do
!Правая часть call set_f_Heat_right (f_Heat_right, nz)
do i = 1, nz
f_heat_right(i) = 0
end do
do i = 1, nz
f_sal_right(i) = 0
end do
do i = 1, nt
flux_sal_surf(i) = 0
end do
do i = 1, nt call set_f_Sal_right (f_Sal_right, nz)
flux_sal_bot(i) = 0
end do
do i = 2, nz-1 call set_f_u_right (f_v_right, V, nz)
TKE(i) = 1 / (C_mu)**(1.0/2.0) * 0.001 * 0.001
end do
do i = 2, nz-1 call set_f_v_right (f_v_right, U, nz)
eps(i) = 0.001 * 0.001 * 0.001 / (kappa_k_eps * z)
end do
open (199, file= 'hml.txt', status ='replace') open (199, file= 'hml.txt', status ='replace')
...@@ -366,6 +181,12 @@ program CALCULATE_OBL ...@@ -366,6 +181,12 @@ program CALCULATE_OBL
time_end = t time_end = t
time_current = time_begin 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
i=1 i=1
!time_count = 0 !time_count = 0
!time_mark = 10 !time_mark = 10
...@@ -373,43 +194,77 @@ program CALCULATE_OBL ...@@ -373,43 +194,77 @@ program CALCULATE_OBL
!write(*,*) "Hello" !write(*,*) "Hello"
do while (time_current <= time_end .and. i<=nt) do while (time_current <= time_end .and. i<=nt)
if (closure_mode.eq.1) then
call solve_state_eq (Rho, Theta, Salin, nz) call solve_state_eq (Rho, Theta, Salin, nz)
call get_N2 (N2, Rho, dz, nz) call get_N2 (N2, Rho, dz, nz)
call get_S2 (S2, U, V, 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_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_lab_hml(lab_hml, flux_u_surf(i), flux_v_surf(i), dt*i, z)
!write (*,*) i, lab_hml
!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_Km(Km, Ri_grad, nz) call get_Km(Km, Ri_grad, nz)
call get_Kh(Kh, Ri_grad, nz) call get_Kh(Kh, Ri_grad, nz)
!call get_TKE_generation(P_TKE, Km, S2, nz) call get_TKE_generation(P_TKE, Km, S2, nz)
!call get_TKE_buoyancy(B_TKE, Kh, N2, 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.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_Rig(Ri_grad, N2, S2, nz)
call get_hml (hml, maxcellnumber, N2, dz, nz, z)
!call get_Km_TKE(Km_TKE, Km, nz) call get_lab_hml(lab_hml, flux_u_surf(i), flux_v_surf(i), dt*i, z)
!call get_Km_eps(Km_eps, Km, nz) call get_Km_plus(Km, Ri_grad, flux_u_surf(i), flux_v_surf(i), hml, nz, dz, z)
!call solve_TKE_eq (TKE, Km, nz, dz, dt, P_TKE, B_TKE, eps) call get_Kh_plus(Kh, Km, nz)
!call solve_TKE_eq_semiimplicit (TKE, Km_TKE, nz, dz, dt, P_TKE, B_TKE, eps) call get_TKE_generation(P_TKE, Km, S2, nz)
!call solve_eps_eq_semiimplicit (eps, Km_eps, nz, dz, dt, P_TKE, B_TKE, TKE) call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
!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 get_Km_plus(Km, Ri_grad, flux_u_surf(i), flux_v_surf(i), hml, nz, dz, z) call solve_scalar_eq (Salin, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right)
!call get_Kh_plus(Kh, Km, nz) call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right)
!call get_Km(Km, Ri_grad, nz) call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
!call get_Kh(Kh, Ri_grad, nz) 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 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 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 solve_scalar_eq (Theta, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right) 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_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 (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) call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf(i), flux_v_bot(i), f_v_right)
endif
Theta_write(:, i) = Theta Theta_write(:, i) = Theta
Salin_write(:, i) = Salin Salin_write(:, i) = Salin
...@@ -456,8 +311,8 @@ program CALCULATE_OBL ...@@ -456,8 +311,8 @@ program CALCULATE_OBL
! Writing 1D arrays (scalar variables over time) ! Writing 1D arrays (scalar variables over time)
call write_1d_real_nc(hml_write, 'output.nc', meta_hml) 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(maxcellnumber_write, 'output.nc', meta_maxcellnumber)
call write_1d_real_nc(tau_x, 'output.nc', meta_tau_u) call write_1d_real_nc(Tau_x_surf, 'output.nc', meta_tau_u)
call write_1d_real_nc(tau_y, 'output.nc', meta_tau_v) 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) call write_1d_real_nc(lab_hml_write, 'output.nc', meta_lab_hml)
! Writing 2D arrays (variables with depth and time) to ASCII files ! Writing 2D arrays (variables with depth and time) to ASCII files
...@@ -476,8 +331,8 @@ program CALCULATE_OBL ...@@ -476,8 +331,8 @@ program CALCULATE_OBL
call write_1d_real_ascii(hml_write, 'output_hml.txt', meta_hml) 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(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(maxcellnumber_write, 'output_maxcellnumber.txt', meta_maxcellnumber)
call write_1d_real_ascii(tau_x, 'output_tau_u.txt', meta_tau_u) call write_1d_real_ascii(Tau_x_surf, 'output_tau_u.txt', meta_tau_u)
call write_1d_real_ascii(tau_y, 'output_tau_v.txt', meta_tau_v) call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)
close(15) close(15)
close(199) close(199)
......
forcing_and_boundary.f90
+ 114
− 4
View file @ b19e6e91
...@@ -63,10 +63,10 @@ module FORCING_AND_BOUNDARY ...@@ -63,10 +63,10 @@ module FORCING_AND_BOUNDARY
close (555) close (555)
end subroutine end subroutine
subroutine set_Tau_x_surf(Tau_x_surf, nt) subroutine set_Tau_x_surf(Tau_x_surf, flux_u_surf, nt)
integer, intent(in) :: nt integer, intent(in) :: nt
real, intent(out) :: Tau_x_surf(nt) real, intent(out) :: Tau_x_surf(nt), flux_u_surf(nt)
real :: Tau_x_surf_tmp real :: Tau_x_surf_tmp
integer :: i, status, num, set_Tau integer :: i, status, num, set_Tau
...@@ -95,12 +95,17 @@ module FORCING_AND_BOUNDARY ...@@ -95,12 +95,17 @@ module FORCING_AND_BOUNDARY
end do end do
endif endif
close (188) close (188)
do i = 1, nt
flux_u_surf(i) = Tau_x_surf(i) / Rho_0
end do
end subroutine end subroutine
subroutine set_Tau_y_surf(Tau_y_surf, nt) subroutine set_Tau_y_surf(Tau_y_surf, flux_v_surf, nt)
integer, intent(in) :: nt integer, intent(in) :: nt
real, intent(out) :: Tau_y_surf(nt) real, intent(out) :: Tau_y_surf(nt), flux_v_surf(nt)
real :: Tau_y_surf_tmp real :: Tau_y_surf_tmp
integer :: i, status, num, set_Tau integer :: i, status, num, set_Tau
...@@ -129,7 +134,112 @@ module FORCING_AND_BOUNDARY ...@@ -129,7 +134,112 @@ module FORCING_AND_BOUNDARY
end do end do
endif endif
close (177) close (177)
do i = 1, nt
flux_v_surf(i) = Tau_y_surf(i) / Rho_0
end do
end subroutine
subroutine set_Tau_x_bot(Tau_x_bot, flux_u_bot, nt)
integer, intent(in) :: nt
real, intent(out) :: Tau_x_bot(nt), flux_u_bot(nt)
real :: Tau_x_bot_tmp
integer :: i, status, num, set_Tau
set_Tau = 1
!open (188, file= 'Tau_x.txt', status ='old')
!status = 0
!num = 0
!if (set_Tau.eq.1) then
do i = 1, nt
Tau_x_bot(i) = 0.0 !flux_v_surf_tmp;
end do
!else if (set_Tau.eq.2) then
! read (188, *, iostat = status) Tau_x_surf_tmp
! do i = 1, nt
! Tau_x_surf(i) = Tau_x_surf_tmp
! end do
!else if (set_Tau.eq.3) then
! do while (status.eq.0)
! read (188, *, iostat = status) Tau_x_surf_tmp
! num = num + 1
! Tau_x_surf(num) = Tau_x_surf_tmp
! if (num >= nt) then
! exit
! endif
! end do
!endif
!close (188)
do i = 1, nt
flux_u_bot(i) = 0.0 !Tau_x_bot(i) / Rho_0
end do
end subroutine
subroutine set_Tau_y_bot(Tau_y_bot, flux_v_bot, nt)
integer, intent(in) :: nt
real, intent(out) :: Tau_y_bot(nt), flux_v_bot(nt)
real :: Tau_y_bot_tmp
integer :: i, status, num, set_Tau
set_Tau = 1
!open (177, file= 'Tau_y.txt', status ='old')
!status = 0
!num = 0
!if (set_Tau.eq.1) then
do i = 1, nt
Tau_y_bot(i) = 0.0
end do
!else if (set_Tau.eq.2) then
! read (177, *, iostat = status) Tau_y_surf_tmp
! do i = 1, nt
! Tau_y_surf(i) = Tau_y_surf_tmp
! end do
!else if (set_Tau.eq.3) then
! do while (status.eq.0)
! read (177, *, iostat = status) Tau_y_surf_tmp
! num = num + 1
! Tau_y_surf(num) = Tau_y_surf_tmp
! if (num >= nt) then
! exit
! endif
! end do
!endif
!close (177)
do i = 1, nt
flux_v_bot(i) = 0.0 !Tau_y_bot(i) / Rho_0
end do
end subroutine end subroutine
subroutine set_Flux_sal_surf(Flux_sal_surf, nt)
integer, intent(in) :: nt
real, intent(out) :: Flux_sal_surf(nt)
integer :: i
do i = 1, nt
Flux_sal_surf(i) = 0.0
end do
end subroutine
subroutine set_Flux_sal_bot(Flux_sal_bot, nt)
integer, intent(in) :: nt
real, intent(out) :: Flux_sal_bot(nt)
integer :: i
do i = 1, nt
Flux_sal_bot(i) = 0.0
end do
end subroutine
end module end module
k_epsilon.f90
+ 22
− 0
View file @ b19e6e91
...@@ -18,6 +18,28 @@ module K_EPSILON ...@@ -18,6 +18,28 @@ module K_EPSILON
real, parameter :: Sch_eps = 1.3 !1.11 real, parameter :: Sch_eps = 1.3 !1.11
contains contains
subroutine TKE_init (TKE, nz)
integer, intent(in) :: nz
real, intent(out) :: TKE(nz)
integer :: i
do i = 2, nz-1
TKE(i) = 1 / (C_mu)**(1.0/2.0) * 0.001 * 0.001
end do
end subroutine
subroutine eps_init (eps, nz, z)
integer, intent(in) :: nz
real, intent(in) :: z
real, intent(out) :: eps(nz)
integer :: i
do i = 2, nz-1
eps(i) = 0.001 * 0.001 * 0.001 / (kappa_k_eps * z)
end do
end subroutine
subroutine get_TKE_generation (P_TKE, Km, S2, nz) subroutine get_TKE_generation (P_TKE, Km, S2, nz)
real, intent(in) :: S2(nz) real, intent(in) :: S2(nz)
real, intent(in) :: Km(nz) real, intent(in) :: Km(nz)
......
rhs.f90 0 → 100644
+ 54
− 0
View file @ b19e6e91
module RHS
implicit none
real, parameter :: f_cor = 0
real, parameter :: Ug = 0
real, parameter :: Vg = 0
contains
subroutine set_f_Heat_right (f_Heat_right, nz)
integer, intent(in) :: nz
real, intent(out) :: f_Heat_right(nz)
integer :: i
do i = 1, nz
f_Heat_right(i) = 0
end do
end subroutine
subroutine set_f_Sal_right (f_Sal_right, nz)
integer, intent(in) :: nz
real, intent(out) :: f_Sal_right(nz)
integer :: i
do i = 1, nz
f_Sal_right(i) = 0
end do
end subroutine
subroutine set_f_u_right (f_u_right, V, nz)
integer, intent(in) :: nz
real, intent(in) :: V(nz)
real, intent(out) :: f_u_right(nz)
integer :: i
do i = 1, nz
f_u_right(i) = - f_cor * (V(i) - Vg)
end do
end subroutine
subroutine set_f_v_right (f_v_right, U, nz)
integer, intent(in) :: nz
real, intent(in) :: U(nz)
real, intent(out) :: f_v_right(nz)
integer :: i
do i = 1, nz
f_v_right(i) = - f_cor * (U(i) - Vg)
end do
end subroutine
end module
\ No newline at end of file
time_and_space.f90 0 → 100644
+ 68
− 0
View file @ b19e6e91
module TIME_AND_SPACE
implicit none
contains
subroutine set_Time_Space(t, nt, z, nz)
real, intent(out) :: t !Time
integer, intent(out) :: nt !Number of steps in t (time)
real, intent(out) :: z !Depth
integer, intent(out) :: nz !Number of steps in z (depth)
integer :: i, status, num
open (1, file= 'time_space_for_calc.txt', status ='old')
status = 0
do while (status.eq.0)
read (1, *, iostat = status) t, nt, z, nz
end do
close (1)
end subroutine
subroutine set_dz_dt(dz, dt, nz, nt, t, z)
integer, intent(in) :: nz, nt
real, intent(in) :: t, z
real, intent(out) :: dz(nz) !z-step
real, intent(out) :: dt
integer :: init_dz, i, status, num
real :: dz_start
dt = t / nt
init_dz = 1 ! 1 - const, 2 - function(z), 3 - file
open (33, file= 'dz.txt', status ='old') !'dz_PAPA.txt'
status = 0
num = 0
if (init_dz.eq.1) then
do i = 1, nz
dz(i) = z / nz
end do
!else if (init_dz.eq.2) then
!read (33, *, iostat = status) dz_start, dz_last !just example, not for use
!dz(1) = dz_start
!dz(nz) = dz_last
!do i = 2, nz-1
!dz(i) = dz(1) + 0.5 * (dz(nz)-dz(1)) * i
!end do
else if (init_dz.eq.3) then
do while (status.eq.0)
read (33, *, iostat = status) dz_start
num = num + 1
dz(num) = dz_start
if (num >= nz) then
exit
endif
end do
endif
close (33)
end subroutine
end module
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment