Newer
Older
Ramil Ahtamyanov
committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
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
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
use obl_pph_dyn
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
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
integer :: nx, ny, nz
integer :: i, j
! Modules types
type(pphParamType) :: pphParams
type(gridDataType) :: inmomGrid
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
nx = size(rit, 1)
ny = size(rit, 2)
nz = size(rit, 3)
! print *, "hello 1"
den = legacy_denp(tt, ss + 35.0, 0.0)
if (richnum_mode == 1) then
! print *, "hello rit calc"
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))
end if
! print *, "hello aft rit"
call sync_xy_border(rit)
do j = 1, ny
do i = 1, nx
if (lu(i, j) > lu_min) then
if (kh_km_mode == 1) then
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,:))
else if (kh_km_mode == 2) then
! call get_Kh(kh(i,j,:), rit(i,j,:), nz)
! call get_Km(km(i,j,:), rit(i,j,:), nz)
call get_eddy_diffusivity(kh(i,j,:), rit(i,j,:), pphParams, inmomGrid)
call get_eddy_viscosity(km(i,j,:), rit(i,j,:), pphParams, inmomGrid)
else if (kh_km_mode == 3) then
! u_dynH =
! mld =
call get_eddy_diffusivity(kh(i,j,:), rit(i,j,:), u_dynH, mld, pphParams, inmomGrid)
call get_eddy_viscosity(km(i,j,:), rit(i,j,:), u_dynH, mld, pphParams, inmomGrid)
end if
end if
end do
end do
if (kh_km_mode == 2) then
kh = kh * 10000.0
km = km * 10000.0
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(array)
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
end subroutine sync_xy_border
end module obl_inmom