1 | MODULE dynldf_lap_blp_lf |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE dynldf_lap_blp *** |
---|
4 | !! Ocean dynamics: lateral viscosity trend (laplacian and bilaplacian) |
---|
5 | !!====================================================================== |
---|
6 | !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian |
---|
7 | !! 4.0 ! 2020-04 (A. Nasser, G. Madec) Add symmetric mixing tensor |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | |
---|
10 | !!---------------------------------------------------------------------- |
---|
11 | !! dyn_ldf_lap : update the momentum trend with the lateral viscosity using an iso-level laplacian operator |
---|
12 | !! dyn_ldf_blp : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator |
---|
13 | !!---------------------------------------------------------------------- |
---|
14 | USE oce ! ocean dynamics and tracers |
---|
15 | USE dom_oce ! ocean space and time domain |
---|
16 | USE domutl, ONLY : is_tile |
---|
17 | USE ldfdyn ! lateral diffusion: eddy viscosity coef. |
---|
18 | USE ldfslp ! iso-neutral slopes |
---|
19 | USE zdf_oce ! ocean vertical physics |
---|
20 | ! |
---|
21 | USE in_out_manager ! I/O manager |
---|
22 | USE lib_mpp |
---|
23 | |
---|
24 | IMPLICIT NONE |
---|
25 | PRIVATE |
---|
26 | |
---|
27 | PUBLIC dyn_ldf_lap_lf ! called by dynldf.F90 |
---|
28 | PUBLIC dyn_ldf_blp_lf ! called by dynldf.F90 |
---|
29 | |
---|
30 | !! * Substitutions |
---|
31 | # include "do_loop_substitute.h90" |
---|
32 | # include "domzgr_substitute.h90" |
---|
33 | !!---------------------------------------------------------------------- |
---|
34 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
---|
35 | !! $Id: dynldf_lap_blp.F90 14757 2021-04-27 15:33:44Z francesca $ |
---|
36 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
37 | !!---------------------------------------------------------------------- |
---|
38 | CONTAINS |
---|
39 | |
---|
40 | SUBROUTINE dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) |
---|
41 | !! |
---|
42 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
43 | INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices |
---|
44 | INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage |
---|
45 | REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pu, pv ! before velocity [m/s] |
---|
46 | REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] |
---|
47 | !! |
---|
48 | CALL dyn_ldf_lap_lf_t( kt, Kbb, Kmm, pu, pv, is_tile(pu), pu_rhs, pv_rhs, is_tile(pu_rhs), kpass ) |
---|
49 | |
---|
50 | END SUBROUTINE dyn_ldf_lap_lf |
---|
51 | |
---|
52 | SUBROUTINE dyn_ldf_lap_lf_t( kt, Kbb, Kmm, pu, pv, ktuv, pu_rhs, pv_rhs, ktuv_rhs, kpass ) |
---|
53 | !!---------------------------------------------------------------------- |
---|
54 | !! *** ROUTINE dyn_ldf_lap *** |
---|
55 | !! |
---|
56 | !! ** Purpose : Compute the before horizontal momentum diffusive |
---|
57 | !! trend and add it to the general trend of momentum equation. |
---|
58 | !! |
---|
59 | !! ** Method : The Laplacian operator apply on horizontal velocity is |
---|
60 | !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) |
---|
61 | !! |
---|
62 | !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. |
---|
63 | !! |
---|
64 | !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ |
---|
65 | !!---------------------------------------------------------------------- |
---|
66 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
67 | INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices |
---|
68 | INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage |
---|
69 | INTEGER , INTENT(in ) :: ktuv, ktuv_rhs |
---|
70 | REAL(wp), DIMENSION(A2D_T(ktuv) ,JPK), INTENT(in ) :: pu, pv ! before velocity [m/s] |
---|
71 | REAL(wp), DIMENSION(A2D_T(ktuv_rhs),JPK), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] |
---|
72 | ! |
---|
73 | INTEGER :: ji, jj, jk ! dummy loop indices |
---|
74 | INTEGER :: iij |
---|
75 | REAL(wp) :: zsign ! local scalars |
---|
76 | REAL(wp) :: zcur, zcur_im1, zcur_jm1 ! local scalars |
---|
77 | REAL(wp) :: zdiv, zdiv_ip1, zdiv_jp1 ! local scalars |
---|
78 | REAL(wp) :: zten, zten_ip1, zten_jp1, zshe, zshe_im1, zshe_jm1 ! tension (diagonal) and shearing (anti-diagonal) terms |
---|
79 | !!---------------------------------------------------------------------- |
---|
80 | ! |
---|
81 | IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile |
---|
82 | IF( kt == nit000 .AND. lwp ) THEN |
---|
83 | WRITE(numout,*) |
---|
84 | WRITE(numout,*) 'dyn_ldf_lf : iso-level harmonic (laplacian) operator, pass=', kpass |
---|
85 | WRITE(numout,*) '~~~~~~~ ' |
---|
86 | ENDIF |
---|
87 | ENDIF |
---|
88 | ! |
---|
89 | ! Define pu_rhs/pv_rhs halo points for multi-point haloes in bilaplacian case |
---|
90 | IF( nldf_dyn == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls |
---|
91 | ELSE ; iij = 1 |
---|
92 | ENDIF |
---|
93 | ! |
---|
94 | IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign |
---|
95 | ELSE ; zsign = -1._wp ! (eddy viscosity coef. >0) |
---|
96 | ENDIF |
---|
97 | ! |
---|
98 | SELECT CASE( nn_dynldf_typ ) |
---|
99 | ! |
---|
100 | CASE ( np_typ_rot ) !== Vorticity-Divergence operator ==! |
---|
101 | ! |
---|
102 | DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) ! Horizontal slab |
---|
103 | ! ! ahm * e3 * curl (warning: computed for ji-1,jj-1) |
---|
104 | zcur = ahmf(ji ,jj ,jk) * e3f(ji ,jj ,jk) * r1_e1e2f(ji ,jj ) & ! ahmf already * by fmask |
---|
105 | & * ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & |
---|
106 | & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) |
---|
107 | zcur_jm1 = ahmf(ji ,jj-1,jk) * e3f(ji ,jj-1,jk) * r1_e1e2f(ji ,jj-1) & ! ahmf already * by fmask |
---|
108 | & * ( e2v(ji+1,jj-1) * pv(ji+1,jj-1,jk) - e2v(ji,jj-1) * pv(ji,jj-1,jk) & |
---|
109 | & - e1u(ji ,jj ) * pu(ji ,jj ,jk) + e1u(ji,jj-1) * pu(ji,jj-1,jk) ) |
---|
110 | zcur_im1 = ahmf(ji-1,jj ,jk) * e3f(ji-1,jj ,jk) * r1_e1e2f(ji-1,jj ) & ! ahmf already * by fmask |
---|
111 | & * ( e2v(ji ,jj ) * pv(ji ,jj ,jk) - e2v(ji-1,jj) * pv(ji-1,jj,jk) & |
---|
112 | & - e1u(ji-1,jj+1) * pu(ji-1,jj+1,jk) + e1u(ji-1,jj) * pu(ji-1,jj,jk) ) |
---|
113 | ! ! ahm * div (warning: computed for ji,jj) |
---|
114 | zdiv = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & ! ahmt already * by tmask |
---|
115 | & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & |
---|
116 | & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) |
---|
117 | zdiv_ip1 = ahmt(ji+1,jj,jk) * r1_e1e2t(ji+1,jj) / e3t(ji+1,jj,jk,Kbb) & ! ahmt already * by tmask |
---|
118 | & * ( e2u(ji+1,jj)*e3u(ji+1,jj,jk,Kbb) * pu(ji+1,jj,jk) - e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) & |
---|
119 | & + e1v(ji+1,jj)*e3v(ji+1,jj,jk,Kbb) * pv(ji+1,jj,jk) - e1v(ji+1,jj-1)*e3v(ji+1,jj-1,jk,Kbb) * pv(ji+1,jj-1,jk) ) |
---|
120 | zdiv_jp1 = ahmt(ji,jj+1,jk) * r1_e1e2t(ji,jj+1) / e3t(ji,jj+1,jk,Kbb) & ! ahmt already * by tmask |
---|
121 | & * ( e2u(ji,jj+1)*e3u(ji,jj+1,jk,Kbb) * pu(ji,jj+1,jk) - e2u(ji-1,jj+1)*e3u(ji-1,jj+1,jk,Kbb) * pu(ji-1,jj+1,jk) & |
---|
122 | & + e1v(ji,jj+1)*e3v(ji,jj+1,jk,Kbb) * pv(ji,jj+1,jk) - e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) ) |
---|
123 | ! ! - curl( curl) + grad( div ) |
---|
124 | pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * ( & ! * by umask is mandatory for dyn_ldf_blp use |
---|
125 | & - ( zcur - zcur_jm1 ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & |
---|
126 | & + ( zdiv_ip1 - zdiv ) * r1_e1u(ji,jj) ) |
---|
127 | ! |
---|
128 | pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * ( & ! * by vmask is mandatory for dyn_ldf_blp use |
---|
129 | & ( zcur - zcur_im1 ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & |
---|
130 | & + ( zdiv_jp1 - zdiv ) * r1_e2v(ji,jj) ) |
---|
131 | END_3D ! End of slab |
---|
132 | ! |
---|
133 | CASE ( np_typ_sym ) !== Symmetric operator ==! |
---|
134 | ! |
---|
135 | DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 ) ! Horizontal slab |
---|
136 | ! ! shearing stress component (F-point) NB : ahmf has already been multiplied by fmask |
---|
137 | zshe = ahmf(ji,jj,jk) & |
---|
138 | & * ( e1f(ji,jj) * r1_e2f(ji,jj) & |
---|
139 | & * ( pu(ji,jj+1,jk) * r1_e1u(ji,jj+1) - pu(ji,jj,jk) * r1_e1u(ji,jj) ) & |
---|
140 | & + e2f(ji,jj) * r1_e1f(ji,jj) & |
---|
141 | & * ( pv(ji+1,jj,jk) * r1_e2v(ji+1,jj) - pv(ji,jj,jk) * r1_e2v(ji,jj) ) ) |
---|
142 | zshe_im1 = ahmf(ji-1,jj,jk) & |
---|
143 | & * ( e1f(ji-1,jj) * r1_e2f(ji-1,jj) & |
---|
144 | & * ( pu(ji-1,jj+1,jk) * r1_e1u(ji-1,jj+1) - pu(ji-1,jj,jk) * r1_e1u(ji-1,jj) ) & |
---|
145 | & + e2f(ji-1,jj) * r1_e1f(ji-1,jj) & |
---|
146 | & * ( pv(ji ,jj,jk) * r1_e2v(ji ,jj) - pv(ji-1,jj,jk) * r1_e2v(ji-1,jj) ) ) |
---|
147 | zshe_jm1 = ahmf(ji,jj-1,jk) & |
---|
148 | & * ( e1f(ji,jj-1) * r1_e2f(ji,jj-1) & |
---|
149 | & * ( pu(ji,jj,jk) * r1_e1u(ji,jj) - pu(ji,jj-1,jk) * r1_e1u(ji,jj-1) ) & |
---|
150 | & + e2f(ji,jj-1) * r1_e1f(ji,jj-1) & |
---|
151 | & * ( pv(ji+1,jj-1,jk) * r1_e2v(ji+1,jj-1) - pv(ji,jj-1,jk) * r1_e2v(ji,jj-1) ) ) |
---|
152 | ! ! tension stress component (T-point) NB : ahmt has already been multiplied by tmask |
---|
153 | zten = ahmt(ji,jj,jk) & |
---|
154 | & * ( e2t(ji,jj) * r1_e1t(ji,jj) & |
---|
155 | & * ( pu(ji,jj,jk) * r1_e2u(ji,jj) - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & |
---|
156 | & - e1t(ji,jj) * r1_e2t(ji,jj) & |
---|
157 | & * ( pv(ji,jj,jk) * r1_e1v(ji,jj) - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) ) |
---|
158 | zten_ip1 = ahmt(ji+1,jj,jk) & |
---|
159 | & * ( e2t(ji+1,jj) * r1_e1t(ji+1,jj) & |
---|
160 | & * ( pu(ji+1,jj,jk) * r1_e2u(ji+1,jj) - pu(ji,jj,jk) * r1_e2u(ji,jj) ) & |
---|
161 | & - e1t(ji+1,jj) * r1_e2t(ji+1,jj) & |
---|
162 | & * ( pv(ji+1,jj,jk) * r1_e1v(ji+1,jj) - pv(ji+1,jj-1,jk) * r1_e1v(ji+1,jj-1) ) ) |
---|
163 | zten_jp1 = ahmt(ji,jj+1,jk) & |
---|
164 | & * ( e2t(ji,jj+1) * r1_e1t(ji,jj+1) & |
---|
165 | & * ( pu(ji,jj+1,jk) * r1_e2u(ji,jj+1) - pu(ji-1,jj+1,jk) * r1_e2u(ji-1,jj+1) ) & |
---|
166 | & - e1t(ji,jj+1) * r1_e2t(ji,jj+1) & |
---|
167 | & * ( pv(ji,jj+1,jk) * r1_e1v(ji,jj+1) - pv(ji,jj,jk) * r1_e1v(ji,jj) ) ) |
---|
168 | ! |
---|
169 | pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & |
---|
170 | & * ( ( zten_ip1 * e2t(ji+1,jj )*e2t(ji+1,jj ) * e3t(ji+1,jj ,jk,Kmm) & |
---|
171 | & - zten * e2t(ji ,jj )*e2t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) * r1_e2u(ji,jj) & |
---|
172 | & + ( zshe * e1f(ji ,jj )*e1f(ji ,jj ) * e3f(ji ,jj ,jk) & |
---|
173 | & - zshe_jm1 * e1f(ji ,jj-1)*e1f(ji ,jj-1) * e3f(ji ,jj-1,jk) ) * r1_e1u(ji,jj) ) |
---|
174 | ! |
---|
175 | pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) & |
---|
176 | & * ( ( zshe * e2f(ji ,jj )*e2f(ji ,jj ) * e3f(ji ,jj ,jk) & |
---|
177 | & - zshe_im1 * e2f(ji-1,jj )*e2f(ji-1,jj ) * e3f(ji-1,jj ,jk) ) * r1_e2v(ji,jj) & |
---|
178 | & - ( zten_jp1 * e1t(ji ,jj+1)*e1t(ji ,jj+1) * e3t(ji ,jj+1,jk,Kmm) & |
---|
179 | & - zten * e1t(ji ,jj )*e1t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) * r1_e1v(ji,jj) ) |
---|
180 | ! |
---|
181 | END_3D |
---|
182 | ! |
---|
183 | END SELECT |
---|
184 | ! |
---|
185 | END SUBROUTINE dyn_ldf_lap_lf_t |
---|
186 | |
---|
187 | |
---|
188 | SUBROUTINE dyn_ldf_blp_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) |
---|
189 | !!---------------------------------------------------------------------- |
---|
190 | !! *** ROUTINE dyn_ldf_blp *** |
---|
191 | !! |
---|
192 | !! ** Purpose : Compute the before lateral momentum viscous trend |
---|
193 | !! and add it to the general trend of momentum equation. |
---|
194 | !! |
---|
195 | !! ** Method : The lateral viscous trends is provided by a bilaplacian |
---|
196 | !! operator applied to before field (forward in time). |
---|
197 | !! It is computed by two successive calls to dyn_ldf_lap routine |
---|
198 | !! |
---|
199 | !! ** Action : pt(:,:,:,:,Krhs) updated with the before rotated bilaplacian diffusion |
---|
200 | !!---------------------------------------------------------------------- |
---|
201 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
---|
202 | INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices |
---|
203 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity fields |
---|
204 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! momentum trend |
---|
205 | ! |
---|
206 | REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zulap, zvlap ! laplacian at u- and v-point |
---|
207 | !!---------------------------------------------------------------------- |
---|
208 | ! |
---|
209 | IF( kt == nit000 ) THEN |
---|
210 | IF(lwp) WRITE(numout,*) |
---|
211 | IF(lwp) WRITE(numout,*) 'dyn_ldf_blp_lf : bilaplacian operator momentum ' |
---|
212 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' |
---|
213 | ENDIF |
---|
214 | ! |
---|
215 | zulap(:,:,:) = 0._wp |
---|
216 | zvlap(:,:,:) = 0._wp |
---|
217 | ! |
---|
218 | CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt (output in zlap,Kbb) |
---|
219 | ! |
---|
220 | CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) |
---|
221 | ! |
---|
222 | END SUBROUTINE dyn_ldf_blp_lf |
---|
223 | |
---|
224 | !!====================================================================== |
---|
225 | END MODULE dynldf_lap_blp_lf |
---|