Newer
Older
Ramil Ahtamyanov
committed
module obl_inmom
implicit none
integer, parameter :: border_shift = 1
private :: border_shift
contains
subroutine init_vermix(richnum_mode, kh_km_mode, uu, vv, lu, dx, dy, dxh, dyh, hhu, hhv, hhq, zw, g, rh0, &
rit, den, rlh, taux, tauy, aice0, tt, ss, anzt, anzu, anumaxt, anumaxu, anubgrt, anubgru)
use obl_legacy
use obl_pph, only: pph_kh => get_eddy_diffusivity_vec, &
pph_km => get_eddy_viscosity_vec, &
pphParamType
use obl_pph_dyn, only: pph_dyn_kh => get_eddy_diffusivity_vec, &
pph_dyn_km => get_eddy_viscosity_vec, &
pphDynParamType
Ramil Ahtamyanov
committed
integer, intent(in) :: richnum_mode, kh_km_mode
real, intent(inout) :: anzt(:,:,:) ! kh
real, intent(inout) :: anzu(:,:,:) ! km
real, intent(inout) :: rit(:,:,:)
real, intent(inout) :: den(:,:,:)
real, intent(in) :: uu(:,:,:), vv(:,:,:)
real, intent(in) :: tt(:,:,:)
real, intent(in) :: ss(:,:,:)
real, intent(in) :: taux(:,:), tauy(:,:)
real, intent(in) :: lu(:,:)
real, intent(in) :: dx(:,:), dy(:,:), dxh(:,:), dyh(:,:)
real, intent(in) :: hhu(:,:), hhv(:,:), hhq(:,:)
real, intent(in) :: aice0(:,:)
real, intent(in) :: rlh(:,:)
real, intent(in) :: zw(:)
real, intent(in) :: g, rh0
real, intent(in) :: anumaxt ! kh_0
real, intent(in) :: anumaxu ! km_0
real, intent(in) :: anubgrt ! kh_b0
real, intent(in) :: anubgru ! km_b0
! Local variables
integer :: nx, ny, nz
integer :: i, j
! Modules types and vars
!! obl_legacy vars
Ramil Ahtamyanov
committed
real, allocatable :: u2(:,:,:), v2(:,:,:)
real, allocatable :: s2(:,:,:)
real, allocatable :: n2(:,:,:)
real, allocatable :: kh_b(:)
real, allocatable :: km_b(:)
real, allocatable :: kh(:,:,:)
real, allocatable :: km(:,:,:)
real :: kh_str, km_str
real :: kh_undimdepth, km_undimdepth
real :: kh_unstable
real :: kh_b0, km_b0
real :: kh_0, km_0
Ramil Ahtamyanov
committed
type(pphParamType) :: pphParams
real :: neutral_mld(size(rit, 1), size(rit, 2))
real :: u_dynH(size(rit, 1), size(rit, 2))
!! pph_dyn vars
type(pphDynParamType) :: pphDynParams
!! allocate obl_legacy vars
Ramil Ahtamyanov
committed
allocate(u2(size(uu, 1), size(uu, 2), size(uu, 3)))
allocate(v2(size(vv, 1), size(vv, 2), size(vv, 3)))
allocate(s2(size(u2, 1), size(u2, 2), size(u2, 3) - 1))
allocate(n2(size(den, 1), size(den, 2), size(den, 3) - 1))
allocate(kh(size(anzt, 1), size(anzt, 2), size(anzt, 3)))
allocate(km(size(anzu, 1), size(anzu, 2), size(anzu, 3)))
allocate(kh_b(size(kh, 3)))
allocate(km_b(size(km, 3)))
kh = anzt
km = anzu
kh_0 = anumaxt
km_0 = anumaxu
kh_b0 = anubgrt
km_b0 = anubgru
Ramil Ahtamyanov
committed
nx = size(rit, 1)
ny = size(rit, 2)
nz = size(rit, 3)
den = legacy_denp(tt, ss + 35.0, 0.0)
if (richnum_mode == 1) then
call legacy_u2(uu, dy, dyh, hhu, hhq, border_shift, lu, u2)
call legacy_v2(vv, dx, dxh, hhv, hhq, border_shift, lu, v2)
call legacy_s2(u2, v2, lu, s2)
call legacy_n2(den, hhq, zw, g, lu, n2)
call legacy_rit(n2, s2, border_shift, lu, rit(:,:,2:size(rit, 3)))
call legacy_rit_top(rlh, taux, tauy, border_shift, lu, rit(:,:,1))
call sync_xy_border_3d(rit)
! neutral mld & u_dynH for pph_dyn
call legacy_neutral_mld(rlh, taux, tauy, border_shift, lu, neutral_mld, u_dynH)
call sync_xy_border_2d(neutral_mld)
call sync_xy_border_2d(u_dynH)
neutral_mld = neutral_mld * 1e-2
u_dynH = u_dynH * 1e-2
!print *, "neutral_mld:", neutral_mld
Ramil Ahtamyanov
committed
end if
! obl_legacy mixing mode
if (kh_km_mode == 1) then
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
Ramil Ahtamyanov
committed
call legacy_str("kh", taux(i,j), tauy(i,j), rh0, kh_str)
call legacy_undimdepth("kh", kh_str, hhq(i,j), aice0(i,j), kh_undimdepth)
call legacy_kh_b(zw(:), hhq(i,j), kh_unstable, kh_undimdepth, kh_b0, kh_b)
call legacy_kh_unstable(tt(i,j,1), kh_unstable)
call legacy_kh(rit(i,j,:), kh_0, kh_b(:), kh_unstable, kh(i,j,:))
call legacy_str("km", taux(i,j), tauy(i,j), rh0, km_str)
call legacy_undimdepth("km", km_str, hhq(i,j), aice0(i,j), km_undimdepth)
call legacy_km_b(zw, km_unstable, km_undimdepth, km_b0, km_b)
call legacy_km(rit(i,j,:), km_0, km_b(:), km(i,j,:))
end if
end do
end do
! obl_pph mixing mode (lake constants)
else if (kh_km_mode == 2) then
pphParams%Kh_unstable = 0.05
pphParams%Km_unstable = 0.05
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
call pph_kh(kh(i,j,:), rit(i,j,:), pphParams, nz)
call pph_km(km(i,j,:), rit(i,j,:), pphParams, nz)
end if
end do
end do
kh = kh * 10000.0
km = km * 10000.0
! obl_pph mixing mode (inmom constants)
else if (kh_km_mode == 3) then
pphParams%Km_0 = 7.0 * 0.01
pphParams%Kh_0 = 5.0 * 0.01
pphParams%alpha = 5.0
pphParams%n = 2.0
pphParams%Km_b = 5e-6
pphParams%Kh_b = 0.00001
pphParams%Kh_unstable = 0.05
pphParams%Km_unstable = 0.05
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
call pph_kh(kh(i,j,:), rit(i,j,:), pphParams, nz)
call pph_km(km(i,j,:), rit(i,j,:), pphParams, nz)
Ramil Ahtamyanov
committed
end if
Ramil Ahtamyanov
committed
end do
kh = kh * 10000.0
km = km * 10000.0
! obl_pph_dyn mixing mode
else if (kh_km_mode == 4) then
pphDynParams%Kh_unstable = 0.05
pphDynParams%Km_unstable = 0.05
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
call pph_dyn_kh(kh(i,j,:), rit(i,j,:), u_dynH(i,j), neutral_mld(i,j), pphdynParams, nz)
call pph_dyn_km(km(i,j,:), rit(i,j,:), u_dynH(i,j), neutral_mld(i,j), pphdynParams, nz)
end if
end do
end do
kh = kh * 10000.0
km = km * 10000.0
! obl_pph mixing mode (dasha constants)
else if (kh_km_mode == 5) then
pphParams%Km_0 = 7.0 * 0.01
pphParams%Kh_0 = 5.0 * 0.01
pphParams%alpha = 25.0 / 7.0 !Nuzhno tak!
pphParams%Kh_unstable = 0.05
pphParams%Km_unstable = 0.05
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
call pph_kh(kh(i,j,:), rit(i,j,:), pphParams, nz)
call pph_km(km(i,j,:), rit(i,j,:), pphParams, nz)
end if
end do
end do
kh = kh * 10000.0
km = km * 10000.0
Ramil Ahtamyanov
committed
end if
! print *, "Kh first:", kh(3,3,1:4)
! print *, "Km first:", km(3,3,1:4)
! print *, "Kh last:", kh(3,3,size(kh,3)-2:size(kh,3))
! print *, "Km last:", km(3,3,size(km,3)-2:size(km,3))
anzt = kh
anzu = km
end subroutine init_vermix
subroutine sync_xy_border_3d(array)
Ramil Ahtamyanov
committed
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
implicit none
real, intent(inout) :: array(:,:,:)
integer :: nx, ny, nz, k
! Determine array dimensions
nx = size(array, 1)
ny = size(array, 2)
nz = size(array, 3)
! Update boundary points along x (first_x and end_x boundaries)
do k = 1, nz
array(1,2:ny-1,k) = array(2,2:ny-1,k) ! first_x boundary
array(nx,2:ny-1,k) = array(nx-1,2:ny-1,k) ! end_x boundary
end do
! Update boundary points along y (first_y and end_y boundaries)
do k = 1, nz
array(2:nx-1,1,k) = array(2:nx-1,2,k) ! first_y boundary
array(2:nx-1,ny,k) = array(2:nx-1,ny-1,k) ! end_y boundary
end do
! Update corner points
do k = 1, nz
array(1,1,k) = (array(2,1,k) + array(1,2,k)) / 2.0
array(1,ny,k) = (array(2,ny,k) + array(1,ny-1,k)) / 2.0
array(nx,1,k) = (array(nx-1,1,k) + array(nx,2,k)) / 2.0
array(nx,ny,k) = (array(nx-1,ny,k) + array(nx,ny-1,k)) / 2.0
end do
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
end subroutine sync_xy_border_3d
subroutine sync_xy_border_2d(array)
implicit none
real, intent(inout) :: array(:,:)
integer :: nx, ny
! Determine array dimensions
nx = size(array, 1)
ny = size(array, 2)
! Update boundary points along x (first_x and end_x boundaries)
array(1,2:ny-1) = array(2,2:ny-1) ! first_x boundary
array(nx,2:ny-1) = array(nx-1,2:ny-1) ! end_x boundary
! Update boundary points along y (first_y and end_y boundaries)
array(2:nx-1,1) = array(2:nx-1,2) ! first_y boundary
array(2:nx-1,ny) = array(2:nx-1,ny-1) ! end_y boundary
! Update corner points
array(1,1) = (array(2,1) + array(1,2)) / 2.0
array(1,ny) = (array(2,ny) + array(1,ny-1)) / 2.0
array(nx,1) = (array(nx-1,1) + array(nx,2)) / 2.0
array(nx,ny) = (array(nx-1,ny) + array(nx,ny-1)) / 2.0
end subroutine sync_xy_border_2d
Ramil Ahtamyanov
committed
end module obl_inmom