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

latest versiom

parent f6d9ac18
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
calculate_obl.f90
+ 7
− 5
View file @ 9094fa4f
......@@ -73,7 +73,7 @@ program CALCULATE_OBL
real, allocatable :: Theta_C_write(:,:)
real, allocatable :: lab_hml_write(:)
closure_mode = 3 ! 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit
closure_mode = 1 ! 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit
call set_Time_Space(t, nt, z, nz)
......@@ -81,6 +81,8 @@ program CALCULATE_OBL
call set_dz_dt(dz, dt, nz, nt, t, z)
write (*,*) t, dt, nt, z, nz
allocate(Theta(1:nz))
allocate(Salin(1:nz))
allocate(U(1:nz))
......@@ -193,7 +195,7 @@ program CALCULATE_OBL
!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 get_N2 (N2, Rho, dz, nz)
......@@ -301,9 +303,9 @@ program CALCULATE_OBL
write(20,*) time_current, Theta(nz)
write(15,*) time_current
do k=1, nz
write(15,*) k, Theta(k), Salin(k), U(k), V(k), Kh(k), Km(k), TKE(k), eps(k)
enddo
!do k=1, nz
!write(15,*) k, Theta(k), Salin(k), U(k), V(k), Kh(k), Km(k), TKE(k), eps(k)
!enddo
i = i+1
time_current = time_current+dt
......
equation.f90
+ 53
− 2
View file @ 9094fa4f
......@@ -73,6 +73,31 @@ module EQUATION
f(nz) = - Field(nz) - f_right(nz) * dt - dt / dz(nz) * G(nz) * &
gamma + dt/dz(nz) * flux_surf * (1-gamma)
!Default case
!do i = 2, nz-1
! a(i) = (dt / dz(i) * Kvisc_minus(i)/dz_minus(i) * (1 - gamma)) / dt
! b(i) = (-1 - dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma) &
! - dt / dz(i) * Kvisc_minus(i)/dz_minus(i) * (1 - gamma)) /dt
! c(i) = (dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma)) / dt
! f(i) =( - Field(i) - f_right(i) * dt - dt / dz(i) * G(i) * gamma)/dt
!end do
!Bottom
! a(1) = 0/dt;
! b(1) = (-1 - dt / dz(1) * Kvisc_plus(1)/dz_plus(1) * (1 - gamma))/dt
! c(1) = (dt / dz(1) * Kvisc_plus(1)/dz_plus(1) * (1 - gamma))/dt
! f(1) = (- Field(1) - f_right(1) * dt - dt / dz(1) * G(1) * &
! gamma - dt/dz(1) * flux_bot * (1-gamma))/dt
!Коэффициенты на поверхности
! a(nz) = (dt / dz(nz) * Kvisc_minus(nz)/dz_minus(nz) * (1 - gamma))/dt
! b(nz) = (-1 - dt / dz(nz) * Kvisc_minus(nz)/dz_minus(nz) * (1 - gamma))/dt
! c(nz) = 0/dt
! f(nz) = (- Field(nz) - f_right(nz) * dt - dt / dz(nz) * G(nz) * &
! gamma + dt/dz(nz) * flux_surf * (1-gamma))/dt
!Call trigiagonal matrix solver
call tma (Field, a, b, c, f, nz)
......@@ -126,7 +151,7 @@ module EQUATION
!Default case
do i = 2, nz-1
a(i) = dt / dz(i) * Kvisc_minus(i)/dz_minus(i) * (1 - gamma)
b(i) = -1. - dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma) &
b(i) = -1 - dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma) &
- dt / dz(i) * Kvisc_minus(i)/dz_minus(i) * (1 - gamma)
c(i) = dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma)
f(i) = - Field(i) - f_right(i) * dt - dt / dz(i) * G(i) * gamma
......@@ -146,6 +171,32 @@ module EQUATION
f(nz) = - Field(nz) - f_right(nz) * dt - dt / dz(nz) * G(nz) * &
gamma + dt/dz(nz) * flux_surf * (1-gamma)
!Default case
! do i = 2, nz-1
! a(i) = (dt / dz(i) * Kvisc_minus(i)/dz_minus(i) * (1 - gamma)) / dt
! b(i) = (-1 - dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma) &
! - dt / dz(i) * Kvisc_minus(i)/dz_minus(i) * (1 - gamma)) /dt
! c(i) = (dt / dz(i) * Kvisc_plus(i)/dz_plus(i) * (1 - gamma)) / dt
! f(i) =( - Field(i) - f_right(i) * dt - dt / dz(i) * G(i) * gamma)/dt
! end do
!Bottom
! a(1) = 0/dt;
! b(1) = (-1 - dt / dz(1) * Kvisc_plus(1)/dz_plus(1) * (1 - gamma))/dt
! c(1) = (dt / dz(1) * Kvisc_plus(1)/dz_plus(1) * (1 - gamma))/dt
! f(1) = (- Field(1) - f_right(1) * dt - dt / dz(1) * G(1) * &
! gamma - dt/dz(1) * flux_bot * (1-gamma))/dt
!Коэффициенты на поверхности
! a(nz) = (dt / dz(nz) * Kvisc_minus(nz)/dz_minus(nz) * (1 - gamma))/dt
! b(nz) = (-1 - dt / dz(nz) * Kvisc_minus(nz)/dz_minus(nz) * (1 - gamma))/dt
! c(nz) = 0/dt
! f(nz) = (- Field(nz) - f_right(nz) * dt - dt / dz(nz) * G(nz) * &
! gamma + dt/dz(nz) * flux_surf * (1-gamma))/dt
!write(*,*) b(nz), c(nz)
!Call trigiagonal matrix solver
......
forcing_and_boundary.f90
+ 4
− 1
View file @ 9094fa4f
......@@ -44,6 +44,7 @@ module FORCING_AND_BOUNDARY
set_Flux = 1
open (555, file= 'Kato-Phillips/Flux_heat_surface.txt', status ='old') !'Flux_heat_surface_PAPA.txt'
!open (555, file= 'PAPA_06_2017/Flux_heat_surface_PAPA.txt', status ='old')
status = 0
num = 0
if (set_Flux.eq.1) then
......@@ -73,6 +74,7 @@ module FORCING_AND_BOUNDARY
set_Tau = 1
open (188, file= 'Kato-Phillips/Tau_x.txt', status ='old')
!open (188, file= 'PAPA_06_2017/Tau_x_PAPA.txt', status ='old')
status = 0
num = 0
if (set_Tau.eq.1) then
......@@ -112,6 +114,7 @@ module FORCING_AND_BOUNDARY
set_Tau = 1
open (177, file= 'Kato-Phillips/Tau_y.txt', status ='old')
!open (177, file= 'PAPA_06_2017/Tau_y_PAPA.txt', status ='old')
status = 0
num = 0
if (set_Tau.eq.1) then
......
initial_conditions.f90
+ 3
− 1
View file @ 9094fa4f
......@@ -25,7 +25,8 @@ module INITIAL_CONDITIONS
init_Theta = 2
open (222, file= 'Kato-Phillips/Theta_start.txt', status ='old') !'Theta_PAPA.txt'
open (222, file= 'Kato-Phillips/Theta_start.txt', status ='old') !'Theta_PAPA.txt' Flux_heat_surface_PAPA.txt
!open (222, file= 'PAPA_06_2017/Theta_PAPA.txt', status ='old')
status = 0
num = 0
if (init_Theta.eq.1) then
......@@ -60,6 +61,7 @@ module INITIAL_CONDITIONS
init_Salin = 1 ! 1 - const, 2 - function(z), 3 - file
open (777, file= 'Kato-Phillips/Salin_start.txt', status ='old') !Salin_PAPA.txt
!open (777, file= 'PAPA_06_2017/Salin_PAPA.txt', status ='old') !Salin_PAPA.txt
status = 0
num = 0
if (init_Salin.eq.1) then
......
k_epsilon.f90
+ 1
− 1
View file @ 9094fa4f
......@@ -7,7 +7,7 @@ module K_EPSILON
implicit none
real, parameter :: C_mu = 0.09
real, parameter :: PrT = 0.8 !0.8 !1.25
real, parameter :: PrT = 1.5 !0.8 !1.25
real, parameter :: eps_min = 0.0000001
real, parameter :: TKE_min = 0.0001
real, parameter :: kappa_k_eps = 0.4
......
pacanowski_philander.f90
+ 3
− 1
View file @ 9094fa4f
......@@ -36,4 +36,6 @@ module PACANOWSKI_PHILANDER
end do
end subroutine
!subroutine get_eps_pacanowski_philander
end module
\ No newline at end of file
pacanowski_philander_plus.f90
+ 1
− 1
View file @ 9094fa4f
......@@ -8,7 +8,7 @@ module PACANOWSKI_PHILANDER_PLUS
real, parameter :: alpha_turb = 5
real, parameter :: n_turb = 2
real, parameter :: Km_b = 10.**(-4.) ! 5 * 10**(-6) inmcm
real, parameter :: PrT = 0.8
real, parameter :: PrT = 1.3 !0.8
real, parameter :: kappa_pph_plus = 0.4
contains
......
richardson.f90
+ 0
− 5
View file @ 9094fa4f
......@@ -7,11 +7,6 @@ module RICHARDSON
real, parameter :: Ri_grad_min = 0.0000001
real, parameter :: Ri_grad_max = 10
!real, parameter :: eps_min = 0.0000001
!real, parameter :: eps_max = 1.0
!real, parameter :: TKE_min = 0.0001
!real, parameter :: TKE_max = 1.0
!real, parameter :: kappa_k_eps = 0.4
contains
......
time_and_space.f90
+ 5
− 1
View file @ 9094fa4f
......@@ -14,6 +14,7 @@ module TIME_AND_SPACE
integer :: i, status, num
!open (1, file= 'PAPA_06_2017/time_space_for_calc_PAPA.txt', status ='old')
open (1, file= 'Kato-Phillips/time_space_for_calc.txt', status ='old')
status = 0
do while (status.eq.0)
......@@ -33,11 +34,13 @@ module TIME_AND_SPACE
integer :: init_dz, i, status, num
real :: dz_start
!dt = 5.0;
dt = t/nt
init_dz = 1 ! 1 - const, 2 - function(z), 3 - file
open (33, file= 'Kato-Phillips/dz.txt', status ='old') !'dz_PAPA.txt'
!open (33, file= 'PAPA_06_2017/dz_PAPA.txt', status ='old') !'dz_PAPA.txt'
status = 0
num = 0
if (init_dz.eq.1) then
......@@ -64,5 +67,6 @@ module TIME_AND_SPACE
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