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_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DYN/sshwzv.F90 @ 13998

Last change on this file since 13998 was 13998, checked in by techene, 3 years ago
branch updated with trunk 13787
Property svn:keywords set to `Id`
File size: 25.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	!! 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	!! ssh_init_rst : ssh set from restart or domcfg.nc file or usr_def_istat_ssh
21	!!----------------------------------------------------------------------
22	USE oce ! ocean dynamics and tracers variables
23	USE isf_oce ! ice shelf
24	USE dom_oce ! ocean space and time domain variables
25	USE sbc_oce ! surface boundary condition: ocean
26	USE domvvl ! Variable volume
27	USE divhor ! horizontal divergence
28	USE phycst ! physical constants
29	USE bdy_oce , ONLY : ln_bdy, bdytmask ! Open BounDarY
30	USE bdydyn2d ! bdy_ssh routine
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	USE wet_dry ! Wetting/Drying flux limiting
44	USE usrdef_istate, ONLY : usr_def_istate_ssh ! user defined ssh initial state
45
46	IMPLICIT NONE
47	PRIVATE
48
49	PUBLIC ssh_nxt ! called by step.F90
50	PUBLIC wzv ! called by step.F90
51	PUBLIC wAimp ! called by step.F90
52	PUBLIC ssh_atf ! called by step.F90
53	PUBLIC ssh_init_rst ! called by domain.F90
54
55	!! * Substitutions
56	# include "do_loop_substitute.h90"
57	# include "domzgr_substitute.h90"
58	!!----------------------------------------------------------------------
59	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
60	!! $Id$
61	!! Software governed by the CeCILL license (see ./LICENSE)
62	!!----------------------------------------------------------------------
63	CONTAINS
64
65	SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
66	!!----------------------------------------------------------------------
67	!! * ROUTINE ssh_nxt *
68	!!
69	!! ** Purpose : compute the after ssh (ssh(Kaa))
70	!!
71	!! ** Method : - Using the incompressibility hypothesis, the ssh increment
72	!! is computed by integrating the horizontal divergence and multiply by
73	!! by the time step.
74	!!
75	!! ** action : ssh(:,:,Kaa), after sea surface height
76	!!
77	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
78	!!----------------------------------------------------------------------
79	INTEGER , INTENT(in ) :: kt ! time step
80	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
81	REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
82	!
83	INTEGER :: jk ! dummy loop index
84	REAL(wp) :: zcoef ! local scalar
85	REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
86	!!----------------------------------------------------------------------
87	!
88	IF( ln_timing ) CALL timing_start('ssh_nxt')
89	!
90	IF( kt == nit000 ) THEN
91	IF(lwp) WRITE(numout,*)
92	IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
93	IF(lwp) WRITE(numout,*) '~~~~~~~ '
94	ENDIF
95	!
96	zcoef = 0.5_wp * r1_rho0
97
98	! !------------------------------!
99	! ! After Sea Surface Height !
100	! !------------------------------!
101	IF(ln_wd_il) THEN
102	CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
103	ENDIF
104
105	CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
106	!
107	zhdiv(:,:) = 0._wp
108	DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports
109	zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
110	END DO
111	! ! Sea surface elevation time stepping
112	! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
113	! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
114	!
115	pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - rDt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)
116	!
117	#if defined key_agrif
118	Kbb_a = Kbb ; Kmm_a = Kmm ; Krhs_a = Kaa
119	CALL agrif_ssh( kt )
120	#endif
121	!
122	IF ( .NOT.ln_dynspg_ts ) THEN
123	IF( ln_bdy ) THEN
124	CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
125	CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
126	ENDIF
127	ENDIF
128	! !------------------------------!
129	! ! outputs !
130	! !------------------------------!
131	!
132	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
133	!
134	IF( ln_timing ) CALL timing_stop('ssh_nxt')
135	!
136	END SUBROUTINE ssh_nxt
137
138
139	SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
140	!!----------------------------------------------------------------------
141	!! * ROUTINE wzv *
142	!!
143	!! ** Purpose : compute the now vertical velocity
144	!!
145	!! ** Method : - Using the incompressibility hypothesis, the vertical
146	!! velocity is computed by integrating the horizontal divergence
147	!! from the bottom to the surface minus the scale factor evolution.
148	!! The boundary conditions are w=0 at the bottom (no flux) and.
149	!!
150	!! ** action : pww : now vertical velocity
151	!!
152	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
153	!!----------------------------------------------------------------------
154	INTEGER , INTENT(in) :: kt ! time step
155	INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices
156	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
157	!
158	INTEGER :: ji, jj, jk ! dummy loop indices
159	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
160	!!----------------------------------------------------------------------
161	!
162	IF( ln_timing ) CALL timing_start('wzv')
163	!
164	IF( kt == nit000 ) THEN
165	IF(lwp) WRITE(numout,*)
166	IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
167	IF(lwp) WRITE(numout,*) '~~~~~ '
168	!
169	pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
170	ENDIF
171	! !------------------------------!
172	! ! Now Vertical Velocity !
173	! !------------------------------!
174	!
175	! !===============================!
176	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
177	! !===============================!
178	ALLOCATE( zhdiv(jpi,jpj,jpk) )
179	!
180	DO jk = 1, jpkm1
181	! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
182	! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
183	DO_2D( 0, 0, 0, 0 )
184	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) )
185	END_2D
186	END DO
187	CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
188	! ! Is it problematic to have a wrong vertical velocity in boundary cells?
189	! ! Same question holds for hdiv. Perhaps just for security
190	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
191	! computation of w
192	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
193	& + zhdiv(:,:,jk) &
194	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
195	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
196	END DO
197	! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
198	DEALLOCATE( zhdiv )
199	! !=================================!
200	ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
201	! !=================================!
202	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
203	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) ) * tmask(:,:,jk)
204	END DO
205	! !==========================================!
206	ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
207	! !==========================================!
208	DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
209	pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) &
210	& + r1_Dt * ( e3t(:,:,jk,Kaa) &
211	& - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
212	END DO
213	ENDIF
214
215	IF( ln_bdy ) THEN
216	DO jk = 1, jpkm1
217	pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
218	END DO
219	ENDIF
220	!
221	#if defined key_agrif
222	IF( .NOT. AGRIF_Root() ) THEN
223	!
224	! Mask vertical velocity at first/last columns/row
225	! inside computational domain (cosmetic)
226	DO jk = 1, jpkm1
227	IF( lk_west ) THEN ! --- West --- !
228	DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
229	DO jj = 1, jpj
230	pww(ji,jj,jk) = 0._wp
231	END DO
232	END DO
233	ENDIF
234	IF( lk_east ) THEN ! --- East --- !
235	DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
236	DO jj = 1, jpj
237	pww(ji,jj,jk) = 0._wp
238	END DO
239	END DO
240	ENDIF
241	IF( lk_south ) THEN ! --- South --- !
242	DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
243	DO ji = 1, jpi
244	pww(ji,jj,jk) = 0._wp
245	END DO
246	END DO
247	ENDIF
248	IF( lk_north ) THEN ! --- North --- !
249	DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
250	DO ji = 1, jpi
251	pww(ji,jj,jk) = 0._wp
252	END DO
253	END DO
254	ENDIF
255	!
256	END DO
257	!
258	ENDIF
259	#endif
260	!
261	IF( ln_timing ) CALL timing_stop('wzv')
262	!
263	END SUBROUTINE wzv
264
265
266	SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh, pssh_f )
267	!!----------------------------------------------------------------------
268	!! * ROUTINE ssh_atf *
269	!!
270	!! ** Purpose : Apply Asselin time filter to now SSH.
271	!!
272	!! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing
273	!! from the filter, see Leclair and Madec 2010) and swap :
274	!! pssh(:,:,Kmm) = pssh(:,:,Kaa) + rn_atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
275	!! - rn_atfp * rn_Dt * ( emp_b - emp ) / rho0
276	!!
277	!! ** action : - pssh(:,:,Kmm) time filtered
278	!!
279	!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
280	!!----------------------------------------------------------------------
281	INTEGER , INTENT(in ) :: kt ! ocean time-step index
282	INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices
283	REAL(wp), DIMENSION(jpi,jpj,jpt) , TARGET, INTENT(inout) :: pssh ! SSH field
284	REAL(wp), DIMENSION(jpi,jpj ), OPTIONAL, TARGET, INTENT( out) :: pssh_f ! filtered SSH field
285	!
286	REAL(wp) :: zcoef ! local scalar
287	REAL(wp), POINTER, DIMENSION(:,:) :: zssh ! pointer for filtered SSH
288	!!----------------------------------------------------------------------
289	!
290	IF( ln_timing ) CALL timing_start('ssh_atf')
291	!
292	IF( kt == nit000 ) THEN
293	IF(lwp) WRITE(numout,*)
294	IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
295	IF(lwp) WRITE(numout,*) '~~~~~~~ '
296	ENDIF
297	! !== Euler time-stepping: no filter, just swap ==!
298	IF ( .NOT.( l_1st_euler ) ) THEN ! Only do time filtering for leapfrog timesteps
299	IF( PRESENT( pssh_f ) ) THEN ; zssh => pssh_f
300	ELSE ; zssh => pssh(:,:,Kmm)
301	ENDIF
302	! ! filtered "now" field
303	pssh(:,:,Kmm) = pssh(:,:,Kmm) + rn_atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
304	!
305	IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed
306	zcoef = rn_atfp * rn_Dt * r1_rho0
307	pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) &
308	& - rnf_b(:,:) + rnf (:,:) &
309	& + fwfisf_cav_b(:,:) - fwfisf_cav(:,:) &
310	& + fwfisf_par_b(:,:) - fwfisf_par(:,:) ) * ssmask(:,:)
311
312	! ice sheet coupling
313	IF( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1 ) &
314	& pssh(:,:,Kbb) = pssh(:,:,Kbb) - rn_atfp * rn_Dt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
315
316	ENDIF
317	ENDIF
318	!
319	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' atf - pssh(:,:,Kmm): ', mask1=tmask )
320	!
321	IF( ln_timing ) CALL timing_stop('ssh_atf')
322	!
323	END SUBROUTINE ssh_atf
324
325
326	SUBROUTINE wAimp( kt, Kmm )
327	!!----------------------------------------------------------------------
328	!! * ROUTINE wAimp *
329	!!
330	!! ** Purpose : compute the Courant number and partition vertical velocity
331	!! if a proportion needs to be treated implicitly
332	!!
333	!! ** Method : -
334	!!
335	!! ** action : ww : now vertical velocity (to be handled explicitly)
336	!! : wi : now vertical velocity (for implicit treatment)
337	!!
338	!! Reference : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
339	!! implicit scheme for vertical advection in oceanic modeling.
340	!! Ocean Modelling, 91, 38-69.
341	!!----------------------------------------------------------------------
342	INTEGER, INTENT(in) :: kt ! time step
343	INTEGER, INTENT(in) :: Kmm ! time level index
344	!
345	INTEGER :: ji, jj, jk ! dummy loop indices
346	REAL(wp) :: zCu, zcff, z1_e3t, zdt ! local scalars
347	REAL(wp) , PARAMETER :: Cu_min = 0.15_wp ! local parameters
348	REAL(wp) , PARAMETER :: Cu_max = 0.30_wp ! local parameters
349	REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters
350	REAL(wp) , PARAMETER :: Fcu = 4._wpCu_max(Cu_max-Cu_min) ! local parameters
351	!!----------------------------------------------------------------------
352	!
353	IF( ln_timing ) CALL timing_start('wAimp')
354	!
355	IF( kt == nit000 ) THEN
356	IF(lwp) WRITE(numout,*)
357	IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
358	IF(lwp) WRITE(numout,*) '~~~~~ '
359	wi(:,:,:) = 0._wp
360	ENDIF
361	!
362	! Calculate Courant numbers
363	zdt = 2._wp * rn_Dt ! 2*rn_Dt and not rDt (for restartability)
364	IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
365	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
366	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
367	Cu_adv(ji,jj,jk) = zdt * &
368	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
369	& + ( MAX( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) &
370	& * uu (ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - &
371	& MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) &
372	& * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) &
373	& * r1_e1e2t(ji,jj) &
374	& + ( MAX( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) &
375	& * vv (ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - &
376	& MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) &
377	& * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) &
378	& * r1_e1e2t(ji,jj) &
379	& ) * z1_e3t
380	END_3D
381	ELSE
382	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
383	z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
384	Cu_adv(ji,jj,jk) = zdt * &
385	& ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) &
386	& + ( MAX( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm)uu(ji ,jj,jk,Kmm), 0._wp ) - &
387	& MIN( e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm)uu(ji-1,jj,jk,Kmm), 0._wp ) ) &
388	& * r1_e1e2t(ji,jj) &
389	& + ( MAX( e1v(ji,jj )e3v(ji,jj ,jk,Kmm)vv(ji,jj ,jk,Kmm), 0._wp ) - &
390	& MIN( e1v(ji,jj-1)e3v(ji,jj-1,jk,Kmm)vv(ji,jj-1,jk,Kmm), 0._wp ) ) &
391	& * r1_e1e2t(ji,jj) &
392	& ) * z1_e3t
393	END_3D
394	ENDIF
395	CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1.0_wp )
396	!
397	CALL iom_put("Courant",Cu_adv)
398	!
399	IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere
400	DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 ) ! or scan Courant criterion and partition ! w where necessary
401	!
402	zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
403	! alt:
404	! IF ( ww(ji,jj,jk) > 0._wp ) THEN
405	! zCu = Cu_adv(ji,jj,jk)
406	! ELSE
407	! zCu = Cu_adv(ji,jj,jk-1)
408	! ENDIF
409	!
410	IF( zCu <= Cu_min ) THEN !<-- Fully explicit
411	zcff = 0._wp
412	ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit
413	zcff = ( zCu - Cu_min )**2
414	zcff = zcff / ( Fcu + zcff )
415	ELSE !<-- Mostly implicit
416	zcff = ( zCu - Cu_max )/ zCu
417	ENDIF
418	zcff = MIN(1._wp, zcff)
419	!
420	wi(ji,jj,jk) = zcff * ww(ji,jj,jk)
421	ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
422	!
423	Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl
424	END_3D
425	Cu_adv(:,:,1) = 0._wp
426	ELSE
427	! Fully explicit everywhere
428	Cu_adv(:,:,:) = 0._wp ! Reuse array to output coefficient below and in stp_ctl
429	wi (:,:,:) = 0._wp
430	ENDIF
431	CALL iom_put("wimp",wi)
432	CALL iom_put("wi_cff",Cu_adv)
433	CALL iom_put("wexp",ww)
434	!
435	IF( ln_timing ) CALL timing_stop('wAimp')
436	!
437	END SUBROUTINE wAimp
438
439
440	SUBROUTINE ssh_init_rst( Kbb, Kmm, Kaa )
441	!!---------------------------------------------------------------------
442	!! * ROUTINE ssh_init_rst *
443	!!
444	!! ** Purpose : ssh initialization of the sea surface height (ssh)
445	!!
446	!! ** Method : set ssh from restart or read configuration, or user_def
447	!! * ln_rstart = T
448	!! USE of IOM library to read ssh in the restart file
449	!! Leap-Frog: Kbb and Kmm are read except for l_1st_euler=T
450	!!
451	!! * otherwise
452	!! call user defined ssh or
453	!! set to -ssh_ref in wet and drying case with domcfg.nc
454	!!
455	!! NB: ssh_b/n are written by restart.F90
456	!!----------------------------------------------------------------------
457	INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices
458	!
459	INTEGER :: ji, jj, jk
460	!!----------------------------------------------------------------------
461	!
462	IF(lwp) THEN
463	WRITE(numout,*)
464	WRITE(numout,*) 'ssh_init_rst : ssh initialization'
465	WRITE(numout,*) '~~~~~~~~~~~~ '
466	ENDIF
467	!
468	! !=============================!
469	IF( ln_rstart ) THEN !== Read the restart file ==!
470	! !=============================!
471	!
472	! !* Read ssh at Kmm
473	IF(lwp) WRITE(numout,*)
474	IF(lwp) WRITE(numout,*) ' Kmm sea surface height read in the restart file'
475	CALL iom_get( numror, jpdom_auto, 'sshn' , ssh(:,:,Kmm), ldxios = lrxios )
476	!
477	IF( l_1st_euler ) THEN !* Euler at first time-step
478	IF(lwp) WRITE(numout,*)
479	IF(lwp) WRITE(numout,*) ' Euler first time step : ssh(Kbb) = ssh(Kmm)'
480	ssh(:,:,Kbb) = ssh(:,:,Kmm)
481	!
482	ELSE !* read ssh at Kbb
483	IF(lwp) WRITE(numout,*)
484	IF(lwp) WRITE(numout,*) ' Kbb sea surface height read in the restart file'
485	CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb), ldxios = lrxios )
486	ENDIF
487	! !============================!
488	ELSE !== Initialize at "rest" ==!
489	! !============================!
490	!
491	IF(lwp) WRITE(numout,*)
492	IF(lwp) WRITE(numout,*) ' initialization at rest'
493	!
494	IF( ll_wd ) THEN !* wet and dry
495	!
496	IF( ln_read_cfg ) THEN ! read configuration : ssh_ref is read in domain_cfg file
497	!!st why ssh is not masked : i.e. ssh(:,:,Kmm) = -ssh_ref*ssmask(:,:),
498	!!st since at the 1st time step lbclnk will be applied on ssh at Kaa but not initially at Kbb and Kmm
499	ssh(:,:,Kbb) = -ssh_ref
500	!
501	DO_2D( 1, 1, 1, 1 )
502	IF( ht_0(ji,jj)-ssh_ref < rn_wdmin1 ) THEN ! if total depth is less than min depth
503	ssh(ji,jj,Kbb) = rn_wdmin1 - ht_0(ji,jj)
504	ENDIF
505	END_2D
506	ELSE ! user define configuration case
507	CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) )
508	ENDIF
509	!
510	ELSE !* user defined configuration
511	CALL usr_def_istate_ssh( tmask, ssh(:,:,Kbb) )
512	!
513	ENDIF
514	!
515	ssh(:,:,Kmm) = ssh(:,:,Kbb) !* set now values from to before ones
516	ssh(:,:,Kaa) = 0._wp
517	ENDIF
518	!
519	END SUBROUTINE ssh_init_rst
520
521	!!======================================================================
522	END MODULE sshwzv

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

Download in other formats: