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/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r5936_INGV1_WAVE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 7171

Last change on this file since 7171 was 7171, checked in by emanuelaclementi, 7 years ago
#1643 Updates and bug fix on development branch 2015/dev_r5936_INGV1_WAVE
Property svn:keywords set to `Id`
File size: 13.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	!! 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 divhor ! horizontal divergence
24	USE phycst ! physical constants
25	USE bdy_oce
26	USE bdy_par
27	USE bdydyn2d ! bdy_ssh routine
28	#if defined key_agrif
29	USE agrif_opa_interp
30	#endif
31	#if defined key_asminc
32	USE asminc ! Assimilation increment
33	#endif
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 wrk_nemo ! Memory Allocation
40	USE timing ! Timing
41	USE sbcwave, ONLY: usd3dt, vsd3dt,wsd3d
42
43	IMPLICIT NONE
44	PRIVATE
45
46	PUBLIC ssh_nxt ! called by step.F90
47	PUBLIC wzv ! called by step.F90
48	PUBLIC ssh_swp ! called by step.F90
49
50	!! * Substitutions
51	# include "domzgr_substitute.h90"
52	# include "vectopt_loop_substitute.h90"
53	!!----------------------------------------------------------------------
54	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
55	!! $Id$
56	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
57	!!----------------------------------------------------------------------
58	CONTAINS
59
60	SUBROUTINE ssh_nxt( kt )
61	!!----------------------------------------------------------------------
62	!! * ROUTINE ssh_nxt *
63	!!
64	!! ** Purpose : compute the after ssh (ssha)
65	!!
66	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
67	!! is computed by integrating the horizontal divergence and multiply by
68	!! by the time step.
69	!!
70	!! ** action : ssha, after sea surface height
71	!!
72	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
73	!!----------------------------------------------------------------------
74	INTEGER, INTENT(in) :: kt ! time step
75	!
76	INTEGER :: jk ! dummy loop indice
77	REAL(wp) :: z2dt, zcoef ! local scalars
78	REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv ! 2D workspace
79	!!----------------------------------------------------------------------
80	!
81	IF( nn_timing == 1 ) CALL timing_start('ssh_nxt')
82	!
83	CALL wrk_alloc( jpi,jpj, zhdiv )
84	!
85	IF( kt == nit000 ) THEN
86	IF(lwp) WRITE(numout,*)
87	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
88	IF(lwp) WRITE(numout,*) '~~~~~~~ '
89	ENDIF
90	!
91	CALL div_hor( kt ) ! Horizontal divergence
92	!
93	z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog)
94	IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
95
96	! !------------------------------!
97	! ! After Sea Surface Height !
98	! !------------------------------!
99	zhdiv(:,:) = 0._wp
100	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
101	zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
102	END DO
103	! ! Sea surface elevation time stepping
104	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
105	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
106	!
107	zcoef = 0.5_wp * r1_rau0
108	ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
109
110	IF ( .NOT.ln_dynspg_ts ) THEN
111	! These lines are not necessary with time splitting since
112	! boundary condition on sea level is set during ts loop
113	# if defined key_agrif
114	CALL agrif_ssh( kt )
115	# endif
116	# if defined key_bdy
117	IF( lk_bdy ) THEN
118	CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary
119	CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries
120	ENDIF
121	# endif
122	ENDIF
123
124	#if defined key_asminc
125	IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN ! Include the IAU weighted SSH increment
126	CALL ssh_asm_inc( kt )
127	ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:)
128	ENDIF
129	#endif
130	! !------------------------------!
131	! ! outputs !
132	! !------------------------------!
133	!
134	IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 )
135	!
136	CALL wrk_dealloc( jpi, jpj, zhdiv )
137	!
138	IF( nn_timing == 1 ) CALL timing_stop('ssh_nxt')
139	!
140	END SUBROUTINE ssh_nxt
141
142
143	SUBROUTINE wzv( kt )
144	!!----------------------------------------------------------------------
145	!! * ROUTINE wzv *
146	!!
147	!! ** Purpose : compute the now vertical velocity
148	!!
149	!! ** Method : - Using the incompressibility hypothesis, the vertical
150	!! velocity is computed by integrating the horizontal divergence
151	!! from the bottom to the surface minus the scale factor evolution.
152	!! The boundary conditions are w=0 at the bottom (no flux) and.
153	!!
154	!! ** action : wn : now vertical velocity
155	!!
156	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
157	!!----------------------------------------------------------------------
158	INTEGER, INTENT(in) :: kt ! time step
159	!
160	INTEGER :: ji, jj, jk ! dummy loop indices
161	REAL(wp) :: z1_2dt ! local scalars
162	REAL(wp), POINTER, DIMENSION(:,: ) :: z2d
163	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv
164	REAL(wp) :: dsshnu, dsshnv
165	!!----------------------------------------------------------------------
166	!
167	IF( nn_timing == 1 ) CALL timing_start('wzv')
168	!
169	IF( kt == nit000 ) THEN
170	IF(lwp) WRITE(numout,*)
171	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
172	IF(lwp) WRITE(numout,*) '~~~~~ '
173	!
174	wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
175	ENDIF
176	! !------------------------------!
177	! ! Now Vertical Velocity !
178	! !------------------------------!
179	z1_2dt = 1. / ( 2. * rdt ) ! set time step size (Euler/Leapfrog)
180	IF( neuler == 0 .AND. kt == nit000 ) z1_2dt = 1. / rdt
181	!
182	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases
183	CALL wrk_alloc( jpi, jpj, jpk, zhdiv )
184	!
185	DO jk = 1, jpkm1
186	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
187	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
188	DO jj = 2, jpjm1
189	DO ji = fs_2, fs_jpim1 ! vector opt.
190	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) )
191	END DO
192	END DO
193	END DO
194	CALL lbc_lnk(zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
195	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
196	! ! Same question holds for hdivn. Perhaps just for security
197	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
198	! computation of w
199	wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) &
200	& + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
201	END DO
202	! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0
203	CALL wrk_dealloc( jpi, jpj, jpk, zhdiv )
204	ELSE ! z_star and linear free surface cases
205	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
206	! computation of w
207	wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) &
208	& + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) ) * tmask(:,:,jk)
209	END DO
210	ENDIF
211	!
212	! In case ln_wave and ln_sdw, the surface vertical velocity is modified
213	! accounting for Sokes Drift velocity
214	!
215	IF ( ln_wave .AND. ln_sdw ) THEN
216	! Compute d(ssh)/dx and d(ssh)/dy
217	! Compute the surface vertical velocity accounting for the Stokes Drift
218	!---------------------------------
219	DO jj = 2 , jpjm1
220	DO ji = 2 , jpim1
221	dsshnu = ( ssha(ji+1,jj) - ssha(ji,jj) ) / e1u(ji,jj)
222	dsshnv = ( ssha(ji,jj+1) - ssha(ji,jj) ) / e2v(ji,jj)
223	wn(ji,jj,1) = wn(ji,jj,1) +( usd3dt(ji,jj,1) * dsshnu &
224	& + vsd3dt(ji,jj,1) * dsshnv &
225	& - ( wsd3d (ji,jj,1) )) * tmask(ji,jj,1)
226	ENDDO
227	ENDDO
228	ENDIF
229	!
230	#if defined key_bdy
231	IF( lk_bdy ) THEN
232	DO jk = 1, jpkm1
233	wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
234	END DO
235	ENDIF
236	#endif
237	!
238	IF( nn_timing == 1 ) CALL timing_stop('wzv')
239	!
240	END SUBROUTINE wzv
241
242
243	SUBROUTINE ssh_swp( kt )
244	!!----------------------------------------------------------------------
245	!! * ROUTINE ssh_nxt *
246	!!
247	!! ** Purpose : achieve the sea surface height time stepping by
248	!! applying Asselin time filter and swapping the arrays
249	!! ssha already computed in ssh_nxt
250	!!
251	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
252	!! from the filter, see Leclair and Madec 2010) and swap :
253	!! sshn = ssha + atfp * ( sshb -2 sshn + ssha )
254	!! - atfp * rdt * ( emp_b - emp ) / rau0
255	!! sshn = ssha
256	!!
257	!! ** action : - sshb, sshn : before & now sea surface height
258	!! ready for the next time step
259	!!
260	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
261	!!----------------------------------------------------------------------
262	INTEGER, INTENT(in) :: kt ! ocean time-step index
263	!!----------------------------------------------------------------------
264	!
265	IF( nn_timing == 1 ) CALL timing_start('ssh_swp')
266	!
267	IF( kt == nit000 ) THEN
268	IF(lwp) WRITE(numout,*)
269	IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
270	IF(lwp) WRITE(numout,*) '~~~~~~~ '
271	ENDIF
272
273	IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN
274	!** Euler time-stepping: no filter
275	sshb(:,:) = sshn(:,:) ! before <-- now
276	sshn(:,:) = ssha(:,:) ! now <-- after (before already = now)
277	!
278	ELSE !** Leap-Frog time-stepping: Asselin filter + swap
279	sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered
280	IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) &
281	& - rnf_b(:,:) + rnf(:,:) &
282	& + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:)
283	sshn(:,:) = ssha(:,:) ! now <-- after
284	ENDIF
285	!
286	IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )
287	!
288	IF( nn_timing == 1 ) CALL timing_stop('ssh_swp')
289	!
290	END SUBROUTINE ssh_swp
291
292	!!======================================================================
293	END MODULE sshwzv

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

Download in other formats: