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_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DYN/sshwzv.F90 @ 14200

Last change on this file since 14200 was 14200, checked in by mcastril, 4 years ago
Merging r14117 through r14199 into dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final
Property svn:keywords set to `Id`
File size: 21.8 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	!! - ! 2020-08 (S. Techene, G. Madec) add here ssh initiatlisation
14	!!----------------------------------------------------------------------
15
16	!!----------------------------------------------------------------------
17	!! ssh_nxt : after ssh
18	!! ssh_atf : time filter the ssh arrays
19	!! wzv : compute now vertical velocity
20	!!----------------------------------------------------------------------
21	USE oce ! ocean dynamics and tracers variables
22	USE isf_oce ! ice shelf
23	USE dom_oce ! ocean space and time domain variables
24	USE sbc_oce ! surface boundary condition: ocean
25	USE domvvl ! Variable volume
26	USE divhor ! horizontal divergence
27	USE phycst ! physical constants
28	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
29	USE bdydyn2d ! bdy_ssh routine
30	USE wet_dry ! Wetting/Drying flux limiting
31	#if defined key_agrif
32	USE agrif_oce
33	USE agrif_oce_interp
34	#endif
35	!
36	USE iom
37	USE in_out_manager ! I/O manager
38	USE restart ! only for lrst_oce
39	USE prtctl ! Print control
40	USE lbclnk ! ocean lateral boundary condition (or mpp link)
41	USE lib_mpp ! MPP library
42	USE timing ! Timing
43
44	IMPLICIT NONE
45	PRIVATE
46
47	PUBLIC ssh_nxt ! called by step.F90
48	PUBLIC wzv ! called by step.F90
49	PUBLIC wAimp ! called by step.F90
50	PUBLIC ssh_atf ! called by step.F90
51
52	!! * Substitutions
53	# include "do_loop_substitute.h90"
54	# include "domzgr_substitute.h90"
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( 0, 0, 0, 0 )
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	IF( lk_west ) THEN ! --- West --- !
225	DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
226	DO jj = 1, jpj
227	pww(ji,jj,jk) = 0._wp
228	END DO
229	END DO
230	ENDIF
231	IF( lk_east ) THEN ! --- East --- !
232	DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
233	DO jj = 1, jpj
234	pww(ji,jj,jk) = 0._wp
235	END DO
236	END DO
237	ENDIF
238	IF( lk_south ) THEN ! --- South --- !
239	DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
240	DO ji = 1, jpi
241	pww(ji,jj,jk) = 0._wp
242	END DO
243	END DO
244	ENDIF
245	IF( lk_north ) THEN ! --- North --- !
246	DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
247	DO ji = 1, jpi
248	pww(ji,jj,jk) = 0._wp
249	END DO
250	END DO
251	ENDIF
252	!
253	END DO
254	!
255	ENDIF
256	#endif
257	!
258	IF( ln_timing ) CALL timing_stop('wzv')
259	!
260	END SUBROUTINE wzv
261
262
263	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f )
264	!!----------------------------------------------------------------------
265	!! * ROUTINE ssh_atf *
266	!!
267	!! ** Purpose : Apply Asselin time filter to now SSH.
268	!!
269	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
270	!! from the filter, see Leclair and Madec 2010) and swap :
271	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
272	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
273	!!
274	!! ** action : - pssh(:,:,Kmm) time filtered
275	!!
276	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
277	!!----------------------------------------------------------------------
278	INTEGER , INTENT(in ) :: kt ! ocean time-step index
279	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
280	REAL(wp), DIMENSION(jpi,jpj,jpt) , TARGET, INTENT(inout) :: pssh ! SSH field
281	REAL(wp), DIMENSION(jpi,jpj ), OPTIONAL, TARGET, INTENT( out) :: pssh_f ! filtered SSH field
282	!
283	REAL(wp) :: zcoef ! local scalar
284	REAL(wp), POINTER, DIMENSION(:,:) :: zssh ! pointer for filtered SSH
285	!!----------------------------------------------------------------------
286	!
287	IF( ln_timing ) CALL timing_start('ssh_atf')
288	!
289	IF( kt == nit000 ) THEN
290	IF(lwp) WRITE(numout,*)
291	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
292	IF(lwp) WRITE(numout,*) '~~~~~~~ '
293	ENDIF
294	! !== Euler time-stepping: no filter, just swap ==!
295	IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps
296	IF( PRESENT( pssh_f ) ) THEN ; zssh => pssh_f
297	ELSE ; zssh => pssh(:,:,Kmm)
298	ENDIF
299	! ! filtered "now" field
300	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
301	!
302	IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed
303	zcoef = rn_atfp * rn_Dt * r1_rho0
304	pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) &
305	& - rnf_b(:,:) + rnf (:,:) &
306	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
307	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
308
309	! ice sheet coupling
310	IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 ) &
311	& pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
312
313	ENDIF
314	ENDIF
315	!
316	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf - pssh(:,:,Kmm): ', mask1=tmask )
317	!
318	IF( ln_timing ) CALL timing_stop('ssh_atf')
319	!
320	END SUBROUTINE ssh_atf
321
322
323	SUBROUTINE wAimp( kt, Kmm )
324	!!----------------------------------------------------------------------
325	!! * ROUTINE wAimp *
326	!!
327	!! ** Purpose : compute the Courant number and partition vertical velocity
328	!! if a proportion needs to be treated implicitly
329	!!
330	!! ** Method : -
331	!!
332	!! ** action : ww : now vertical velocity (to be handled explicitly)
333	!! : wi : now vertical velocity (for implicit treatment)
334	!!
335	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
336	!! implicit scheme for vertical advection in oceanic modeling.
337	!! Ocean Modelling, 91, 38-69.
338	!!----------------------------------------------------------------------
339	INTEGER, INTENT(in) :: kt ! time step
340	INTEGER, INTENT(in) :: Kmm ! time level index
341	!
342	INTEGER :: ji, jj, jk ! dummy loop indices
343	REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars
344	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
345	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
346	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
347	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
348	!!----------------------------------------------------------------------
349	!
350	IF( ln_timing ) CALL timing_start('wAimp')
351	!
352	IF( kt == nit000 ) THEN
353	IF(lwp) WRITE(numout,*)
354	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
355	IF(lwp) WRITE(numout,*) '~~~~~ '
356	wi(:,:,:) = 0._wp
357	ENDIF
358	!
359	! Calculate Courant numbers
360	zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability)
361	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
362	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
363	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
364	Cu_adv(ji,jj,jk) = zdt * &
365	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
366	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
367	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
368	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
369	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
370	& * r1_e1e2t(ji,jj) &
371	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
372	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
373	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
374	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
375	& * r1_e1e2t(ji,jj) &
376	& ) * z1_e3t
377	END_3D
378	ELSE
379	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
380	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
381	Cu_adv(ji,jj,jk) = zdt * &
382	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
383	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
384	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
385	& * r1_e1e2t(ji,jj) &
386	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
387	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
388	& * r1_e1e2t(ji,jj) &
389	& ) * z1_e3t
390	END_3D
391	ENDIF
392	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
393	!
394	CALL iom_put("Courant",Cu_adv)
395	!
396	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
397	DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary
398	!
399	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
400	! alt:
401	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
402	! zCu = Cu_adv(ji,jj,jk)
403	! ELSE
404	! zCu = Cu_adv(ji,jj,jk-1)
405	! ENDIF
406	!
407	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
408	zcff = 0._wp
409	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
410	zcff = ( zCu - Cu_min )**2
411	zcff = zcff / ( Fcu + zcff )
412	ELSE !<-- Mostly implicit
413	zcff = ( zCu - Cu_max )/ zCu
414	ENDIF
415	zcff = MIN(1._wp, zcff)
416	!
417	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
418	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
419	!
420	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
421	END_3D
422	Cu_adv(:,:,1) = 0._wp
423	ELSE
424	! Fully explicit everywhere
425	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
426	wi (:,:,:) = 0._wp
427	ENDIF
428	CALL iom_put("wimp",wi)
429	CALL iom_put("wi_cff",Cu_adv)
430	CALL iom_put("wexp",ww)
431	!
432	IF( ln_timing ) CALL timing_stop('wAimp')
433	!
434	END SUBROUTINE wAimp
435
436	!!======================================================================
437	END MODULE sshwzv

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

Download in other formats: