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_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/DYN/sshwzv.F90 @ 13220

Last change on this file since 13220 was 13220, checked in by orioltp, 4 years ago
dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: updating from trunk r13218
Property svn:keywords set to `Id`
File size: 19.4 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	!!----------------------------------------------------------------------
14
15	!!----------------------------------------------------------------------
16	!! ssh_nxt : after ssh
17	!! ssh_atf : time filter the ssh arrays
18	!! wzv : compute now vertical velocity
19	!!----------------------------------------------------------------------
20	USE oce ! ocean dynamics and tracers variables
21	USE isf_oce ! ice shelf
22	USE dom_oce ! ocean space and time domain variables
23	USE sbc_oce ! surface boundary condition: ocean
24	USE domvvl ! Variable volume
25	USE divhor ! horizontal divergence
26	USE phycst ! physical constants
27	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
28	USE bdydyn2d ! bdy_ssh routine
29	#if defined key_agrif
30	USE agrif_oce_interp
31	#endif
32	!
33	USE iom
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 timing ! Timing
40	USE wet_dry ! Wetting/Drying flux limiting
41
42	IMPLICIT NONE
43	PRIVATE
44
45	PUBLIC ssh_nxt ! called by step.F90
46	PUBLIC wzv ! called by step.F90
47	PUBLIC wAimp ! called by step.F90
48	PUBLIC ssh_atf ! called by step.F90
49
50	!! * Substitutions
51	# include "do_loop_substitute.h90"
52	!!----------------------------------------------------------------------
53	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
54	!! $Id$
55	!! Software governed by the CeCILL license (see ./LICENSE)
56	!!----------------------------------------------------------------------
57	CONTAINS
58
59	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
60	!!----------------------------------------------------------------------
61	!! * ROUTINE ssh_nxt *
62	!!
63	!! ** Purpose : compute the after ssh (ssh(Kaa))
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 : ssh(:,:,Kaa), after sea surface height
70	!!
71	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
72	!!----------------------------------------------------------------------
73	INTEGER , INTENT(in ) :: kt ! time step
74	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
75	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
76	!
77	INTEGER :: jk ! dummy loop index
78	REAL(wp) :: zcoef ! local scalar
79	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
80	!!----------------------------------------------------------------------
81	!
82	IF( ln_timing ) CALL timing_start('ssh_nxt')
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	zcoef = 0.5_wp * r1_rho0
91
92	! !------------------------------!
93	! ! After Sea Surface Height !
94	! !------------------------------!
95	IF(ln_wd_il) THEN
96	CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
97	ENDIF
98
99	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
100	!
101	zhdiv(:,:) = 0._wp
102	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
103	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
104	END DO
105	! ! Sea surface elevation time stepping
106	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
107	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
108	!
109	pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
110	!
111	#if defined key_agrif
112	Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )
113	#endif
114	!
115	IF ( .NOT.ln_dynspg_ts ) THEN
116	IF( ln_bdy ) THEN
117	CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
118	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
119	ENDIF
120	ENDIF
121	! !------------------------------!
122	! ! outputs !
123	! !------------------------------!
124	!
125	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
126	!
127	IF( ln_timing ) CALL timing_stop('ssh_nxt')
128	!
129	END SUBROUTINE ssh_nxt
130
131
132	SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
133	!!----------------------------------------------------------------------
134	!! * ROUTINE wzv *
135	!!
136	!! ** Purpose : compute the now vertical velocity
137	!!
138	!! ** Method : - Using the incompressibility hypothesis, the vertical
139	!! velocity is computed by integrating the horizontal divergence
140	!! from the bottom to the surface minus the scale factor evolution.
141	!! The boundary conditions are w=0 at the bottom (no flux) and.
142	!!
143	!! ** action : pww : now vertical velocity
144	!!
145	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
146	!!----------------------------------------------------------------------
147	INTEGER , INTENT(in) :: kt ! time step
148	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
149	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! now vertical velocity
150	!
151	INTEGER :: ji, jj, jk ! dummy loop indices
152	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
153	!!----------------------------------------------------------------------
154	!
155	IF( ln_timing ) CALL timing_start('wzv')
156	!
157	IF( kt == nit000 ) THEN
158	IF(lwp) WRITE(numout,*)
159	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
160	IF(lwp) WRITE(numout,*) '~~~~~ '
161	!
162	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
163	ENDIF
164	! !------------------------------!
165	! ! Now Vertical Velocity !
166	! !------------------------------!
167	!
168	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases
169	ALLOCATE( zhdiv(jpi,jpj,jpk) )
170	!
171	DO jk = 1, jpkm1
172	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
173	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
174	DO_2D_00_00
175	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) )
176	END_2D
177	END DO
178	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
179	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
180	! ! Same question holds for hdiv. Perhaps just for security
181	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
182	! computation of w
183	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk) &
184	& + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
185	END DO
186	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
187	DEALLOCATE( zhdiv )
188	ELSE ! z_star and linear free surface cases
189	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
190	! computation of w
191	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
192	& + r1_Dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
193	END DO
194	ENDIF
195
196	IF( ln_bdy ) THEN
197	DO jk = 1, jpkm1
198	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
199	END DO
200	ENDIF
201	!
202	#if defined key_agrif
203	IF( .NOT. AGRIF_Root() ) THEN
204	! Mask vertical velocity at first/last columns/row
205	! inside computational domain (cosmetic)
206	! --- West --- !
207	IF( lk_west) THEN
208	DO ji = mi0(2), mi1(2)
209	DO jj = 1, jpj
210	pww(ji,jj,:) = 0._wp
211	ENDDO
212	ENDDO
213	ENDIF
214	!
215	! --- East --- !
216	IF( lk_east) THEN
217	DO ji = mi0(jpiglo-1), mi1(jpiglo-1)
218	DO jj = 1, jpj
219	pww(ji,jj,:) = 0._wp
220	ENDDO
221	ENDDO
222	ENDIF
223	!
224	! --- South --- !
225	IF( lk_south) THEN
226	DO jj = mj0(2), mj1(2)
227	DO ji = 1, jpi
228	pww(ji,jj,:) = 0._wp
229	ENDDO
230	ENDDO
231	ENDIF
232	!
233	! --- North --- !
234	IF( lk_north) THEN
235	DO jj = mj0(jpjglo-1), mj1(jpjglo-1)
236	DO ji = 1, jpi
237	pww(ji,jj,:) = 0._wp
238	ENDDO
239	ENDDO
240	ENDIF
241	!
242	ENDIF
243	#endif
244	!
245	IF( ln_timing ) CALL timing_stop('wzv')
246	!
247	END SUBROUTINE wzv
248
249
250	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
251	!!----------------------------------------------------------------------
252	!! * ROUTINE ssh_atf *
253	!!
254	!! ** Purpose : Apply Asselin time filter to now SSH.
255	!!
256	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
257	!! from the filter, see Leclair and Madec 2010) and swap :
258	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
259	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
260	!!
261	!! ** action : - pssh(:,:,Kmm) time filtered
262	!!
263	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
264	!!----------------------------------------------------------------------
265	INTEGER , INTENT(in ) :: kt ! ocean time-step index
266	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
267	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field
268	!
269	REAL(wp) :: zcoef ! local scalar
270	!!----------------------------------------------------------------------
271	!
272	IF( ln_timing ) CALL timing_start('ssh_atf')
273	!
274	IF( kt == nit000 ) THEN
275	IF(lwp) WRITE(numout,*)
276	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
277	IF(lwp) WRITE(numout,*) '~~~~~~~ '
278	ENDIF
279	! !== Euler time-stepping: no filter, just swap ==!
280	IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps
281	! ! filtered "now" field
282	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
283	IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed
284	zcoef = rn_atfp * rn_Dt * r1_rho0
285	pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) &
286	& - rnf_b(:,:) + rnf (:,:) &
287	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
288	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
289
290	! ice sheet coupling
291	IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
292
293	ENDIF
294	ENDIF
295	!
296	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm) - : ', mask1=tmask )
297	!
298	IF( ln_timing ) CALL timing_stop('ssh_atf')
299	!
300	END SUBROUTINE ssh_atf
301
302	SUBROUTINE wAimp( kt, Kmm )
303	!!----------------------------------------------------------------------
304	!! * ROUTINE wAimp *
305	!!
306	!! ** Purpose : compute the Courant number and partition vertical velocity
307	!! if a proportion needs to be treated implicitly
308	!!
309	!! ** Method : -
310	!!
311	!! ** action : ww : now vertical velocity (to be handled explicitly)
312	!! : wi : now vertical velocity (for implicit treatment)
313	!!
314	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
315	!! implicit scheme for vertical advection in oceanic modeling.
316	!! Ocean Modelling, 91, 38-69.
317	!!----------------------------------------------------------------------
318	INTEGER, INTENT(in) :: kt ! time step
319	INTEGER, INTENT(in) :: Kmm ! time level index
320	!
321	INTEGER :: ji, jj, jk ! dummy loop indices
322	REAL(wp) :: zCu, zcff, z1_e3t ! local scalars
323	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
324	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
325	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
326	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
327	!!----------------------------------------------------------------------
328	!
329	IF( ln_timing ) CALL timing_start('wAimp')
330	!
331	IF( kt == nit000 ) THEN
332	IF(lwp) WRITE(numout,*)
333	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
334	IF(lwp) WRITE(numout,*) '~~~~~ '
335	wi(:,:,:) = 0._wp
336	ENDIF
337	!
338	! Calculate Courant numbers
339	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
340	DO_3D_00_00( 1, jpkm1 )
341	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
342	! 2*rn_Dt and not rDt (for restartability)
343	Cu_adv(ji,jj,jk) = 2._wp * rn_Dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
344	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
345	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
346	& * r1_e1e2t(ji,jj) &
347	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
348	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
349	& * r1_e1e2t(ji,jj) &
350	& ) * z1_e3t
351	END_3D
352	ELSE
353	DO_3D_00_00( 1, jpkm1 )
354	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
355	! 2*rn_Dt and not rDt (for restartability)
356	Cu_adv(ji,jj,jk) = 2._wp * rn_Dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
357	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
358	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
359	& * r1_e1e2t(ji,jj) &
360	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
361	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
362	& * r1_e1e2t(ji,jj) &
363	& ) * z1_e3t
364	END_3D
365	ENDIF
366	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
367	!
368	CALL iom_put("Courant",Cu_adv)
369	!
370	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
371	DO_3DS_11_11( jpkm1, 2, -1 )
372	!
373	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
374	! alt:
375	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
376	! zCu = Cu_adv(ji,jj,jk)
377	! ELSE
378	! zCu = Cu_adv(ji,jj,jk-1)
379	! ENDIF
380	!
381	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
382	zcff = 0._wp
383	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
384	zcff = ( zCu - Cu_min )**2
385	zcff = zcff / ( Fcu + zcff )
386	ELSE !<-- Mostly implicit
387	zcff = ( zCu - Cu_max )/ zCu
388	ENDIF
389	zcff = MIN(1._wp, zcff)
390	!
391	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
392	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
393	!
394	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
395	END_3D
396	Cu_adv(:,:,1) = 0._wp
397	ELSE
398	! Fully explicit everywhere
399	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
400	wi (:,:,:) = 0._wp
401	ENDIF
402	CALL iom_put("wimp",wi)
403	CALL iom_put("wi_cff",Cu_adv)
404	CALL iom_put("wexp",ww)
405	!
406	IF( ln_timing ) CALL timing_stop('wAimp')
407	!
408	END SUBROUTINE wAimp
409	!!======================================================================
410	END MODULE sshwzv

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

Download in other formats: