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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

sshwzv.F90 in NEMO/trunk/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/DYN/sshwzv.F90 @ 13237

Last change on this file since 13237 was 13237, checked in by smasson, 4 years ago
trunk: Mid-year merge, merge back KERNEL-06_techene_e3
Property svn:keywords set to `Id`
File size: 21.2 KB

Line
1	MODULE sshwzv
2	!!==============================================================================
3	!! * MODULE sshwzv *
4	!! Ocean dynamics : sea surface height and vertical velocity
5	!!==============================================================================
6	!! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code
7	!! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA
8	!! - ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9	!! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module
10	!! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work
11	!! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection
12	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering.
13	!!----------------------------------------------------------------------
14
15	!!----------------------------------------------------------------------
16	!! ssh_nxt : after ssh
17	!! ssh_atf : time filter the ssh arrays
18	!! wzv : compute now vertical velocity
19	!!----------------------------------------------------------------------
20	USE oce ! ocean dynamics and tracers variables
21	USE isf_oce ! ice shelf
22	USE dom_oce ! ocean space and time domain variables
23	USE sbc_oce ! surface boundary condition: ocean
24	USE domvvl ! Variable volume
25	USE divhor ! horizontal divergence
26	USE phycst ! physical constants
27	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
28	USE bdydyn2d ! bdy_ssh routine
29	#if defined key_agrif
30	USE agrif_oce_interp
31	#endif
32	!
33	USE iom
34	USE in_out_manager ! I/O manager
35	USE restart ! only for lrst_oce
36	USE prtctl ! Print control
37	USE lbclnk ! ocean lateral boundary condition (or mpp link)
38	USE lib_mpp ! MPP library
39	USE timing ! Timing
40	USE wet_dry ! Wetting/Drying flux limiting
41
42	IMPLICIT NONE
43	PRIVATE
44
45	PUBLIC ssh_nxt ! called by step.F90
46	PUBLIC wzv ! called by step.F90
47	PUBLIC wAimp ! called by step.F90
48	PUBLIC ssh_atf ! called by step.F90
49
50	!! * Substitutions
51	# include "do_loop_substitute.h90"
52	# include "domzgr_substitute.h90"
53
54	!!----------------------------------------------------------------------
55	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
56	!! $Id$
57	!! Software governed by the CeCILL license (see ./LICENSE)
58	!!----------------------------------------------------------------------
59	CONTAINS
60
61	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
62	!!----------------------------------------------------------------------
63	!! * ROUTINE ssh_nxt *
64	!!
65	!! ** Purpose : compute the after ssh (ssh(Kaa))
66	!!
67	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
68	!! is computed by integrating the horizontal divergence and multiply by
69	!! by the time step.
70	!!
71	!! ** action : ssh(:,:,Kaa), after sea surface height
72	!!
73	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
74	!!----------------------------------------------------------------------
75	INTEGER , INTENT(in ) :: kt ! time step
76	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
77	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
78	!
79	INTEGER :: jk ! dummy loop index
80	REAL(wp) :: zcoef ! local scalar
81	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
82	!!----------------------------------------------------------------------
83	!
84	IF( ln_timing ) CALL timing_start('ssh_nxt')
85	!
86	IF( kt == nit000 ) THEN
87	IF(lwp) WRITE(numout,*)
88	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
89	IF(lwp) WRITE(numout,*) '~~~~~~~ '
90	ENDIF
91	!
92	zcoef = 0.5_wp * r1_rho0
93
94	! !------------------------------!
95	! ! After Sea Surface Height !
96	! !------------------------------!
97	IF(ln_wd_il) THEN
98	CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
99	ENDIF
100
101	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
102	!
103	zhdiv(:,:) = 0._wp
104	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
105	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
106	END DO
107	! ! Sea surface elevation time stepping
108	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
109	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
110	!
111	pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
112	!
113	#if defined key_agrif
114	Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
115	CALL agrif_ssh( kt )
116	#endif
117	!
118	IF ( .NOT.ln_dynspg_ts ) THEN
119	IF( ln_bdy ) THEN
120	CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
121	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
122	ENDIF
123	ENDIF
124	! !------------------------------!
125	! ! outputs !
126	! !------------------------------!
127	!
128	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
129	!
130	IF( ln_timing ) CALL timing_stop('ssh_nxt')
131	!
132	END SUBROUTINE ssh_nxt
133
134
135	SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
136	!!----------------------------------------------------------------------
137	!! * ROUTINE wzv *
138	!!
139	!! ** Purpose : compute the now vertical velocity
140	!!
141	!! ** Method : - Using the incompressibility hypothesis, the vertical
142	!! velocity is computed by integrating the horizontal divergence
143	!! from the bottom to the surface minus the scale factor evolution.
144	!! The boundary conditions are w=0 at the bottom (no flux) and.
145	!!
146	!! ** action : pww : now vertical velocity
147	!!
148	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
149	!!----------------------------------------------------------------------
150	INTEGER , INTENT(in) :: kt ! time step
151	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
152	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
153	!
154	INTEGER :: ji, jj, jk ! dummy loop indices
155	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
156	!!----------------------------------------------------------------------
157	!
158	IF( ln_timing ) CALL timing_start('wzv')
159	!
160	IF( kt == nit000 ) THEN
161	IF(lwp) WRITE(numout,*)
162	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
163	IF(lwp) WRITE(numout,*) '~~~~~ '
164	!
165	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
166	ENDIF
167	! !------------------------------!
168	! ! Now Vertical Velocity !
169	! !------------------------------!
170	!
171	! !===============================!
172	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
173	! !===============================!
174	ALLOCATE( zhdiv(jpi,jpj,jpk) )
175	!
176	DO jk = 1, jpkm1
177	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
178	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
179	DO_2D_00_00
180	zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
181	END_2D
182	END DO
183	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
184	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
185	! ! Same question holds for hdiv. Perhaps just for security
186	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
187	! computation of w
188	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
189	& + zhdiv(:,:,jk) &
190	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
191	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
192	END DO
193	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
194	DEALLOCATE( zhdiv )
195	! !=================================!
196	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
197	! !=================================!
198	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
199	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk)
200	END DO
201	! !==========================================!
202	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
203	! !==========================================!
204	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
205	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
206	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
207	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
208	END DO
209	ENDIF
210
211	IF( ln_bdy ) THEN
212	DO jk = 1, jpkm1
213	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
214	END DO
215	ENDIF
216	!
217	#if defined key_agrif
218	IF( .NOT. AGRIF_Root() ) THEN
219	! Mask vertical velocity at first/last columns/row
220	! inside computational domain (cosmetic)
221	! --- West --- !
222	IF( lk_west) THEN
223	DO ji = mi0(2), mi1(2)
224	DO jj = 1, jpj
225	pww(ji,jj,:) = 0._wp
226	ENDDO
227	ENDDO
228	ENDIF
229	!
230	! --- East --- !
231	IF( lk_east) THEN
232	DO ji = mi0(jpiglo-1), mi1(jpiglo-1)
233	DO jj = 1, jpj
234	pww(ji,jj,:) = 0._wp
235	ENDDO
236	ENDDO
237	ENDIF
238	!
239	! --- South --- !
240	IF( lk_south) THEN
241	DO jj = mj0(2), mj1(2)
242	DO ji = 1, jpi
243	pww(ji,jj,:) = 0._wp
244	ENDDO
245	ENDDO
246	ENDIF
247	!
248	! --- North --- !
249	IF( lk_north) THEN
250	DO jj = mj0(jpjglo-1), mj1(jpjglo-1)
251	DO ji = 1, jpi
252	pww(ji,jj,:) = 0._wp
253	ENDDO
254	ENDDO
255	ENDIF
256	!
257	ENDIF
258	#endif
259	!
260	IF( ln_timing ) CALL timing_stop('wzv')
261	!
262	END SUBROUTINE wzv
263
264
265	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f )
266	!!----------------------------------------------------------------------
267	!! * ROUTINE ssh_atf *
268	!!
269	!! ** Purpose : Apply Asselin time filter to now SSH.
270	!!
271	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
272	!! from the filter, see Leclair and Madec 2010) and swap :
273	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
274	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
275	!!
276	!! ** action : - pssh(:,:,Kmm) time filtered
277	!!
278	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
279	!!----------------------------------------------------------------------
280	INTEGER , INTENT(in ) :: kt ! ocean time-step index
281	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
282	REAL(wp), DIMENSION(jpi,jpj,jpt) , TARGET, INTENT(inout) :: pssh ! SSH field
283	REAL(wp), DIMENSION(jpi,jpj ), OPTIONAL, TARGET, INTENT( out) :: pssh_f ! filtered SSH field
284	!
285	REAL(wp) :: zcoef ! local scalar
286	REAL(wp), POINTER, DIMENSION(:,:) :: zssh ! pointer for filtered SSH
287	!!----------------------------------------------------------------------
288	!
289	IF( ln_timing ) CALL timing_start('ssh_atf')
290	!
291	IF( kt == nit000 ) THEN
292	IF(lwp) WRITE(numout,*)
293	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
294	IF(lwp) WRITE(numout,*) '~~~~~~~ '
295	ENDIF
296	! !== Euler time-stepping: no filter, just swap ==!
297	IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps
298	IF( PRESENT( pssh_f ) ) THEN ; zssh => pssh_f
299	ELSE ; zssh => pssh(:,:,Kmm)
300	ENDIF
301	! ! filtered "now" field
302	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
303	IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed
304	zcoef = rn_atfp * rn_Dt * r1_rho0
305	pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) &
306	& - rnf_b(:,:) + rnf (:,:) &
307	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
308	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
309
310	! ice sheet coupling
311	IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
312
313	ENDIF
314	ENDIF
315	!
316	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm) - : ', mask1=tmask )
317	!
318	IF( ln_timing ) CALL timing_stop('ssh_atf')
319	!
320	END SUBROUTINE ssh_atf
321
322
323	SUBROUTINE wAimp( kt, Kmm )
324	!!----------------------------------------------------------------------
325	!! * ROUTINE wAimp *
326	!!
327	!! ** Purpose : compute the Courant number and partition vertical velocity
328	!! if a proportion needs to be treated implicitly
329	!!
330	!! ** Method : -
331	!!
332	!! ** action : ww : now vertical velocity (to be handled explicitly)
333	!! : wi : now vertical velocity (for implicit treatment)
334	!!
335	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
336	!! implicit scheme for vertical advection in oceanic modeling.
337	!! Ocean Modelling, 91, 38-69.
338	!!----------------------------------------------------------------------
339	INTEGER, INTENT(in) :: kt ! time step
340	INTEGER, INTENT(in) :: Kmm ! time level index
341	!
342	INTEGER :: ji, jj, jk ! dummy loop indices
343	REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars
344	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
345	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
346	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
347	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
348	!!----------------------------------------------------------------------
349	!
350	IF( ln_timing ) CALL timing_start('wAimp')
351	!
352	IF( kt == nit000 ) THEN
353	IF(lwp) WRITE(numout,*)
354	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
355	IF(lwp) WRITE(numout,*) '~~~~~ '
356	wi(:,:,:) = 0._wp
357	ENDIF
358	!
359	! Calculate Courant numbers
360	zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability)
361	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
362	DO_3D_00_00( 1, jpkm1 )
363	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
364	Cu_adv(ji,jj,jk) = zdt * &
365	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
366	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
367	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
368	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
369	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
370	& * r1_e1e2t(ji,jj) &
371	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
372	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
373	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
374	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
375	& * r1_e1e2t(ji,jj) &
376	& ) * z1_e3t
377	END_3D
378	ELSE
379	DO_3D_00_00( 1, jpkm1 )
380	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
381	Cu_adv(ji,jj,jk) = zdt * &
382	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
383	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
384	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
385	& * r1_e1e2t(ji,jj) &
386	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
387	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
388	& * r1_e1e2t(ji,jj) &
389	& ) * z1_e3t
390	END_3D
391	ENDIF
392	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
393	!
394	CALL iom_put("Courant",Cu_adv)
395	!
396	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
397	DO_3DS_11_11( jpkm1, 2, -1 )
398	!
399	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
400	! alt:
401	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
402	! zCu = Cu_adv(ji,jj,jk)
403	! ELSE
404	! zCu = Cu_adv(ji,jj,jk-1)
405	! ENDIF
406	!
407	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
408	zcff = 0._wp
409	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
410	zcff = ( zCu - Cu_min )**2
411	zcff = zcff / ( Fcu + zcff )
412	ELSE !<-- Mostly implicit
413	zcff = ( zCu - Cu_max )/ zCu
414	ENDIF
415	zcff = MIN(1._wp, zcff)
416	!
417	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
418	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
419	!
420	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
421	END_3D
422	Cu_adv(:,:,1) = 0._wp
423	ELSE
424	! Fully explicit everywhere
425	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
426	wi (:,:,:) = 0._wp
427	ENDIF
428	CALL iom_put("wimp",wi)
429	CALL iom_put("wi_cff",Cu_adv)
430	CALL iom_put("wexp",ww)
431	!
432	IF( ln_timing ) CALL timing_stop('wAimp')
433	!
434	END SUBROUTINE wAimp
435	!!======================================================================
436	END MODULE sshwzv

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

Download in other formats: