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 branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 7099

Last change on this file since 7099 was 7099, checked in by jcastill, 7 years ago
Changes as in r7078 of the original ING branch: branches/2015/dev_r5936_INGV1_WAVE
Property svn:keywords set to `Id`
File size: 14.7 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	!! 3.6 ! 2014-10 (E. Clementi, P. Oddo) add wave contribution to surface vertical velocity
12	!!----------------------------------------------------------------------
13
14	!!----------------------------------------------------------------------
15	!! ssh_nxt : after ssh
16	!! ssh_swp : filter ans swap the ssh arrays
17	!! wzv : compute now vertical velocity
18	!!----------------------------------------------------------------------
19	USE oce ! ocean dynamics and tracers variables
20	USE dom_oce ! ocean space and time domain variables
21	USE sbc_oce ! surface boundary condition: ocean
22	USE domvvl ! Variable volume
23	USE divcur ! hor. divergence and curl (div & cur routines)
24	USE restart ! only for lrst_oce
25	USE in_out_manager ! I/O manager
26	USE prtctl ! Print control
27	USE phycst
28	USE lbclnk ! ocean lateral boundary condition (or mpp link)
29	USE lib_mpp ! MPP library
30	USE bdy_oce
31	USE bdy_par
32	USE bdydyn2d ! bdy_ssh routine
33	#if defined key_agrif
34	USE agrif_opa_update
35	USE agrif_opa_interp
36	#endif
37	#if defined key_asminc
38	USE asminc ! Assimilation increment
39	#endif
40	USE wrk_nemo ! Memory Allocation
41	USE timing ! Timing
42	USE sbcwave, ONLY: zusd2dt, zvsd2dt,wsd3d
43
44	IMPLICIT NONE
45	PRIVATE
46
47	PUBLIC ssh_nxt ! called by step.F90
48	PUBLIC wzv ! called by step.F90
49	PUBLIC ssh_swp ! called by step.F90
50
51	!! * Substitutions
52	# include "domzgr_substitute.h90"
53	# include "vectopt_loop_substitute.h90"
54	!!----------------------------------------------------------------------
55	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
56	!! $Id$
57	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
58	!!----------------------------------------------------------------------
59	CONTAINS
60
61	SUBROUTINE ssh_nxt( kt )
62	!!----------------------------------------------------------------------
63	!! * ROUTINE ssh_nxt *
64	!!
65	!! ** Purpose : compute the after ssh (ssha)
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 : ssha : after sea surface height
72	!!
73	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
74	!!----------------------------------------------------------------------
75	!
76	REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv
77	INTEGER, INTENT(in) :: kt ! time step
78	!
79	INTEGER :: jk ! dummy loop indice
80	REAL(wp) :: z2dt, z1_rau0 ! local scalars
81	!!----------------------------------------------------------------------
82	!
83	IF( nn_timing == 1 ) CALL timing_start('ssh_nxt')
84	!
85	CALL wrk_alloc( jpi, jpj, zhdiv )
86	!
87	IF( kt == nit000 ) THEN
88	!
89	IF(lwp) WRITE(numout,*)
90	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
91	IF(lwp) WRITE(numout,*) '~~~~~~~ '
92	!
93	ENDIF
94	!
95	CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity
96	!
97	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
98	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
99
100	! !------------------------------!
101	! ! After Sea Surface Height !
102	! !------------------------------!
103	zhdiv(:,:) = 0._wp
104	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
105	zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,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	z1_rau0 = 0.5_wp * r1_rau0
112	ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
113
114	#if ! defined key_dynspg_ts
115	! These lines are not necessary with time splitting since
116	! boundary condition on sea level is set during ts loop
117	#if defined key_agrif
118	CALL agrif_ssh( kt )
119	#endif
120	#if defined key_bdy
121	IF (lk_bdy) THEN
122	CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
123	CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
124	ENDIF
125	#endif
126	#endif
127
128	#if defined key_asminc
129	! ! Include the IAU weighted SSH increment
130	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
131	CALL ssh_asm_inc( kt )
132	ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
133	ENDIF
134	#endif
135
136	! !------------------------------!
137	! ! outputs !
138	! !------------------------------!
139	!
140	IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 )
141	!
142	CALL wrk_dealloc( jpi, jpj, zhdiv )
143	!
144	IF( nn_timing == 1 ) CALL timing_stop('ssh_nxt')
145	!
146	END SUBROUTINE ssh_nxt
147
148
149	SUBROUTINE wzv( kt )
150	!!----------------------------------------------------------------------
151	!! * ROUTINE wzv *
152	!!
153	!! ** Purpose : compute the now vertical velocity
154	!!
155	!! ** Method : - Using the incompressibility hypothesis, the vertical
156	!! velocity is computed by integrating the horizontal divergence
157	!! from the bottom to the surface minus the scale factor evolution.
158	!! The boundary conditions are w=0 at the bottom (no flux) and.
159	!!
160	!! ** action : wn : now vertical velocity
161	!!
162	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
163	!!----------------------------------------------------------------------
164	!
165	INTEGER, INTENT(in) :: kt ! time step
166	REAL(wp), POINTER, DIMENSION(:,: ) :: z2d
167	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv
168	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: sshnu, sshnv, dsshnu, dsshnv
169	!
170	INTEGER :: ji, jj, jk ! dummy loop indices
171	REAL(wp) :: z1_2dt ! local scalars
172	!!----------------------------------------------------------------------
173
174	IF( nn_timing == 1 ) CALL timing_start('wzv')
175	!
176	IF( kt == nit000 ) THEN
177	!
178	IF(lwp) WRITE(numout,*)
179	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
180	IF(lwp) WRITE(numout,*) '~~~~~ '
181	!
182	wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
183	!
184	ENDIF
185	! !------------------------------!
186	! ! Now Vertical Velocity !
187	! !------------------------------!
188	z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)
189	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt
190	!
191	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases
192	CALL wrk_alloc( jpi, jpj, jpk, zhdiv )
193	!
194	DO jk = 1, jpkm1
195	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
196	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
197	DO jj = 2, jpjm1
198	DO ji = fs_2, fs_jpim1 ! vector opt.
199	zhdiv(ji,jj,jk) = r1_e12t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
200	END DO
201	END DO
202	END DO
203	CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
204	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
205	! ! Same question holds for hdivn. Perhaps just for security
206	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
207	! computation of w
208	wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) &
209	& + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
210	END DO
211	! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
212	CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )
213	ELSE ! z_star and linear free surface cases
214	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
215	! computation of w
216	wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) &
217	& + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
218	END DO
219	ENDIF
220
221	!
222	! In case ln_wave and ln_sdw the surface vertical velocity is modified
223	! accounting for Sokes Drift velocity
224	IF ( ln_wave .AND. ln_sdw ) THEN
225	ALLOCATE(sshnu(jpi,jpj),sshnv(jpi,jpj),dsshnu(jpi,jpj),dsshnv(jpi,jpj))
226	sshnu (:,:) = 0._wp
227	sshnv (:,:) = 0._wp
228	dsshnu(:,:) = 0._wp
229	dsshnv(:,:) = 0._wp
230	! sshn interpolated on U-V grid
231	!---------------------------------
232	DO jj = 1 , jpjm1
233	DO ji = 1 , jpim1
234	sshnu(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * &
235	& ( sshn(ji ,jj) * tmask(ji ,jj,1) &
236	& + sshn(ji+1,jj) * tmask(ji+1,jj,1) )
237	sshnv(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * &
238	& ( sshn(ji,jj ) * tmask(ji,jji ,1) &
239	& + sshn(ji,jj+1) * tmask(ji,jj+1,1) )
240	ENDDO
241	ENDDO
242	! Compute d(ssh)/dx and d(ssh)/dy
243	!---------------------------------
244	DO jj = 1 , jpjm1
245	DO ji = 1 , jpim1
246	dsshnu(ji,jj) = ( sshnu(ji+1,jj) - sshnu(ji,jj) ) / e1u(ji,jj)
247	dsshnv(ji,jj) = ( sshnv(ji,jj+1) - sshnv(ji,jj) ) / e2v(ji,jj)
248	ENDDO
249	ENDDO
250	! Compute the surface vertical velocity accounting for the Stokes Drift
251	!---------------------------------------------------------------------
252	wn(:,:,1) = wn(:,:,1) + zusd2dt(:,:) * dsshnu(:,:) &
253	& + zvsd2dt(:,:) * dsshnv(:,:) &
254	& - ( wsd3d (:,:,1) ) * tmask(:,:,1)
255	ENDIF
256	!
257	#if defined key_bdy
258	IF (lk_bdy) THEN
259	DO jk = 1, jpkm1
260	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
261	END DO
262	ENDIF
263	#endif
264	!
265	IF( nn_timing == 1 ) CALL timing_stop('wzv')
266
267
268	END SUBROUTINE wzv
269
270	SUBROUTINE ssh_swp( kt )
271	!!----------------------------------------------------------------------
272	!! * ROUTINE ssh_nxt *
273	!!
274	!! ** Purpose : achieve the sea surface height time stepping by
275	!! applying Asselin time filter and swapping the arrays
276	!! ssha already computed in ssh_nxt
277	!!
278	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
279	!! from the filter, see Leclair and Madec 2010) and swap :
280	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
281	!! - atfp * rdt * ( emp_b - emp ) / rau0
282	!! sshn = ssha
283	!!
284	!! ** action : - sshb, sshn : before & now sea surface height
285	!! ready for the next time step
286	!!
287	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
288	!!----------------------------------------------------------------------
289	INTEGER, INTENT(in) :: kt ! ocean time-step index
290	!!----------------------------------------------------------------------
291	!
292	IF( nn_timing == 1 ) CALL timing_start('ssh_swp')
293	!
294	IF( kt == nit000 ) THEN
295	IF(lwp) WRITE(numout,*)
296	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
297	IF(lwp) WRITE(numout,*) '~~~~~~~ '
298	ENDIF
299
300	# if defined key_dynspg_ts
301	IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter
302	# else
303	IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter
304	#endif
305	sshb(:,:) = sshn(:,:) ! before <-- now
306	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
307	ELSE !** Leap-Frog time-stepping: Asselin filter + swap
308	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered
309	IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) - rnf_b(:,:) + rnf(:,:) ) * ssmask(:,:)
310	sshn(:,:) = ssha(:,:) ! now <-- after
311	ENDIF
312	!
313	! Update velocity at AGRIF zoom boundaries
314	#if defined key_agrif
315	IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt )
316	#endif
317	!
318	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
319	!
320	IF( nn_timing == 1 ) CALL timing_stop('ssh_swp')
321	!
322	END SUBROUTINE ssh_swp
323
324	!!======================================================================
325	END MODULE sshwzv

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

Download in other formats: