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

interpolation added

parent a56d26d6
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Kato-Phillips/time_space_for_forcing.txt 0 → 100644
+ 1
− 0
View file @ 0d97260d
10 1
PAPA_06_2017/time_space_for_calc_PAPA.txt
+ 1
− 1
View file @ 0d97260d
1551600 431 200 14 1555200 1555200 200 14
PAPA_06_2017/time_space_for_forcing_PAPA.txt 0 → 100644
+ 1
− 0
View file @ 0d97260d
432 3600
obl_forcing_and_boundary.f90
+ 81
− 66
View file @ 0d97260d
...@@ -11,11 +11,13 @@ module obl_forcing_and_boundary ...@@ -11,11 +11,13 @@ module obl_forcing_and_boundary
contains contains
subroutine set_Flux_heat_bot (Flux_heat_bot, nt) subroutine set_Flux_heat_bot (Flux_heat_bot, Time_flux_heat_bot, nf, df)
!< @brief set bottom heat flux !< @brief set bottom heat flux
integer, intent(in) :: nt !< number of timesteps integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: Flux_heat_bot(nt) !< bottom heat flux, input: [W/m**2], output: [K*m/s] real, intent(in) :: df
real, intent(out) :: Flux_heat_bot(nf) !< bottom heat flux, input: [W/m**2], output: [K*m/s]
real, intent(out) :: Time_flux_heat_bot(nf) !< Time array
real :: Flux_heat_bot_tmp !< temporary bottom heat flux real :: Flux_heat_bot_tmp !< temporary bottom heat flux
integer :: i !< counter integer :: i !< counter
integer :: status, num !< for file input integer :: status, num !< for file input
...@@ -27,15 +29,17 @@ module obl_forcing_and_boundary ...@@ -27,15 +29,17 @@ module obl_forcing_and_boundary
status = 0 status = 0
num = 0 num = 0
if (set_Flux.eq.1) then if (set_Flux.eq.1) then
do i = 1, nt do i = 1, nf
Flux_heat_bot(i) = 0 Flux_heat_bot(i) = 0
Time_flux_heat_bot(i) = i * df
end do end do
else if (set_Flux.eq.3) then else if (set_Flux.eq.3) then
do while (status.eq.0) do while (status.eq.0)
read (444, *, iostat = status) Flux_heat_bot_tmp read (444, *, iostat = status) Flux_heat_bot_tmp
num = num + 1 num = num + 1
Flux_heat_bot(num) = Flux_heat_bot_tmp/Rho_ref/cp_water Flux_heat_bot(num) = Flux_heat_bot_tmp/Rho_ref/cp_water
if (num >= nt) then Time_flux_heat_bot(num) = num * df
if (num >= nf) then
exit exit
endif endif
end do end do
...@@ -43,25 +47,18 @@ module obl_forcing_and_boundary ...@@ -43,25 +47,18 @@ module obl_forcing_and_boundary
close (444) close (444)
end subroutine end subroutine
subroutine set_Flux_heat_surf (Flux_heat_surf, Flux_heat_surf_time_forcing, nt) subroutine set_Flux_heat_surf (Flux_heat_surf, Time_flux_heat_surf, nf, df)
!< @brief set surface heat flux !< @brief set surface heat flux
integer, intent(in) :: nf !< number of timesteps
integer, intent(in) :: nt !< number of timesteps real, intent(in) :: df
real, intent(out) :: Flux_heat_surf(nt) !< surface heat flux, output: [K*m/s] real, intent(out) :: Flux_heat_surf(nf) !< surface heat flux, output: [K*m/s]
real, intent(out) :: Time_flux_heat_surf(nf) !< Time array
real :: Flux_heat_surf_tmp !< temporary surface heat flux real :: Flux_heat_surf_tmp !< temporary surface heat flux
integer :: i !< counter integer :: i !< counter
integer :: status, num !< for file input integer :: status, num !< for file input
real :: Flux_heat_surf_time_forcing(nt) !< surface heat flux, input: [W/m**2]
integer :: set_Flux !< how to set: 1 - const || 2 - function || 3 - file integer :: set_Flux !< how to set: 1 - const || 2 - function || 3 - file
real :: df set_Flux = 1
df = 3600.0
!integer :: nf
set_Flux = 3
!open (555, file= 'Kato-Phillips/Flux_heat_surface.txt', status ='old') !'Flux_heat_surface_PAPA.txt' !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') open (555, file= 'PAPA_06_2017/Flux_heat_surface_PAPA.txt', status ='old')
...@@ -69,49 +66,50 @@ module obl_forcing_and_boundary ...@@ -69,49 +66,50 @@ module obl_forcing_and_boundary
status = 0 status = 0
num = 0 num = 0
if (set_Flux.eq.1) then if (set_Flux.eq.1) then
do i = 1, nt do i = 1, nf
Flux_heat_surf(i) = 0.0 Flux_heat_surf(i) = 0.0
Time_flux_heat_surf(i) = i * df
end do end do
else if (set_Flux.eq.3) then else if (set_Flux.eq.3) then
do while (status.eq.0) do while (status.eq.0)
read (555, *, iostat = status) Flux_heat_surf_tmp read (555, *, iostat = status) Flux_heat_surf_tmp
num = num + 1 num = num + 1
Flux_heat_surf(num) = - Flux_heat_surf_tmp/Rho_ref/cp_water Flux_heat_surf(num) = - Flux_heat_surf_tmp/Rho_ref/cp_water
Flux_heat_surf_time_forcing(num) = num * df Time_flux_heat_surf(num) = num * df
if (num >= nt) then if (num >= nf) then
exit exit
endif endif
end do end do
!call interpolate(Flux_heat_surf(nt), nt, Flux_heat_surf_in(num), num)
!write(*,*) num
endif endif
close (555) close (555)
end subroutine end subroutine
subroutine set_Tau_x_surf(Tau_x_surf, flux_u_surf, nt) subroutine set_Tau_x_surf(Tau_x_surf, flux_u_surf, Time_flux_u_surf, nf, df)
!< @brief set surface u-momentum flux !< @brief set surface u-momentum flux
integer, intent(in) :: nt !< number of timesteps integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: Tau_x_surf(nt) !< [N/m**2] real, intent(in) :: df
real, intent(out) :: flux_u_surf(nt) !< u_dyn**2, [(m/s)**2] (?) real, intent(out) :: Tau_x_surf(nf) !< [N/m**2]
real, intent(out) :: flux_u_surf(nf) !< u_dyn**2, [(m/s)**2] (?)
real, intent(out) :: Time_flux_u_surf(nf) !< Time array
real :: Tau_x_surf_tmp !< temporary surface u-momentum flux real :: Tau_x_surf_tmp !< temporary surface u-momentum flux
integer :: i !< counter integer :: i !< counter
integer :: status, num !< for file input integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 3 set_Tau = 1
!open (188, file= 'Kato-Phillips/Tau_x.txt', status ='old') !open (188, file= 'Kato-Phillips/Tau_x.txt', status ='old')
open (188, file= 'PAPA_06_2017/Tau_x_PAPA.txt', status ='old') open (188, file= 'PAPA_06_2017/Tau_x_PAPA.txt', status ='old')
status = 0 status = 0
num = 0 num = 0
if (set_Tau.eq.1) then if (set_Tau.eq.1) then
do i = 1, nt do i = 1, nf
Tau_x_surf(i) = 0.1 !0.1 Tau_x_surf(i) = 0.1 !0.1
end do end do
else if (set_Tau.eq.2) then else if (set_Tau.eq.2) then
read (188, *, iostat = status) Tau_x_surf_tmp read (188, *, iostat = status) Tau_x_surf_tmp
do i = 1, nt do i = 1, nf
Tau_x_surf(i) = Tau_x_surf_tmp Tau_x_surf(i) = Tau_x_surf_tmp
end do end do
else if (set_Tau.eq.3) then else if (set_Tau.eq.3) then
...@@ -119,43 +117,47 @@ module obl_forcing_and_boundary ...@@ -119,43 +117,47 @@ module obl_forcing_and_boundary
read (188, *, iostat = status) Tau_x_surf_tmp read (188, *, iostat = status) Tau_x_surf_tmp
num = num + 1 num = num + 1
Tau_x_surf(num) = Tau_x_surf_tmp Tau_x_surf(num) = Tau_x_surf_tmp
if (num >= nt) then if (num >= nf) then
exit exit
endif endif
end do end do
endif endif
close (188) close (188)
do i = 1, nt do i = 1, nf
flux_u_surf(i) = Tau_x_surf(i) / Rho_ref flux_u_surf(i) = Tau_x_surf(i) / Rho_ref
Time_flux_u_surf(i) = i * df
end do end do
end subroutine end subroutine
subroutine set_Tau_y_surf(Tau_y_surf, flux_v_surf, nt) subroutine set_Tau_y_surf(Tau_y_surf, flux_v_surf, Time_flux_v_surf, nf, df)
!< @brief set surface v-momentum flux !< @brief set surface v-momentum flux
integer, intent(in) :: nt !< number of timesteps
real, intent(out) :: Tau_y_surf(nt) !< [N/m**2] integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: flux_v_surf(nt) !< v_dyn**2, [(m/s)**2] real, intent(in) :: df
real, intent(out) :: Tau_y_surf(nf) !< [N/m**2]
real, intent(out) :: flux_v_surf(nf) !< v_dyn**2, [(m/s)**2]
real, intent(out) :: Time_flux_v_surf(nf) !< Time array
real :: Tau_y_surf_tmp !< temporary surface v-momentum flux real :: Tau_y_surf_tmp !< temporary surface v-momentum flux
integer :: i !< counter integer :: i !< counter
integer :: status, num !< for file input integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 3 set_Tau = 1
!open (177, file= 'Kato-Phillips/Tau_y.txt', status ='old') !open (177, file= 'Kato-Phillips/Tau_y.txt', status ='old')
open (177, file= 'PAPA_06_2017/Tau_y_PAPA.txt', status ='old') open (177, file= 'PAPA_06_2017/Tau_y_PAPA.txt', status ='old')
status = 0 status = 0
num = 0 num = 0
if (set_Tau.eq.1) then if (set_Tau.eq.1) then
do i = 1, nt do i = 1, nf
Tau_y_surf(i) = 0.0 Tau_y_surf(i) = 0.0
end do end do
else if (set_Tau.eq.2) then else if (set_Tau.eq.2) then
read (177, *, iostat = status) Tau_y_surf_tmp read (177, *, iostat = status) Tau_y_surf_tmp
do i = 1, nt do i = 1, nf
Tau_y_surf(i) = Tau_y_surf_tmp Tau_y_surf(i) = Tau_y_surf_tmp
end do end do
else if (set_Tau.eq.3) then else if (set_Tau.eq.3) then
...@@ -163,85 +165,98 @@ module obl_forcing_and_boundary ...@@ -163,85 +165,98 @@ module obl_forcing_and_boundary
read (177, *, iostat = status) Tau_y_surf_tmp read (177, *, iostat = status) Tau_y_surf_tmp
num = num + 1 num = num + 1
Tau_y_surf(num) = Tau_y_surf_tmp Tau_y_surf(num) = Tau_y_surf_tmp
if (num >= nt) then if (num >= nf) then
exit exit
endif endif
end do end do
endif endif
close (177) close (177)
do i = 1, nt do i = 1, nf
flux_v_surf(i) = Tau_y_surf(i) / Rho_ref flux_v_surf(i) = Tau_y_surf(i) / Rho_ref
Time_flux_v_surf(i) = i * df
end do end do
end subroutine end subroutine
subroutine set_Tau_x_bot(Tau_x_bot, flux_u_bot, nt) subroutine set_Tau_x_bot(Tau_x_bot, flux_u_bot, Time_flux_u_bot, nf, df)
!< @brief set bottom u-momentum flux !< @brief set bottom u-momentum flux
integer, intent(in) :: nt !< number of timesteps integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: Tau_x_bot(nt) !< [N/m**2] real, intent(in) :: df
real, intent(out) :: flux_u_bot(nt) !< u_dyn**2, [(m/s)**2] real, intent(out) :: Tau_x_bot(nf) !< [N/m**2]
real :: Tau_x_bot_tmp !< temporary bottom u-momentum flux real, intent(out) :: flux_u_bot(nf) !< u_dyn**2, [(m/s)**2] (?)
real, intent(out) :: Time_flux_u_bot(nf) !< Time array
real :: Tau_x_bot_tmp !< temporary surface u-momentum flux
integer :: i !< counter integer :: i !< counter
integer :: status, num !< for file input integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 1 set_Tau = 1
do i = 1, nt do i = 1, nf
Tau_x_bot(i) = 0.0 Tau_x_bot(i) = 0.0
end do end do
do i = 1, nt do i = 1, nf
flux_u_bot(i) = Tau_x_bot(i) / Rho_ref flux_u_bot(i) = Tau_x_bot(i) / Rho_ref
Time_flux_u_bot(i) = i * df
end do end do
end subroutine end subroutine
subroutine set_Tau_y_bot(Tau_y_bot, flux_v_bot, nt) subroutine set_Tau_y_bot(Tau_y_bot, flux_v_bot, Time_flux_v_bot, nf, df)
!< @brief set bottom v-momentum flux !< @brief set bottom v-momentum flux
integer, intent(in) :: nt !< number of timesteps integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: Tau_y_bot(nt) !< [N/m**2] real, intent(in) :: df
real, intent(out) :: flux_v_bot(nt) !< v_dyn**2, [(m/s)**2] real, intent(out) :: Tau_y_bot(nf) !< [N/m**2]
real :: Tau_y_bot_tmp !< temporary bottom v-momentum flux real, intent(out) :: flux_v_bot(nf) !< u_dyn**2, [(m/s)**2] (?)
real, intent(out) :: Time_flux_v_bot(nf) !< Time array
real :: Tau_y_bot_tmp !< temporary surface u-momentum flux
integer :: i !< counter integer :: i !< counter
integer :: status, num !< for file input integer :: status, num !< for file input
integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
set_Tau = 1 set_Tau = 1
do i = 1, nt do i = 1, nf
Tau_y_bot(i) = 0.0 Tau_y_bot(i) = 0.0
end do end do
do i = 1, nt do i = 1, nf
flux_v_bot(i) = Tau_y_bot(i) / Rho_ref flux_v_bot(i) = Tau_y_bot(i) / Rho_ref
Time_flux_v_bot(i) = i * df
end do end do
end subroutine end subroutine
subroutine set_Flux_sal_surf(Flux_sal_surf, nt) subroutine set_Flux_sal_surf(Flux_sal_surf, Time_flux_sal_surf, nf, df)
!< @brief setsurface salinity flux !< @brief setsurface salinity flux
integer, intent(in) :: nt !< number of timesteps integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: Flux_sal_surf(nt) ! [PSU*m/s] ??? real, intent(in) :: df
real, intent(out) :: Flux_sal_surf(nf) ! [PSU*m/s] ???
real, intent(out) :: Time_flux_sal_surf(nf) !< Time array
integer :: i !< counter integer :: i !< counter
do i = 1, nt do i = 1, nf
Flux_sal_surf(i) = 0.0 Flux_sal_surf(i) = 0.0
Time_flux_sal_surf(i) = i * df
end do end do
end subroutine end subroutine
subroutine set_Flux_sal_bot(Flux_sal_bot, nt) subroutine set_Flux_sal_bot(Flux_sal_bot, Time_flux_sal_bot, nf, df)
!< @brief set bottom salinity flux !< @brief set bottom salinity flux
integer, intent(in) :: nt !< number of timesteps integer, intent(in) :: nf !< number of timesteps
real, intent(out) :: Flux_sal_bot(nt) ! [???] real, intent(in) :: df
real, intent(out) :: Flux_sal_bot(nf) ! [???]
real, intent(out) :: Time_flux_sal_bot(nf) !< Time array
integer :: i !< counter integer :: i !< counter
do i = 1, nt do i = 1, nf
Flux_sal_bot(i) = 0.0 Flux_sal_bot(i) = 0.0
Time_flux_sal_bot(i) = i * df
end do end do
end subroutine end subroutine
......
obl_main.f90
+ 171
− 130
View file @ 0d97260d
...@@ -47,19 +47,22 @@ program obl_main ...@@ -47,19 +47,22 @@ program obl_main
real, allocatable :: Km(:) !< eddy viscosity for momentum, [m**2 / s] real, allocatable :: Km(:) !< eddy viscosity for momentum, [m**2 / s]
real, allocatable :: Kh(:) !eddy diffusivity for tracers, [m**2 / s] real, allocatable :: Kh(:) !eddy diffusivity for tracers, [m**2 / s]
real, allocatable :: Flux_heat_surf(:), Flux_heat_bot(:) !< heat fluxes, [K*m/s] real, allocatable :: Flux_heat_surf(:), Flux_heat_bot(:) !< heat fluxes, [K*m/s]
real, allocatable :: Flux_heat_surf_time_forcing(:) real, allocatable :: Times_flux_Heat_surf(:), Times_flux_Heat_bot(:)
!real :: Flux_heat_surf_time_forcing(:)
real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:) !< salt fluxes, [PSU*m/s] ??? real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:) !< salt fluxes, [PSU*m/s] ???
real, allocatable :: Times_flux_Sal_surf(:), Times_flux_Sal_bot(:)
real, allocatable :: Tau_x_surf(:), Tau_y_surf(:), Tau_x_bot(:), Tau_y_bot(:) !< [N/m**2] real, allocatable :: Tau_x_surf(:), Tau_y_surf(:), Tau_x_bot(:), Tau_y_bot(:) !< [N/m**2]
real, allocatable :: flux_u_surf(:), flux_u_bot(:) !< [(m/s)**2] real, allocatable :: flux_u_surf(:), flux_u_bot(:) !< [(m/s)**2]
real, allocatable :: Times_flux_u_surf(:), Times_flux_u_bot(:)
real, allocatable :: flux_v_surf(:), flux_v_bot(:) !< [(m/s)**2] real, allocatable :: flux_v_surf(:), flux_v_bot(:) !< [(m/s)**2]
real, allocatable :: Times_flux_v_surf(:), Times_flux_v_bot(:)
real, allocatable :: f_Heat_right(:), f_Sal_right(:), f_u_right(:), f_v_right(:) !< RHS real, allocatable :: f_Heat_right(:), f_Sal_right(:), f_u_right(:), f_v_right(:) !< RHS
integer :: closure_mode !< closure type: integer :: closure_mode !< closure type:
!1 - pacanowski-philander, 2 - pacanowski-philander+, !1 - pacanowski-philander, 2 - pacanowski-philander+,
!3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom !3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
real, allocatable :: Rho(:) !< density, [kg / m**3] real, allocatable :: Rho(:) !< density, [kg / m**3]
real, allocatable :: Rho_dev(:) !< density, [kg / m**3] real, allocatable :: Rho_dev(:) !< density, [kg / m**3]
real, allocatable :: Rho_dyn_surf(:), Rho_dyn_bot(:) !< density, [kg / m**3] real :: Rho_dyn_surf, Rho_dyn_bot !< density dynamic, [kg / m**3]
real, allocatable :: N2(:) !< square of buoyancy frequency, [s**(-2)] real, allocatable :: N2(:) !< square of buoyancy frequency, [s**(-2)]
real, allocatable :: S2(:) !< square of shear frequency, [s**(-2)] real, allocatable :: S2(:) !< square of shear frequency, [s**(-2)]
real, allocatable :: Ri_grad(:) !< Richardson gradient number real, allocatable :: Ri_grad(:) !< Richardson gradient number
...@@ -73,6 +76,11 @@ program obl_main ...@@ -73,6 +76,11 @@ program obl_main
real, allocatable :: Km_TKE(:) !eddy viscosity for momentum adjusted for Schmidt TKE number, [m**2 / s] real, allocatable :: Km_TKE(:) !eddy viscosity for momentum adjusted for Schmidt TKE number, [m**2 / s]
real, allocatable :: Km_eps(:) !eddy viscosity for momentum adjusted for Schmidt eps number, [m**2 / s] real, allocatable :: Km_eps(:) !eddy viscosity for momentum adjusted for Schmidt eps number, [m**2 / s]
real :: flux_heat_surf_res, flux_heat_bot_res, flux_sal_surf_res, flux_sal_bot_res
real :: flux_u_surf_res, flux_u_bot_res, flux_v_surf_res, flux_v_bot_res
! arrays for output ! arrays for output
real, allocatable :: Theta_write(:,:) real, allocatable :: Theta_write(:,:)
real, allocatable :: Salin_write(:,:) real, allocatable :: Salin_write(:,:)
...@@ -89,12 +97,13 @@ program obl_main ...@@ -89,12 +97,13 @@ program obl_main
real, allocatable :: Theta_C_write(:,:) real, allocatable :: Theta_C_write(:,:)
real, allocatable :: lab_hml_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 closure_mode = 4 !< 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
! set time, depth, steps and number of space ! set time, depth, steps and number of space
call set_Time_Space(t, nt, z, nz, nf) call set_Time_Space(t, nt, z, nz)
allocate(dz(1:nz)) allocate(dz(1:nz))
call set_dz_dt(dz, dt, nz, nt, t, z, df) call set_dz_dt(dz, dt, nz, nt, t, z)
call set_Time_Space_Forcing(nf, df)
write (*,*) t, dt, nt, z, nz, nf, df write (*,*) t, dt, nt, z, nz, nf, df
! allocation ! allocation
...@@ -105,28 +114,33 @@ program obl_main ...@@ -105,28 +114,33 @@ program obl_main
allocate(U(1:nz)) allocate(U(1:nz))
allocate(V(1:nz)) allocate(V(1:nz))
allocate(flux_Heat_surf(1:nf)) allocate(flux_Heat_surf(1:nf))
!allocate(flux_Heat_surf_time_forcing(1:nt)) allocate(Times_flux_Heat_surf(1:nf))
allocate(flux_Heat_bot(1:nt)) allocate(flux_Heat_bot(1:nf))
allocate(flux_Sal_surf(1:nt)) allocate(Times_flux_Heat_bot(1:nf))
allocate(flux_Sal_bot(1:nt)) allocate(flux_Sal_surf(1:nf))
allocate(Tau_x_surf(1:nt)) allocate(Times_flux_Sal_surf(1:nf))
allocate(Tau_y_surf(1:nt)) allocate(flux_Sal_bot(1:nf))
allocate(Tau_x_bot(1:nt)) allocate(Times_flux_Sal_bot(1:nf))
allocate(Tau_y_bot(1:nt)) allocate(Tau_x_surf(1:nf))
allocate(flux_u_surf(1:nt)) allocate(Tau_y_surf(1:nf))
allocate(flux_u_bot(1:nt)) allocate(Tau_x_bot(1:nf))
allocate(flux_v_surf(1:nt)) allocate(Tau_y_bot(1:nf))
allocate(flux_v_bot(1:nt)) allocate(flux_u_surf(1:nf))
allocate(f_Heat_right(1:nt)) allocate(Times_flux_u_surf(1:nf))
allocate(f_Sal_right(1:nt)) allocate(flux_u_bot(1:nf))
allocate(Times_flux_u_bot(1:nf))
allocate(flux_v_surf(1:nf))
allocate(Times_flux_v_surf(1:nf))
allocate(flux_v_bot(1:nf))
allocate(Times_flux_v_bot(1:nf))
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(Km(1:nz))
allocate(Kh(1:nz)) allocate(Kh(1:nz))
allocate(Rho(1:nz)) allocate(Rho(1:nz))
allocate(Rho_dev(1:nz)) allocate(Rho_dev(1:nz))
allocate(Rho_dyn_surf(nt)) ! вызывать после интерполяций!!!
allocate(Rho_dyn_bot(nt)) ! вызывать после интерполяций!!!
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))
...@@ -159,24 +173,23 @@ program obl_main ...@@ -159,24 +173,23 @@ program obl_main
call V_init(V, nz) call V_init(V, nz)
!set boundary conditions !set boundary conditions
!call set_Flux_heat_surf (Flux_heat_surf, Flux_heat_surf_time_forcing, n) call set_Flux_heat_surf (Flux_heat_surf, Times_flux_heat_surf, nf, df)
call set_Flux_heat_bot (Flux_heat_bot, nf) call set_Flux_heat_bot (Flux_heat_bot, Times_flux_heat_bot, nf, df)
call set_Flux_sal_surf (Flux_sal_surf, nf) call set_Flux_sal_surf (Flux_sal_surf, Times_flux_sal_surf, nf, df)
call set_Flux_sal_bot (Flux_sal_bot, nf) call set_Flux_sal_bot (Flux_sal_bot, Times_flux_sal_bot, nf, df)
call set_Tau_x_surf(Tau_x_surf, flux_u_surf, nf) call set_Tau_x_surf(Tau_x_surf, flux_u_surf, Times_flux_u_surf, nf, df)
call set_Tau_y_surf(Tau_y_surf, flux_v_surf, nf) call set_Tau_y_surf(Tau_y_surf, flux_v_surf, Times_flux_v_surf, nf, df)
call set_Tau_x_bot(Tau_x_bot, flux_u_bot, nf) call set_Tau_x_bot(Tau_x_bot, flux_u_bot, Times_flux_u_bot, nf, df)
call set_Tau_y_bot(Tau_y_bot, flux_v_bot, nf) call set_Tau_y_bot(Tau_y_bot, flux_v_bot, Times_flux_v_bot, nf, df)
call set_f_Heat_right (f_Heat_right, nz) call set_f_Heat_right (f_Heat_right, nz)
call set_f_Sal_right (f_Sal_right, nz) call set_f_Sal_right (f_Sal_right, nz)
call set_f_u_right (f_v_right, V, nz) call set_f_u_right (f_v_right, V, nz)
call set_f_v_right (f_v_right, U, nz) call set_f_v_right (f_v_right, U, nz)
call get_rho_dyn(Rho_dyn_surf, flux_u_surf, flux_v_surf, flux_heat_surf, flux_sal_surf, nt) ! как быть с числом шагов???
call get_rho_dyn(Rho_dyn_bot, flux_u_bot, flux_v_bot, flux_heat_bot, flux_sal_bot, nt)
open (199, file= 'output_Daria/hml.txt', status ='replace')
open (20, file= 'output_Daria/surf_temp.txt', status ='replace') !open (199, file= 'output_Daria/hml.txt', status ='replace')
!open (20, file= 'output_Daria/surf_temp.txt', status ='replace')
status = 0 status = 0
num = 0 num = 0
...@@ -198,75 +211,43 @@ program obl_main ...@@ -198,75 +211,43 @@ program obl_main
i=1 i=1
do while (time_current < time_end ) !.and. i<=nt do while (time_current < time_end )
if (closure_mode.eq.1) then if (closure_mode.ne.5) then
call get_value_interpolate(flux_heat_surf_res, flux_heat_surf, Times_flux_heat_surf, nf, time_current)
call get_value_interpolate(flux_heat_bot_res, flux_heat_bot, Times_flux_heat_bot, nf, time_current)
call get_value_interpolate(flux_sal_surf_res, flux_sal_surf, Times_flux_sal_surf, nf, time_current)
call get_value_interpolate(flux_sal_bot_res, flux_sal_bot, Times_flux_sal_bot, nf, time_current)
call get_value_interpolate(flux_u_surf_res, flux_u_surf, Times_flux_u_surf, nf, time_current)
call get_value_interpolate(flux_v_surf_res, flux_v_surf, Times_flux_v_surf, nf, time_current)
call get_value_interpolate(flux_u_bot_res, flux_u_bot, Times_flux_u_surf, nf, time_current)
call get_value_interpolate(flux_v_bot_res, flux_v_bot, Times_flux_v_bot, nf, time_current)
call get_rho_dyn(Rho_dyn_surf, flux_u_surf_res, flux_v_surf_res, flux_heat_surf_res, flux_sal_surf_res)
call get_rho_dyn(Rho_dyn_bot, flux_u_bot_res, flux_v_bot_res, flux_heat_bot_res, flux_sal_bot_res)
call solve_state_eq (Rho, Theta_dev, Salin_dev, nz) call solve_state_eq (Rho, Theta_dev, Salin_dev, nz)
!call get_value_interpolate(flux_heat_surf(time_current), flux_heat_surf, Flux_heat_surf_time_forcing, nf, time_current) call get_N2(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, dz, nz)
call get_N2(N2, Rho, Rho_dyn_surf(i), Rho_dyn_bot(i), dz, nz) call get_S2 (S2, U, V, flux_u_surf_res, flux_v_surf_res, flux_u_bot_res, flux_v_bot_res, dz, nz)
call get_S2 (S2, U, V, flux_u_surf(i), flux_v_surf(i), flux_u_bot(i), flux_v_bot(i), 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_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_res, flux_v_surf_res, time_current, z)
if (closure_mode.eq.1) then
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 solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq (Salin_dev, 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_dev, Salin_dev, nz)
call get_N2(N2, Rho, Rho_dyn_surf(i), Rho_dyn_bot(i), dz, nz)
call get_S2 (S2, U, V, flux_u_surf(i), flux_v_surf(i), flux_u_bot(i), flux_v_bot(i), 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_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_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right) else if (closure_mode.eq.2) then
call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right) call get_Km_plus(Km, Ri_grad, flux_u_surf_res, flux_v_surf_res, hml, nz, dz, z)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right) call get_Kh_plus(Kh, Km, nz)
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_dev, Salin_dev, nz)
call get_N2(N2, Rho, Rho_dyn_surf(i), Rho_dyn_bot(i), dz, nz)
call get_S2 (S2, U, V, flux_u_surf(i), flux_v_surf(i), flux_u_bot(i), flux_v_bot(i), 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_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 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 limit_min_array(TKE, TKE_min, nz)
call solve_eps_eq (eps, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
call limit_min_array(eps, eps_min, nz)
call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_scalar_eq (Salin_dev, 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 else if (closure_mode.eq.4) then
!call solve_state_eq_dev (Rho_dev, Theta_dev, Salin_dev, nz) call TKE_bc(TKE, flux_u_surf_res, flux_v_surf_res, nz)
call solve_state_eq (Rho, Theta_dev, Salin_dev, nz) call eps_bc(eps, flux_u_surf_res, flux_v_surf_res, nz, dz)
call get_N2(N2, Rho, Rho_dyn_surf(i), Rho_dyn_bot(i), dz, nz)
call get_S2 (S2, U, V, flux_u_surf(i), flux_v_surf(i), flux_u_bot(i), flux_v_bot(i), 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_Km_k_eps (Km, TKE, eps, nz)
call get_Kh_k_eps (Kh, Km, nz) call get_Kh_k_eps (Kh, Km, 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 TKE_bc(TKE, flux_u_surf(i), flux_v_surf(i), nz) call TKE_bc(TKE, flux_u_surf_res, flux_v_surf_res, nz)
call eps_bc(eps, flux_u_surf(i), flux_v_surf(i), nz, dz) call eps_bc(eps, flux_u_surf_res, flux_v_surf_res, nz, dz)
call get_Km_TKE(Km_TKE, Km, nz) call get_Km_TKE(Km_TKE, Km, nz)
call get_Km_eps(Km_eps, 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_TKE_eq_semiimplicit (TKE, Km_TKE, nz, dz, dt, P_TKE, B_TKE, eps)
...@@ -274,17 +255,18 @@ program obl_main ...@@ -274,17 +255,18 @@ program obl_main
call limit_min_array(TKE, TKE_min, nz) call limit_min_array(TKE, TKE_min, nz)
call solve_eps_eq_semiimplicit (eps, Km_eps, nz, dz, dt, P_TKE, B_TKE, TKE) 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 limit_min_array(eps, eps_min, nz)
call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right) endif
call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right) call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf_res, flux_heat_bot_res, f_heat_right)
call solve_vector_eq (U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right) call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf_res, flux_sal_bot_res, f_sal_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 (U, Km, nz, dz, dt, flux_u_surf_res, flux_u_bot_res, f_u_right)
else if (closure_mode.eq.5) then call solve_vector_eq (V, Km, nz, dz, dt, flux_v_surf_res, flux_v_bot_res, f_v_right)
call solve_state_eq (Rho, Theta_dev, Salin_dev, nz) !else if (closure_mode.eq.5) then
call vermix_init(rho, u, v, Tau_x_surf(i), Tau_y_surf(i), Ri_grad, Kh, Km) ! call solve_state_eq (Rho, Theta_dev, Salin_dev, nz)
call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right) ! call vermix_init(rho, u, v, Tau_x_surf(i), Tau_y_surf(i), Ri_grad, Kh, Km)
call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_right) ! call solve_scalar_eq (Theta_dev, Kh, nz, dz, dt, flux_heat_surf(i), flux_heat_bot(i), f_heat_right)
call solve_vector_eq(U, Km, nz, dz, dt, flux_u_surf(i), flux_u_bot(i), f_u_right) ! call solve_scalar_eq (Salin_dev, Kh, nz, dz, dt, flux_sal_surf(i), flux_sal_bot(i), f_sal_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(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)
endif endif
!do k = 1, nz !do k = 1, nz
...@@ -305,12 +287,11 @@ program obl_main ...@@ -305,12 +287,11 @@ program obl_main
maxcellnumber_write(i) = maxcellnumber maxcellnumber_write(i) = maxcellnumber
lab_hml_write(i) = lab_hml lab_hml_write(i) = lab_hml
write(199,*) time_current, hml, lab_hml !write(199,*) time_current, hml, lab_hml
write(20,*) time_current, Theta_dev(nz) + Theta_ref !write(20,*) time_current, Theta_dev(nz) + Theta_ref
i = i+1 i = i+1
time_current = time_current+dt time_current = time_current+dt
enddo enddo
Theta_C_write = Theta_write - 273.15 Theta_C_write = Theta_write - 273.15
...@@ -334,8 +315,8 @@ program obl_main ...@@ -334,8 +315,8 @@ program obl_main
! 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_surf, '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_surf, '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
...@@ -354,11 +335,71 @@ program obl_main ...@@ -354,11 +335,71 @@ program obl_main
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_surf, '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_surf, 'output_tau_v.txt', meta_tau_v) !call write_1d_real_ascii(Tau_y_surf, 'output_tau_v.txt', meta_tau_v)
!close(199)
!close(20)
close(199) ! deallocation
close(20) deallocate(dz)
deallocate(Theta)
deallocate(Theta_dev)
deallocate(Salin)
deallocate(Salin_dev)
deallocate(U)
deallocate(V)
deallocate(flux_Heat_surf)
deallocate(Times_flux_Heat_surf)
deallocate(flux_Heat_bot)
deallocate(Times_flux_Heat_bot)
deallocate(flux_Sal_surf)
deallocate(Times_flux_Sal_surf)
deallocate(flux_Sal_bot)
deallocate(Times_flux_Sal_bot)
deallocate(Tau_x_surf)
deallocate(Tau_y_surf)
deallocate(Tau_x_bot)
deallocate(Tau_y_bot)
deallocate(flux_u_surf)
deallocate(Times_flux_u_surf)
deallocate(flux_u_bot)
deallocate(Times_flux_u_bot)
deallocate(flux_v_surf)
deallocate(Times_flux_v_surf)
deallocate(flux_v_bot)
deallocate(Times_flux_v_bot)
deallocate(f_Heat_right)
deallocate(f_Sal_right)
deallocate(f_u_right)
deallocate(f_v_right)
deallocate(Km)
deallocate(Kh)
deallocate(Rho)
deallocate(Rho_dev)
deallocate(N2)
deallocate(S2)
deallocate(Ri_grad)
deallocate(TKE)
deallocate(eps)
deallocate(P_TKE)
deallocate(B_TKE)
deallocate(Km_TKE)
deallocate(Km_eps)
deallocate(Theta_write)
deallocate(Salin_write)
deallocate(U_write)
deallocate(V_write)
deallocate(Rho_write)
deallocate(N2_write)
deallocate(S2_write)
deallocate(Ri_grad_write)
deallocate(Kh_write)
deallocate(Km_write)
deallocate(hml_write)
deallocate(maxcellnumber_write)
deallocate(Theta_C_write)
deallocate(lab_hml_write)
end program end program
\ No newline at end of file
obl_richardson.f90
+ 0
− 5
View file @ 0d97260d
...@@ -17,11 +17,6 @@ module obl_richardson ...@@ -17,11 +17,6 @@ module obl_richardson
contains contains
!subroutine get_u_star(u_star, flux_u_surf, flux_v_surf, nz)
! real, intent(out) :: u_star(nz) !< square of buoyancy frequency, [s**(-2)]
!
!end subroutine
subroutine get_N2(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, dz, nz) subroutine get_N2(N2, Rho, Rho_dyn_surf, Rho_dyn_bot, dz, nz)
!< @brief calculate square of buoyancy frequency !< @brief calculate square of buoyancy frequency
......
obl_state_eq.f90
+ 9
− 13
View file @ 0d97260d
...@@ -46,22 +46,18 @@ module obl_state_eq ...@@ -46,22 +46,18 @@ module obl_state_eq
end subroutine end subroutine
subroutine get_rho_dyn(Rho_dyn, flux_u, flux_v, flux_heat, flux_sal, nt) subroutine get_rho_dyn(Rho_dyn, flux_u, flux_v, flux_heat, flux_sal)
!< @brief rho_dyn for N2 calculation !< @brief rho_dyn for N2 calculation
real, intent(out) :: Rho_dyn(nt) !< density, [kg / m**3] real, intent(out) :: Rho_dyn !< density, [kg / m**3]
real, intent(inout) :: flux_u(nt), flux_v(nt), flux_heat(nt), flux_sal(nt) !< [(m/s)**2] real, intent(inout) :: flux_u, flux_v, flux_heat, flux_sal !< [(m/s)**2]
integer, intent(in) :: nt real :: Theta_dyn, Salin_dyn
real :: Theta_dyn(nt), Salin_dyn(nt)
integer :: i call limit_min(flux_u, flux_min)
call limit_min(flux_v, flux_min)
do i=1, nt Theta_dyn = abs(flux_heat) / (flux_u**(2.0)+flux_v**(2.0))**(1.0/4.0)
call limit_min(flux_u(i), flux_min) Salin_dyn = abs(flux_sal) / (flux_u**(2.0)+flux_v**(2.0))**(1.0/4.0)
call limit_min(flux_v(i), flux_min) Rho_dyn = - alpha_state * Theta_dyn + beta_state * Salin_dyn
Theta_dyn(i) = abs(flux_heat(i)) / (flux_u(i)**(2.0)+flux_v(i)**(2.0))**(1.0/4.0)
Salin_dyn(i) = abs(flux_sal(i)) / (flux_u(i)**(2.0)+flux_v(i)**(2.0))**(1.0/4.0)
Rho_dyn(i) = - alpha_state * Theta_dyn(i) + beta_state * Salin_dyn(i)
enddo
end subroutine end subroutine
......
obl_time_and_space.f90
+ 21
− 8
View file @ 0d97260d
...@@ -7,26 +7,41 @@ module obl_time_and_space ...@@ -7,26 +7,41 @@ module obl_time_and_space
contains contains
subroutine set_Time_Space(t, nt, z, nz, nf) subroutine set_Time_Space(t, nt, z, nz)
!< @brief set time, depth and number of steps !< @brief set time, depth and number of steps
real, intent(out) :: t !< time, [s] real, intent(out) :: t !< time, [s]
integer, intent(out) :: nt !< number of time steps integer, intent(out) :: nt !< number of time steps
real, intent(out) :: z !< depth, [m] real, intent(out) :: z !< depth, [m]
integer, intent(out) :: nz !< number of steps in z (depth) integer, intent(out) :: nz !< number of steps in z (depth)
integer, intent(out) :: nf !< number of time steps for forcing
integer :: status, num !< for file input integer :: status, num !< for file input
integer :: i !< counter integer :: i !< counter
!open (1, file= 'PAPA_06_2017/time_space_for_calc_PAPA.txt', status ='old') !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') open (1, file= 'Kato-Phillips/time_space_for_calc.txt', status ='old')
status = 0 status = 0
do while (status.eq.0) do while (status.eq.0)
read (1, *, iostat = status) t, nt, z, nz, nf read (1, *, iostat = status) t, nt, z, nz
end do end do
close (1) close (1)
end subroutine end subroutine
subroutine set_dz_dt(dz, dt, nz, nt, t, z, df) subroutine set_Time_Space_Forcing(nf, df)
!< @brief set number of steps and step for forcing
integer, intent(out) :: nf !< number of time steps for forcing
real, intent(out) :: df !< time step for forcing [s]
integer :: status
!open (1, file= 'PAPA_06_2017/time_space_for_calc_PAPA.txt', status ='old')
open (115, file= 'Kato-Phillips/time_space_for_forcing.txt', status ='old')
status = 0
do while (status.eq.0)
read (115, *, iostat = status) nf, df
end do
close (115)
end subroutine
subroutine set_dz_dt(dz, dt, nz, nt, t, z)
!< @brief set time and depth steps !< @brief set time and depth steps
integer, intent(in) :: nz, nt !< number of depth and time steps integer, intent(in) :: nz, nt !< number of depth and time steps
...@@ -34,7 +49,7 @@ module obl_time_and_space ...@@ -34,7 +49,7 @@ module obl_time_and_space
real, intent(in) :: t !< time, [s] real, intent(in) :: t !< time, [s]
real, intent(out) :: dz(nz) !< depth(z) step [m] real, intent(out) :: dz(nz) !< depth(z) step [m]
real, intent(out) :: dt !< time step [s] real, intent(out) :: dt !< time step [s]
real, intent(out) :: df !< time step for forcing [s]
integer :: status, num !< for file input integer :: status, num !< for file input
real :: dz_start !< temporary dz, [m] real :: dz_start !< temporary dz, [m]
integer :: init_dz !< how to set: 1 - const || 3 - file integer :: init_dz !< how to set: 1 - const || 3 - file
...@@ -42,9 +57,7 @@ module obl_time_and_space ...@@ -42,9 +57,7 @@ module obl_time_and_space
dt = t/nt dt = t/nt
df = 3600.0 init_dz = 1 ! 1 - const, 2 - function(z), 3 - file
init_dz = 3 ! 1 - const, 2 - function(z), 3 - file
!open (33, file= 'Kato-Phillips/dz.txt', status ='old') !'dz_PAPA.txt' !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' open (33, file= 'PAPA_06_2017/dz_PAPA.txt', status ='old') !'dz_PAPA.txt'
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment