Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynzdf_imp.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago
The Dr Hook changes from my perl code.
File size: 17.8 KB

Line
1	MODULE dynzdf_imp
2	!!======================================================================
3	!! * MODULE dynzdf_imp *
4	!! Ocean dynamics: vertical component(s) of the momentum mixing trend
5	!!======================================================================
6	!! History : OPA ! 1990-10 (B. Blanke) Original code
7	!! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal
8	!! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module
9	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps
10	!! 3.4 ! 2012-01 (H. Liu) Semi-implicit bottom friction
11	!!----------------------------------------------------------------------
12
13	!!----------------------------------------------------------------------
14	!! dyn_zdf_imp : update the momentum trend with the vertical diffusion using a implicit time-stepping
15	!!----------------------------------------------------------------------
16	USE oce ! ocean dynamics and tracers
17	USE dom_oce ! ocean space and time domain
18	USE domvvl ! variable volume
19	USE sbc_oce ! surface boundary condition: ocean
20	USE zdf_oce ! ocean vertical physics
21	USE phycst ! physical constants
22	USE in_out_manager ! I/O manager
23	USE lib_mpp ! MPP library
24	USE zdfbfr ! Bottom friction setup
25	USE wrk_nemo ! Memory Allocation
26	USE timing ! Timing
27	USE dynadv ! dynamics: vector invariant versus flux form
28	USE dynspg_oce, ONLY: lk_dynspg_ts
29
30	USE yomhook, ONLY: lhook, dr_hook
31	USE parkind1, ONLY: jprb, jpim
32
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC dyn_zdf_imp ! called by step.F90
37
38	REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise
39
40	!! * Substitutions
41	# include "domzgr_substitute.h90"
42	# include "vectopt_loop_substitute.h90"
43	!!----------------------------------------------------------------------
44	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
45	!! $Id$
46	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
47	!!----------------------------------------------------------------------
48	CONTAINS
49
50	SUBROUTINE dyn_zdf_imp( kt, p2dt )
51	!!----------------------------------------------------------------------
52	!! * ROUTINE dyn_zdf_imp *
53	!!
54	!! ** Purpose : Compute the trend due to the vert. momentum diffusion
55	!! and the surface forcing, and add it to the general trend of
56	!! the momentum equations.
57	!!
58	!! ** Method : The vertical momentum mixing trend is given by :
59	!! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
60	!! backward time stepping
61	!! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2)
62	!! Bottom boundary conditions : bottom stress (cf zdfbfr.F)
63	!! Add this trend to the general trend ua :
64	!! ua = ua + dz( avmu dz(u) )
65	!!
66	!! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.
67	!!---------------------------------------------------------------------
68	INTEGER , INTENT(in) :: kt ! ocean time-step index
69	REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step
70	!!
71	INTEGER :: ji, jj, jk ! dummy loop indices
72	INTEGER :: ikbu, ikbv ! local integers
73	REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars
74	REAL(wp) :: ze3ua, ze3va
75	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws
76	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
77	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
78	REAL(KIND=jprb) :: zhook_handle
79
80	CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_ZDF_IMP'
81
82	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
83
84	!!----------------------------------------------------------------------
85	!
86	IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_imp')
87	!
88	CALL wrk_alloc( jpi,jpj,jpk, zwi, zwd, zws )
89	!
90	IF( kt == nit000 ) THEN
91	IF(lwp) WRITE(numout,*)
92	IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
93	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
94	!
95	IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator
96	ELSE ; r_vvl = 0._wp
97	ENDIF
98	ENDIF
99
100	! 0. Local constant initialization
101	! --------------------------------
102	z1_p2dt = 1._wp / p2dt ! inverse of the timestep
103
104	! 1. Apply semi-implicit bottom friction
105	! --------------------------------------
106	! Only needed for semi-implicit bottom friction setup. The explicit
107	! bottom friction has been included in "u(v)a" which act as the R.H.S
108	! column vector of the tri-diagonal matrix equation
109	!
110
111	IF( ln_bfrimp ) THEN
112	DO jj = 2, jpjm1
113	DO ji = 2, jpim1
114	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
115	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
116	avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)
117	avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)
118	END DO
119	END DO
120	IF ( ln_isfcav ) THEN
121	DO jj = 2, jpjm1
122	DO ji = 2, jpim1
123	ikbu = miku(ji,jj) ! ocean top level at u- and v-points
124	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
125	IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)
126	IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)
127	END DO
128	END DO
129	END IF
130	ENDIF
131
132	#if defined key_dynspg_ts
133	IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity
134	DO jk = 1, jpkm1
135	ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk)
136	va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk)
137	END DO
138	ELSE ! applied on thickness weighted velocity
139	DO jk = 1, jpkm1
140	ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) &
141	& + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) &
142	& / fse3u_a(:,:,jk) * umask(:,:,jk)
143	va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) &
144	& + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) &
145	& / fse3v_a(:,:,jk) * vmask(:,:,jk)
146	END DO
147	ENDIF
148
149	IF ( ln_bfrimp .AND.lk_dynspg_ts ) THEN
150	! remove barotropic velocities:
151	DO jk = 1, jpkm1
152	ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk)
153	va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk)
154	END DO
155	! Add bottom/top stress due to barotropic component only:
156	DO jj = 2, jpjm1
157	DO ji = fs_2, fs_jpim1 ! vector opt.
158	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
159	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
160	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)
161	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)
162	ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua
163	va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va
164	END DO
165	END DO
166	IF ( ln_isfcav ) THEN
167	DO jj = 2, jpjm1
168	DO ji = fs_2, fs_jpim1 ! vector opt.
169	ikbu = miku(ji,jj) ! top ocean level at u- and v-points
170	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
171	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)
172	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)
173	ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua
174	va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va
175	END DO
176	END DO
177	END IF
178	ENDIF
179	#endif
180
181	! 2. Vertical diffusion on u
182	! ---------------------------
183	! Matrix and second member construction
184	! bottom boundary condition: both zwi and zws must be masked as avmu can take
185	! non zero value at the ocean bottom depending on the bottom friction used.
186	!
187	DO jk = 1, jpkm1 ! Matrix
188	DO jj = 2, jpjm1
189	DO ji = fs_2, fs_jpim1 ! vector opt.
190	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point
191	zcoef = - p2dt / ze3ua
192	zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )
193	zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk )
194	zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
195	zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1)
196	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
197	END DO
198	END DO
199	END DO
200	DO jj = 2, jpjm1 ! Surface boundary conditions
201	DO ji = fs_2, fs_jpim1 ! vector opt.
202	zwi(ji,jj,1) = 0._wp
203	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
204	END DO
205	END DO
206
207	! Matrix inversion starting from the first level
208	!-----------------------------------------------------------------------
209	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
210	!
211	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
212	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
213	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
214	! ( ... )( ... ) ( ... )
215	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
216	!
217	! m is decomposed in the product of an upper and a lower triangular matrix
218	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
219	! The solution (the after velocity) is in ua
220	!-----------------------------------------------------------------------
221	!
222	!== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
223	DO jk = 2, jpkm1
224	DO jj = 2, jpjm1
225	DO ji = fs_2, fs_jpim1 ! vector opt.
226	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
227	END DO
228	END DO
229	END DO
230	!
231	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
232	DO ji = fs_2, fs_jpim1 ! vector opt.
233	#if defined key_dynspg_ts
234	ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1)
235	ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
236	& / ( ze3ua * rau0 ) * umask(ji,jj,1)
237	#else
238	ua(ji,jj,1) = ub(ji,jj,1) &
239	& + p2dt (ua(ji,jj,1) + 0.5_wp ( utau_b(ji,jj) + utau(ji,jj) ) &
240	& / ( fse3u(ji,jj,1) * rau0 ) * umask(ji,jj,1) )
241	#endif
242	END DO
243	END DO
244	DO jk = 2, jpkm1
245	DO jj = 2, jpjm1
246	DO ji = fs_2, fs_jpim1
247	#if defined key_dynspg_ts
248	zrhs = ua(ji,jj,jk) ! zrhs=right hand side
249	#else
250	zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)
251	#endif
252	ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)
253	END DO
254	END DO
255	END DO
256	!
257	DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==
258	DO ji = fs_2, fs_jpim1 ! vector opt.
259	ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
260	END DO
261	END DO
262	DO jk = jpk-2, 1, -1
263	DO jj = 2, jpjm1
264	DO ji = fs_2, fs_jpim1
265	ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)
266	END DO
267	END DO
268	END DO
269
270	#if ! defined key_dynspg_ts
271	! Normalization to obtain the general momentum trend ua
272	DO jk = 1, jpkm1
273	DO jj = 2, jpjm1
274	DO ji = fs_2, fs_jpim1 ! vector opt.
275	ua(ji,jj,jk) = ( ua(ji,jj,jk) - ub(ji,jj,jk) ) * z1_p2dt
276	END DO
277	END DO
278	END DO
279	#endif
280
281	! 3. Vertical diffusion on v
282	! ---------------------------
283	! Matrix and second member construction
284	! bottom boundary condition: both zwi and zws must be masked as avmv can take
285	! non zero value at the ocean bottom depending on the bottom friction used
286	!
287	DO jk = 1, jpkm1 ! Matrix
288	DO jj = 2, jpjm1
289	DO ji = fs_2, fs_jpim1 ! vector opt.
290	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk) + r_vvl * fse3v_a(ji,jj,jk) ! after scale factor at T-point
291	zcoef = - p2dt / ze3va
292	zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )
293	zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk)
294	zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
295	zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1)
296	zwd(ji,jj,jk) = 1._wp - zzwi - zzws
297	END DO
298	END DO
299	END DO
300	DO jj = 2, jpjm1 ! Surface boundary conditions
301	DO ji = fs_2, fs_jpim1 ! vector opt.
302	zwi(ji,jj,1) = 0._wp
303	zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
304	END DO
305	END DO
306
307	! Matrix inversion
308	!-----------------------------------------------------------------------
309	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
310	!
311	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
312	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
313	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
314	! ( ... )( ... ) ( ... )
315	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
316	!
317	! m is decomposed in the product of an upper and lower triangular matrix
318	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
319	! The solution (after velocity) is in 2d array va
320	!-----------------------------------------------------------------------
321	!
322	!== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
323	DO jk = 2, jpkm1
324	DO jj = 2, jpjm1
325	DO ji = fs_2, fs_jpim1 ! vector opt.
326	zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
327	END DO
328	END DO
329	END DO
330	!
331	DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==
332	DO ji = fs_2, fs_jpim1 ! vector opt.
333	#if defined key_dynspg_ts
334	ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1)
335	va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
336	& / ( ze3va * rau0 ) * vmask(ji,jj,1)
337	#else
338	va(ji,jj,1) = vb(ji,jj,1) &
339	& + p2dt (va(ji,jj,1) + 0.5_wp ( vtau_b(ji,jj) + vtau(ji,jj) ) &
340	& / ( fse3v(ji,jj,1) * rau0 ) * vmask(ji,jj,1) )
341	#endif
342	END DO
343	END DO
344	DO jk = 2, jpkm1
345	DO jj = 2, jpjm1
346	DO ji = fs_2, fs_jpim1 ! vector opt.
347	#if defined key_dynspg_ts
348	zrhs = va(ji,jj,jk) ! zrhs=right hand side
349	#else
350	zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)
351	#endif
352	va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)
353	END DO
354	END DO
355	END DO
356	!
357	DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==
358	DO ji = fs_2, fs_jpim1 ! vector opt.
359	va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
360	END DO
361	END DO
362	DO jk = jpk-2, 1, -1
363	DO jj = 2, jpjm1
364	DO ji = fs_2, fs_jpim1
365	va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)
366	END DO
367	END DO
368	END DO
369
370	! Normalization to obtain the general momentum trend va
371	#if ! defined key_dynspg_ts
372	DO jk = 1, jpkm1
373	DO jj = 2, jpjm1
374	DO ji = fs_2, fs_jpim1 ! vector opt.
375	va(ji,jj,jk) = ( va(ji,jj,jk) - vb(ji,jj,jk) ) * z1_p2dt
376	END DO
377	END DO
378	END DO
379	#endif
380
381	! J. Chanut: Lines below are useless ?
382	!! restore bottom layer avmu(v)
383	IF( ln_bfrimp ) THEN
384	DO jj = 2, jpjm1
385	DO ji = 2, jpim1
386	ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points
387	ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
388	avmu(ji,jj,ikbu+1) = 0.e0
389	avmv(ji,jj,ikbv+1) = 0.e0
390	END DO
391	END DO
392	IF (ln_isfcav) THEN
393	DO jj = 2, jpjm1
394	DO ji = 2, jpim1
395	ikbu = miku(ji,jj) ! ocean top level at u- and v-points
396	ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points)
397	IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0
398	IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0
399	END DO
400	END DO
401	END IF
402	ENDIF
403	!
404	CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)
405	!
406	IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')
407	!
408	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
409	END SUBROUTINE dyn_zdf_imp
410
411	!!==============================================================================
412	END MODULE dynzdf_imp

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90 @ 11738

Download in other formats: