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/2017/dev_r8789_sbc/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2017/dev_r8789_sbc/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 8795

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