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/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/sshwzv.F90 @ 13248

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

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

Download in other formats: