interp

Kato-Phillips/time_space_for_calc.txt

-1080000	1080000	100	32
\ No newline at end of file
+1080000 1080000	100	32 10800
\ No newline at end of file

obl_forcing_and_boundary.f90

@@ -4,6 +4,7 @@ module obl_forcing_and_boundary
     ! modules used
     use obl_common_phys_parameters
     use obl_state_eq
+    use obl_mathematics
     
     ! directives list
     implicit none
@@ -42,20 +43,29 @@ module obl_forcing_and_boundary
         close (444)
     end subroutine
 
-    subroutine set_Flux_heat_surf (Flux_heat_surf, nt)
+    subroutine set_Flux_heat_surf (Flux_heat_surf, Flux_heat_surf_time_forcing, nt)
         !< @brief set surface heat flux
 
         integer, intent(in) :: nt !< number of timesteps
-        real, intent(out) :: Flux_heat_surf(nt) !< surface heat flux, input: [W/m**2], output: [K*m/s] 
+        real, intent(out) :: Flux_heat_surf(nt) !< surface heat flux, output: [K*m/s] 
         real :: Flux_heat_surf_tmp !< temporary surface heat flux
         integer :: i !< counter 
         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
 
-        set_Flux = 1
+        real :: df
+
+        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= 'PAPA_06_2017/Flux_heat_surface_PAPA.txt', status ='old')
+
         status = 0
         num = 0
         if (set_Flux.eq.1) then
@@ -67,10 +77,13 @@ module obl_forcing_and_boundary
                 read (555, *, iostat = status) Flux_heat_surf_tmp
                 num = num + 1    
                 Flux_heat_surf(num) = - Flux_heat_surf_tmp/Rho_ref/cp_water
+                Flux_heat_surf_time_forcing(num) = num * df 
                     if (num >= nt) then
                         exit
                     endif
             end do
+            !call interpolate(Flux_heat_surf(nt), nt, Flux_heat_surf_in(num), num)
+            !write(*,*) num
         endif
         close (555)
     end subroutine
@@ -86,7 +99,7 @@ module obl_forcing_and_boundary
         integer :: status, num !< for file input
         integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
 
-        set_Tau = 1
+        set_Tau = 3
 
         !open (188, file= 'Kato-Phillips/Tau_x.txt', status ='old')
         open (188, file= 'PAPA_06_2017/Tau_x_PAPA.txt', status ='old')
@@ -130,7 +143,7 @@ module obl_forcing_and_boundary
         integer :: status, num !< for file input
         integer :: set_Tau !< how to set: 1 - const || 2 - function || 3 - file
 
-        set_Tau = 1
+        set_Tau = 3
 
         !open (177, file= 'Kato-Phillips/Tau_y.txt', status ='old')
         open (177, file= 'PAPA_06_2017/Tau_y_PAPA.txt', status ='old')

obl_main.f90

@@ -29,9 +29,11 @@ program obl_main
 
     real :: t !< time, [s]
     integer :: nt !< number of time steps 
+    integer :: nf !< number of time steps for forcing
     real :: z !< depth, [m]
     integer :: nz !< number of steps in z (depth), [m]
     real :: dt !< time step [s] 
+    real :: df !< time step for forcing [s] 
     real, allocatable :: dz(:) !< depth(z) step [m] 
     integer :: i, k  !< counters          
     integer :: status, num !< for file input/output
@@ -45,6 +47,8 @@ program obl_main
     real, allocatable :: Km(:) !< eddy viscosity for momentum, [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_time_forcing(:)
+    !real :: Flux_heat_surf_time_forcing(:)
     real, allocatable :: Flux_sal_surf(:), Flux_sal_bot(:) !< salt fluxes, [PSU*m/s] ???
     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]          
@@ -88,10 +92,10 @@ program obl_main
     closure_mode = 1 !< 1 - pacanowski-philander, 2 - pacanowski-philander+, 3 - k-epsilon explicit, 4 - k-epsilon semiimplicit, 5 - inmom
     
     ! set time, depth, steps and number of space
-    call set_Time_Space(t, nt, z, nz) 
+    call set_Time_Space(t, nt, z, nz, nf) 
     allocate(dz(1:nz))
-    call set_dz_dt(dz, dt, nz, nt, t, z)
-    write (*,*) t, dt, nt, z, nz
+    call set_dz_dt(dz, dt, nz, nt, t, z, df)
+    write (*,*) t, dt, nt, z, nz, nf, df
 
     ! allocation
     allocate(Theta(1:nz))
@@ -100,7 +104,8 @@ program obl_main
     allocate(Salin_dev(1:nz))
     allocate(U(1:nz))
     allocate(V(1:nz))
-    allocate(flux_Heat_surf(1:nt))
+    allocate(flux_Heat_surf(1:nf))
+    !allocate(flux_Heat_surf_time_forcing(1:nt))
     allocate(flux_Heat_bot(1:nt))
     allocate(flux_Sal_surf(1:nt))
     allocate(flux_Sal_bot(1:nt))
@@ -112,16 +117,16 @@ program obl_main
     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_Heat_right(1:nt))
+    allocate(f_Sal_right(1:nt))
     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(Rho_dev(1:nz))
-    allocate(Rho_dyn_surf(nt))
-    allocate(Rho_dyn_bot(nt))
+    allocate(Rho_dyn_surf(nt))  ! вызывать после интерполяций!!!
+    allocate(Rho_dyn_bot(nt)) ! вызывать после интерполяций!!!
     allocate(N2(1:nz))
     allocate(S2(1:nz))  
     allocate(Ri_grad(1:nz)) 
@@ -154,20 +159,20 @@ program obl_main
     call V_init(V, nz)
 
     !set boundary conditions
-    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_Flux_sal_bot (Flux_sal_bot, 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_Flux_heat_surf (Flux_heat_surf, Flux_heat_surf_time_forcing, n)
+    call set_Flux_heat_bot (Flux_heat_bot, nf) 
+    call set_Flux_sal_surf (Flux_sal_surf, nf)
+    call set_Flux_sal_bot (Flux_sal_bot, nf)
+    call set_Tau_x_surf(Tau_x_surf, flux_u_surf, nf)
+    call set_Tau_y_surf(Tau_y_surf, flux_v_surf, nf)
+    call set_Tau_x_bot(Tau_x_bot, flux_u_bot, nf)
+    call set_Tau_y_bot(Tau_y_bot, flux_v_bot, nf)
     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)
     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_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')
@@ -196,6 +201,7 @@ program obl_main
     do while (time_current < time_end ) !.and. i<=nt
         if (closure_mode.eq.1) then
             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(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)

obl_mathematics.f90

@@ -6,6 +6,30 @@ module obl_mathematics
 
     contains
 
+    subroutine get_value_interpolate(Value_res, Value, Time, n, time_current) 
+        !< @brief linear interpolation
+        real, intent(in) :: time_current
+        integer, intent(in) :: n !< number of z-steps
+        real, intent(in) :: Value(n) !< input array
+        real, intent(in) :: Time(n) !< input array
+        real, intent(out) :: Value_res !< output array
+
+        integer :: i
+
+        if (time_current < Time(1)) then
+            Value_res = Value(1)
+        else if (time_current > Time(n)) then
+            Value_res = Value(n)
+        else
+            do i=1, n-1 
+                if (time_current >= Time(i).and.time_current <= Time(i+1)) then
+                    Value_res = Value(i) + (Value(i+1) - Value(i))/(Time(i+1) - Time(i)) * (time_current - Time(i))
+                    exit
+                endif
+            enddo
+        endif
+    end subroutine get_value_interpolate
+
     subroutine limit_min_array(Value, Value_min, nz) 
         !< @brief limitation function (minimum) for arrays
         integer, intent(in) :: nz !< number of z-steps

obl_time_and_space.f90

@@ -7,28 +7,26 @@ module obl_time_and_space
      
     contains
 
-    subroutine set_Time_Space(t, nt, z, nz) 
+    subroutine set_Time_Space(t, nt, z, nz, nf) 
         !< @brief set time, depth and number of steps
 
         real, intent(out) :: t !< time, [s]
         integer, intent(out) :: nt !< number of time steps 
         real, intent(out) :: z !< depth, [m]
-        integer, intent(out) :: nz !< number of steps in z (depth), [m]
+        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 :: i !< counter
-
-
         !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)
-            read (1, *, iostat = status) t, nt, z, nz
+            read (1, *, iostat = status) t, nt, z, nz, nf
         end do
         close (1)
-
     end subroutine
 
-    subroutine set_dz_dt(dz, dt, nz, nt, t, z)
+    subroutine set_dz_dt(dz, dt, nz, nt, t, z, df)
         !< @brief set time and depth steps
         
         integer, intent(in) :: nz, nt !< number of depth and time steps 
@@ -36,6 +34,7 @@ module obl_time_and_space
         real, intent(in) :: t  !< time, [s]
         real, intent(out) :: dz(nz) !< depth(z) step [m] 
         real, intent(out) :: dt !< time step [s] 
+        real, intent(out) :: df !< time step for forcing [s]
         integer :: status, num !< for file input
         real :: dz_start !< temporary dz, [m]
         integer :: init_dz  !< how to set: 1 - const || 3 - file
@@ -43,10 +42,12 @@ module obl_time_and_space
 
         dt = t/nt
 
-        init_dz = 1 ! 1 - const, 2 - function(z), 3 - file
+        df = 3600.0
     
-        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'
+        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= 'PAPA_06_2017/dz_PAPA.txt', status ='old') !'dz_PAPA.txt'
         status = 0
         num = 0
         if (init_dz.eq.1) then
@@ -64,6 +65,5 @@ module obl_time_and_space
             end do
         endif
         close (33)
-
     end subroutine
 end module
\ No newline at end of file