-
Notifications
You must be signed in to change notification settings - Fork 1
/
dynldf_iso.f90
201 lines (201 loc) · 10.3 KB
/
dynldf_iso.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
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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
MODULE dynldf_iso
USE oce
USE dom_oce
USE ldfdyn
USE ldftra
USE zdf_oce
USE ldfslp
USE in_out_manager
USE lib_mpp
USE lbclnk
USE prtctl
IMPLICIT NONE
PRIVATE
PUBLIC :: dyn_ldf_iso
PUBLIC :: dyn_ldf_iso_alloc
REAL(KIND = wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:, :, :) :: akzu, akzv
REAL(KIND = wp), ALLOCATABLE, SAVE, DIMENSION(:, :) :: zfuw, zdiu, zdju, zdj1u
REAL(KIND = wp), ALLOCATABLE, SAVE, DIMENSION(:, :) :: zfvw, zdiv, zdjv, zdj1v
CONTAINS
INTEGER FUNCTION dyn_ldf_iso_alloc()
ALLOCATE(akzu(jpi, jpj, jpk), zfuw(jpi, jpk), zdiu(jpi, jpk), zdju(jpi, jpk), zdj1u(jpi, jpk), akzv(jpi, jpj, jpk), zfvw(jpi, jpk), zdiv(jpi, jpk), zdjv(jpi, jpk), zdj1v(jpi, jpk), STAT = dyn_ldf_iso_alloc)
IF (dyn_ldf_iso_alloc /= 0) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
END FUNCTION dyn_ldf_iso_alloc
SUBROUTINE dyn_ldf_iso(kt)
INTEGER, INTENT( IN ) :: kt
INTEGER :: ji, jj, jk
REAL(KIND = wp) :: zabe1, zmskt, zmkt, zuav, zuwslpi, zuwslpj
REAL(KIND = wp) :: zabe2, zmskf, zmkf, zvav, zvwslpi, zvwslpj
REAL(KIND = wp) :: zcof0, zcof1, zcof2, zcof3, zcof4, zaht_0
REAL(KIND = wp), DIMENSION(jpi, jpj) :: ziut, zivf, zdku, zdk1u
REAL(KIND = wp), DIMENSION(jpi, jpj) :: zjuf, zjvt, zdkv, zdk1v
IF (kt == nit000) THEN
IF (lwp) WRITE(numout, FMT = *)
IF (lwp) WRITE(numout, FMT = *) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
IF (lwp) WRITE(numout, FMT = *) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
IF (dyn_ldf_iso_alloc() /= 0) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
END IF
IF (ln_dynldf_hor .AND. ln_traldf_iso) THEN
!$ACC KERNELS
DO jk = 1, jpk
DO jj = 2, jpjm1
DO ji = 2, jpim1
uslp(ji, jj, jk) = - (gdept_b(ji + 1, jj, jk) - gdept_b(ji, jj, jk)) * r1_e1u(ji, jj) * umask(ji, jj, jk)
vslp(ji, jj, jk) = - (gdept_b(ji, jj + 1, jk) - gdept_b(ji, jj, jk)) * r1_e2v(ji, jj) * vmask(ji, jj, jk)
wslpi(ji, jj, jk) = - (gdepw_b(ji + 1, jj, jk) - gdepw_b(ji - 1, jj, jk)) * r1_e1t(ji, jj) * tmask(ji, jj, jk) * 0.5
wslpj(ji, jj, jk) = - (gdepw_b(ji, jj + 1, jk) - gdepw_b(ji, jj - 1, jk)) * r1_e2t(ji, jj) * tmask(ji, jj, jk) * 0.5
END DO
END DO
END DO
!$ACC END KERNELS
CALL lbc_lnk_multi(uslp, 'U', - 1., vslp, 'V', - 1., wslpi, 'W', - 1., wslpj, 'W', - 1.)
END IF
zaht_0 = 0.5_wp * rn_Ud * rn_Ld
DO jk = 1, jpkm1
!$ACC KERNELS
zdk1u(:, :) = (ub(:, :, jk) - ub(:, :, jk + 1)) * umask(:, :, jk + 1)
zdk1v(:, :) = (vb(:, :, jk) - vb(:, :, jk + 1)) * vmask(:, :, jk + 1)
!$ACC END KERNELS
IF (jk == 1) THEN
!$ACC KERNELS
zdku(:, :) = zdk1u(:, :)
zdkv(:, :) = zdk1v(:, :)
!$ACC END KERNELS
ELSE
!$ACC KERNELS
zdku(:, :) = (ub(:, :, jk - 1) - ub(:, :, jk)) * umask(:, :, jk)
zdkv(:, :) = (vb(:, :, jk - 1) - vb(:, :, jk)) * vmask(:, :, jk)
!$ACC END KERNELS
END IF
IF (ln_zps) THEN
!$ACC KERNELS
DO jj = 2, jpjm1
DO ji = 2, jpi
zabe1 = (ahmt(ji, jj, jk) + rn_ahm_b) * e2t(ji, jj) * MIN(e3u_n(ji, jj, jk), e3u_n(ji - 1, jj, jk)) * r1_e1t(ji, jj)
zmskt = 1._wp / MAX(umask(ji - 1, jj, jk) + umask(ji, jj, jk + 1) + umask(ji - 1, jj, jk + 1) + umask(ji, jj, jk), 1._wp)
zcof1 = - zaht_0 * e2t(ji, jj) * zmskt * 0.5 * (uslp(ji - 1, jj, jk) + uslp(ji, jj, jk))
ziut(ji, jj) = (zabe1 * (ub(ji, jj, jk) - ub(ji - 1, jj, jk)) + zcof1 * (zdku(ji, jj) + zdk1u(ji - 1, jj) + zdk1u(ji, jj) + zdku(ji - 1, jj))) * tmask(ji, jj, jk)
END DO
END DO
!$ACC END KERNELS
ELSE
!$ACC KERNELS
DO jj = 2, jpjm1
DO ji = 2, jpi
zabe1 = (ahmt(ji, jj, jk) + rn_ahm_b) * e2t(ji, jj) * e3t_n(ji, jj, jk) * r1_e1t(ji, jj)
zmskt = 1._wp / MAX(umask(ji - 1, jj, jk) + umask(ji, jj, jk + 1) + umask(ji - 1, jj, jk + 1) + umask(ji, jj, jk), 1._wp)
zcof1 = - zaht_0 * e2t(ji, jj) * zmskt * 0.5 * (uslp(ji - 1, jj, jk) + uslp(ji, jj, jk))
ziut(ji, jj) = (zabe1 * (ub(ji, jj, jk) - ub(ji - 1, jj, jk)) + zcof1 * (zdku(ji, jj) + zdk1u(ji - 1, jj) + zdk1u(ji, jj) + zdku(ji - 1, jj))) * tmask(ji, jj, jk)
END DO
END DO
!$ACC END KERNELS
END IF
!$ACC KERNELS
DO jj = 1, jpjm1
DO ji = 1, jpim1
zabe2 = (ahmf(ji, jj, jk) + rn_ahm_b) * e1f(ji, jj) * e3f_n(ji, jj, jk) * r1_e2f(ji, jj)
zmskf = 1._wp / MAX(umask(ji, jj + 1, jk) + umask(ji, jj, jk + 1) + umask(ji, jj + 1, jk + 1) + umask(ji, jj, jk), 1._wp)
zcof2 = - zaht_0 * e1f(ji, jj) * zmskf * 0.5 * (vslp(ji + 1, jj, jk) + vslp(ji, jj, jk))
zjuf(ji, jj) = (zabe2 * (ub(ji, jj + 1, jk) - ub(ji, jj, jk)) + zcof2 * (zdku(ji, jj + 1) + zdk1u(ji, jj) + zdk1u(ji, jj + 1) + zdku(ji, jj))) * fmask(ji, jj, jk)
END DO
END DO
DO jj = 2, jpjm1
DO ji = 1, jpim1
zabe1 = (ahmf(ji, jj, jk) + rn_ahm_b) * e2f(ji, jj) * e3f_n(ji, jj, jk) * r1_e1f(ji, jj)
zmskf = 1._wp / MAX(vmask(ji + 1, jj, jk) + vmask(ji, jj, jk + 1) + vmask(ji + 1, jj, jk + 1) + vmask(ji, jj, jk), 1._wp)
zcof1 = - zaht_0 * e2f(ji, jj) * zmskf * 0.5 * (uslp(ji, jj + 1, jk) + uslp(ji, jj, jk))
zivf(ji, jj) = (zabe1 * (vb(ji + 1, jj, jk) - vb(ji, jj, jk)) + zcof1 * (zdkv(ji, jj) + zdk1v(ji + 1, jj) + zdk1v(ji, jj) + zdkv(ji + 1, jj))) * fmask(ji, jj, jk)
END DO
END DO
!$ACC END KERNELS
IF (ln_zps) THEN
!$ACC KERNELS
DO jj = 2, jpj
DO ji = 1, jpim1
zabe2 = (ahmt(ji, jj, jk) + rn_ahm_b) * e1t(ji, jj) * MIN(e3v_n(ji, jj, jk), e3v_n(ji, jj - 1, jk)) * r1_e2t(ji, jj)
zmskt = 1._wp / MAX(vmask(ji, jj - 1, jk) + vmask(ji, jj, jk + 1) + vmask(ji, jj - 1, jk + 1) + vmask(ji, jj, jk), 1._wp)
zcof2 = - zaht_0 * e1t(ji, jj) * zmskt * 0.5 * (vslp(ji, jj - 1, jk) + vslp(ji, jj, jk))
zjvt(ji, jj) = (zabe2 * (vb(ji, jj, jk) - vb(ji, jj - 1, jk)) + zcof2 * (zdkv(ji, jj - 1) + zdk1v(ji, jj) + zdk1v(ji, jj - 1) + zdkv(ji, jj))) * tmask(ji, jj, jk)
END DO
END DO
!$ACC END KERNELS
ELSE
!$ACC KERNELS
DO jj = 2, jpj
DO ji = 1, jpim1
zabe2 = (ahmt(ji, jj, jk) + rn_ahm_b) * e1t(ji, jj) * e3t_n(ji, jj, jk) * r1_e2t(ji, jj)
zmskt = 1. / MAX(vmask(ji, jj - 1, jk) + vmask(ji, jj, jk + 1) + vmask(ji, jj - 1, jk + 1) + vmask(ji, jj, jk), 1.)
zcof2 = - zaht_0 * e1t(ji, jj) * zmskt * 0.5 * (vslp(ji, jj - 1, jk) + vslp(ji, jj, jk))
zjvt(ji, jj) = (zabe2 * (vb(ji, jj, jk) - vb(ji, jj - 1, jk)) + zcof2 * (zdkv(ji, jj - 1) + zdk1v(ji, jj) + zdk1v(ji, jj - 1) + zdkv(ji, jj))) * tmask(ji, jj, jk)
END DO
END DO
!$ACC END KERNELS
END IF
!$ACC KERNELS
DO jj = 2, jpjm1
DO ji = 2, jpim1
ua(ji, jj, jk) = ua(ji, jj, jk) + (ziut(ji + 1, jj) - ziut(ji, jj) + zjuf(ji, jj) - zjuf(ji, jj - 1)) * r1_e1e2u(ji, jj) / e3u_n(ji, jj, jk)
va(ji, jj, jk) = va(ji, jj, jk) + (zivf(ji, jj) - zivf(ji - 1, jj) + zjvt(ji, jj + 1) - zjvt(ji, jj)) * r1_e1e2v(ji, jj) / e3v_n(ji, jj, jk)
END DO
END DO
!$ACC END KERNELS
END DO
IF (ln_ctl) CALL prt_ctl(tab3d_1 = ua, clinfo1 = ' ldfh - Ua: ', mask1 = umask, tab3d_2 = va, clinfo2 = ' Va: ', mask2 = vmask, clinfo3 = 'dyn')
!$ACC KERNELS
DO jj = 2, jpjm1
DO jk = 1, jpk
DO ji = 2, jpi
zdiu(ji, jk) = tmask(ji, jj, jk) * (ub(ji, jj, jk) - ub(ji - 1, jj, jk))
zdju(ji, jk) = fmask(ji, jj, jk) * (ub(ji, jj + 1, jk) - ub(ji, jj, jk))
zdjv(ji, jk) = tmask(ji, jj, jk) * (vb(ji, jj, jk) - vb(ji, jj - 1, jk))
zdj1u(ji, jk) = fmask(ji, jj - 1, jk) * (ub(ji, jj, jk) - ub(ji, jj - 1, jk))
zdj1v(ji, jk) = tmask(ji, jj + 1, jk) * (vb(ji, jj + 1, jk) - vb(ji, jj, jk))
END DO
END DO
DO jk = 1, jpk
DO ji = 1, jpim1
zdiv(ji, jk) = fmask(ji, jj, jk) * (vb(ji + 1, jj, jk) - vb(ji, jj, jk))
END DO
END DO
DO ji = 1, jpi
zfuw(ji, 1) = 0.E0
zfvw(ji, 1) = 0.E0
zfuw(ji, jpk) = 0.E0
zfvw(ji, jpk) = 0.E0
END DO
DO jk = 2, jpkm1
DO ji = 2, jpim1
zcof0 = 0.5_wp * zaht_0 * umask(ji, jj, jk)
zuwslpi = zcof0 * (wslpi(ji + 1, jj, jk) + wslpi(ji, jj, jk))
zuwslpj = zcof0 * (wslpj(ji + 1, jj, jk) + wslpj(ji, jj, jk))
zmkt = 1. / MAX(tmask(ji, jj, jk - 1) + tmask(ji + 1, jj, jk - 1) + tmask(ji, jj, jk) + tmask(ji + 1, jj, jk), 1.)
zmkf = 1. / MAX(fmask(ji, jj - 1, jk - 1) + fmask(ji, jj, jk - 1) + fmask(ji, jj - 1, jk) + fmask(ji, jj, jk), 1.)
zcof3 = - e2u(ji, jj) * zmkt * zuwslpi
zcof4 = - e1u(ji, jj) * zmkf * zuwslpj
zfuw(ji, jk) = zcof3 * (zdiu(ji, jk - 1) + zdiu(ji + 1, jk - 1) + zdiu(ji, jk) + zdiu(ji + 1, jk)) + zcof4 * (zdj1u(ji, jk - 1) + zdju(ji, jk - 1) + zdj1u(ji, jk) + zdju(ji, jk))
akzu(ji, jj, jk) = (zuwslpi * zuwslpi + zuwslpj * zuwslpj) / zaht_0
END DO
END DO
DO jk = 2, jpkm1
DO ji = 2, jpim1
zcof0 = 0.5_wp * zaht_0 * vmask(ji, jj, jk)
zvwslpi = zcof0 * (wslpi(ji, jj + 1, jk) + wslpi(ji, jj, jk))
zvwslpj = zcof0 * (wslpj(ji, jj + 1, jk) + wslpj(ji, jj, jk))
zmkf = 1. / MAX(fmask(ji - 1, jj, jk - 1) + fmask(ji, jj, jk - 1) + fmask(ji - 1, jj, jk) + fmask(ji, jj, jk), 1.)
zmkt = 1. / MAX(tmask(ji, jj, jk - 1) + tmask(ji, jj + 1, jk - 1) + tmask(ji, jj, jk) + tmask(ji, jj + 1, jk), 1.)
zcof3 = - e2v(ji, jj) * zmkf * zvwslpi
zcof4 = - e1v(ji, jj) * zmkt * zvwslpj
zfvw(ji, jk) = zcof3 * (zdiv(ji, jk - 1) + zdiv(ji - 1, jk - 1) + zdiv(ji, jk) + zdiv(ji - 1, jk)) + zcof4 * (zdjv(ji, jk - 1) + zdj1v(ji, jk - 1) + zdjv(ji, jk) + zdj1v(ji, jk))
akzv(ji, jj, jk) = (zvwslpi * zvwslpi + zvwslpj * zvwslpj) / zaht_0
END DO
END DO
DO jk = 1, jpkm1
DO ji = 2, jpim1
ua(ji, jj, jk) = ua(ji, jj, jk) + (zfuw(ji, jk) - zfuw(ji, jk + 1)) * r1_e1e2u(ji, jj) / e3u_n(ji, jj, jk)
va(ji, jj, jk) = va(ji, jj, jk) + (zfvw(ji, jk) - zfvw(ji, jk + 1)) * r1_e1e2v(ji, jj) / e3v_n(ji, jj, jk)
END DO
END DO
END DO
!$ACC END KERNELS
END SUBROUTINE dyn_ldf_iso
END MODULE dynldf_iso