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/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90 @ 6748

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

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

Download in other formats: