-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdynkeg.f90
122 lines (122 loc) · 4.34 KB
/
dynkeg.f90
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
MODULE dynkeg
USE oce
USE dom_oce
USE trd_oce
USE trddyn
USE in_out_manager
USE lbclnk
USE lib_mpp
USE prtctl
USE timing
USE bdy_oce
IMPLICIT NONE
PRIVATE
PUBLIC :: dyn_keg
INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0
INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1
REAL(KIND = wp) :: r1_48 = 1._wp / 48._wp
CONTAINS
SUBROUTINE dyn_keg(kt, kscheme)
INTEGER, INTENT( IN ) :: kt
INTEGER, INTENT( IN ) :: kscheme
INTEGER :: ji, jj, jk, jb
INTEGER :: ii, ifu, ib_bdy
INTEGER :: ij, ifv, igrd
REAL(KIND = wp) :: zu, zv
REAL(KIND = wp), DIMENSION(jpi, jpj, jpk) :: zhke
REAL(KIND = wp), ALLOCATABLE, DIMENSION(:, :, :) :: ztrdu, ztrdv
IF (ln_timing) CALL timing_start('dyn_keg')
IF (kt == nit000) THEN
IF (lwp) WRITE(numout, FMT = *)
IF (lwp) WRITE(numout, FMT = *) 'dyn_keg : kinetic energy gradient trend, scheme number=', kscheme
IF (lwp) WRITE(numout, FMT = *) '~~~~~~~'
END IF
IF (l_trddyn) THEN
ALLOCATE(ztrdu(jpi, jpj, jpk), ztrdv(jpi, jpj, jpk))
!$ACC KERNELS
ztrdu(:, :, :) = ua(:, :, :)
ztrdv(:, :, :) = va(:, :, :)
!$ACC END KERNELS
END IF
!$ACC KERNELS
zhke(:, :, jpk) = 0._wp
!$ACC END KERNELS
IF (ln_bdy) THEN
DO ib_bdy = 1, nb_bdy
IF (cn_dyn3d(ib_bdy) /= 'none') THEN
igrd = 2
DO jb = 1, idx_bdy(ib_bdy) % nblenrim(igrd)
DO jk = 1, jpkm1
ii = idx_bdy(ib_bdy) % nbi(jb, igrd)
ij = idx_bdy(ib_bdy) % nbj(jb, igrd)
ifu = NINT(idx_bdy(ib_bdy) % flagu(jb, igrd))
un(ii - ifu, ij, jk) = un(ii, ij, jk) * umask(ii, ij, jk)
END DO
END DO
igrd = 3
DO jb = 1, idx_bdy(ib_bdy) % nblenrim(igrd)
DO jk = 1, jpkm1
ii = idx_bdy(ib_bdy) % nbi(jb, igrd)
ij = idx_bdy(ib_bdy) % nbj(jb, igrd)
ifv = NINT(idx_bdy(ib_bdy) % flagv(jb, igrd))
vn(ii, ij - ifv, jk) = vn(ii, ij, jk) * vmask(ii, ij, jk)
END DO
END DO
END IF
END DO
END IF
SELECT CASE (kscheme)
CASE (nkeg_C2)
!$ACC KERNELS
DO jk = 1, jpkm1
DO jj = 2, jpj
DO ji = 2, jpi
zu = un(ji - 1, jj, jk) * un(ji - 1, jj, jk) + un(ji, jj, jk) * un(ji, jj, jk)
zv = vn(ji, jj - 1, jk) * vn(ji, jj - 1, jk) + vn(ji, jj, jk) * vn(ji, jj, jk)
zhke(ji, jj, jk) = 0.25_wp * (zv + zu)
END DO
END DO
END DO
!$ACC END KERNELS
CASE (nkeg_HW)
!$ACC KERNELS
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = 2, jpim1
zu = 8._wp * (un(ji - 1, jj, jk) * un(ji - 1, jj, jk) + un(ji, jj, jk) * un(ji, jj, jk)) + (un(ji - 1, jj - 1, jk) + un(ji - 1, jj + 1, jk)) * (un(ji - 1, jj - 1, jk) + un(ji - 1, jj + 1, jk)) + (un(ji, jj - 1, jk) + un(ji, jj + 1, jk)) * (un(ji, jj - 1, jk) + un(ji, jj + 1, jk))
zv = 8._wp * (vn(ji, jj - 1, jk) * vn(ji, jj - 1, jk) + vn(ji, jj, jk) * vn(ji, jj, jk)) + (vn(ji - 1, jj - 1, jk) + vn(ji + 1, jj - 1, jk)) * (vn(ji - 1, jj - 1, jk) + vn(ji + 1, jj - 1, jk)) + (vn(ji - 1, jj, jk) + vn(ji + 1, jj, jk)) * (vn(ji - 1, jj, jk) + vn(ji + 1, jj, jk))
zhke(ji, jj, jk) = r1_48 * (zv + zu)
END DO
END DO
END DO
!$ACC END KERNELS
CALL lbc_lnk(zhke, 'T', 1.)
END SELECT
IF (ln_bdy) THEN
!$ACC KERNELS
un(:, :, :) = un(:, :, :) * umask(:, :, :)
vn(:, :, :) = vn(:, :, :) * vmask(:, :, :)
!$ACC END KERNELS
END IF
!$ACC KERNELS
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = 2, jpim1
ua(ji, jj, jk) = ua(ji, jj, jk) - (zhke(ji + 1, jj, jk) - zhke(ji, jj, jk)) / e1u(ji, jj)
va(ji, jj, jk) = va(ji, jj, jk) - (zhke(ji, jj + 1, jk) - zhke(ji, jj, jk)) / e2v(ji, jj)
END DO
END DO
END DO
!$ACC END KERNELS
IF (l_trddyn) THEN
!$ACC KERNELS
ztrdu(:, :, :) = ua(:, :, :) - ztrdu(:, :, :)
ztrdv(:, :, :) = va(:, :, :) - ztrdv(:, :, :)
!$ACC END KERNELS
CALL trd_dyn(ztrdu, ztrdv, jpdyn_keg, kt)
DEALLOCATE(ztrdu, ztrdv)
END IF
IF (ln_ctl) CALL prt_ctl(tab3d_1 = ua, clinfo1 = ' keg - Ua: ', mask1 = umask, tab3d_2 = va, clinfo2 = ' Va: ', mask2 = vmask, clinfo3 = 'dyn')
IF (ln_timing) CALL timing_stop('dyn_keg')
END SUBROUTINE dyn_keg
END MODULE dynkeg