Newer
Older
program CALCULATE_OBL
use EQUATION
use RICHARDSON
use PACANOWSKI_PHILANDER
use PACANOWSKI_PHILANDER_PLUS
use K_EPSILON
use io
use io_metadata
implicit none
integer :: nz !Number of steps in z (depth)
integer :: nt !Number of steps in t (time)
real :: dt !Timestep
real, allocatable :: dz(:) !z-step
real :: t !Time
real :: z !Depth
integer :: i, k
integer :: status, num
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 :: 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 :: f_heat_right(:)
real, allocatable :: Flux_heat_surf(:), flux_heat_bot(:)
real, allocatable :: Salin(:)
real, allocatable :: f_sal_right(:)
real, allocatable :: flux_sal_surf(:), flux_sal_bot(:)
real, allocatable :: Kh(:)
real, allocatable :: U(:)
real, allocatable :: V(:)
real, allocatable :: f_u_right(:)
real, allocatable :: f_v_right(:)
real, allocatable :: flux_u_surf(:), flux_u_bot(:)
real, allocatable :: flux_v_surf(:), flux_v_bot(:)
real, allocatable :: Km(:)
real :: f_cor, Ug, Vg
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
integer :: maxcellnumber
real :: ustar
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(:,:)
!call set_time_space (t, nt, z, nz, dt, dz)
!write (*,*) t, nt, z, nz, dt, dz
!Time and space
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
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(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(f_sal_right(1:nz))
allocate(flux_sal_surf(1:nt))
allocate(flux_sal_bot(1:nt))
allocate(U(1:nz))
allocate(V(1:nz))
allocate(Km(1:nz))
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_u_right(1:nz))
allocate(f_v_right(1:nz))
allocate(Rho(1:nz))
allocate(N2(1:nz))
allocate(S2(1:nz))
allocate(Ri_grad(1:nz))
allocate(tau(1:nt))
allocate(tau_x(1:nt))
allocate(tau_y(1:nt))
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))
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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
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 Salin_init(Salin, nz, dz)
!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_Tau_x_surf(Tau_x, 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
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
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)
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
do i = 1, nt
flux_u_surf(i) = tau_x(i) / Rho_0
flux_v_surf(i) = tau_y(i) / Rho_0
end do
f_cor = 0
Ug = 0
Vg = 0
do i = 1, nz
f_u_right(i) = - f_cor * (V(i) - Vg)
end do
do i = 1, nz
f_v_right(i) = + f_cor * (U(i) - Ug)
end do
!Правая часть
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
flux_sal_bot(i) = 0
end do
do i = 2, nz-1
TKE(i) = 1 / (C_mu)**(1.0/2.0) * 0.001 * 0.001
end do
do i = 2, nz-1
eps(i) = 0.001 * 0.001 * 0.001 / (kappa_k_eps * z)
end do
open (199, file= 'hml.txt', status ='replace')
open (20, file= 'surf_temp.txt', status ='replace')
open (15, file= 'output.txt', status ='replace')
status = 0
num = 0
time_begin = 0.0
time_end = t
time_current = time_begin
i=1
!time_count = 0
!time_mark = 10
do while (time_current <= time_end .and. i<=nt)
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_hml (hml, maxcellnumber, N2, dz, nz, 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_Kh(Kh, Ri_grad, nz)
!call get_TKE_generation(P_TKE, Km, S2, nz)
!call get_TKE_buoyancy(B_TKE, Kh, N2, nz)
call get_Rig(Ri_grad, N2, S2, nz)
!call get_Km_TKE(Km_TKE, Km, nz)
!call get_Km_eps(Km_eps, Km, nz)
!call solve_TKE_eq (TKE, Km, nz, dz, dt, P_TKE, B_TKE, eps)
!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_eps_eq (eps, Km, nz, dz, dt, P_TKE, B_TKE, TKE)
!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_Km(Km, Ri_grad, nz)
!call get_Kh(Kh, Ri_grad, 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)
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)
write(15,*) time_current
do k=1, nz
write(15,*) k, Theta(k), Salin(k), U(k), V(k), Kh(k), Km(k), maxcellnumber, hml, lab_hml, N2(k)
enddo
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, 'output.nc', meta_tau_u)
call write_1d_real_nc(tau_y, '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, 'output_tau_u.txt', meta_tau_u)
call write_1d_real_ascii(tau_y, 'output_tau_v.txt', meta_tau_v)
close(15)
close(199)
close(20)
end program