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.

stepLF.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/stepLF.F90 @ 12724

Last change on this file since 12724 was 12724, checked in by techene, 4 years ago
branch KERNEL-06 : merge with trunk@12698 #2385 - in duplcated files : changes to comply to the new trunk variables and some loop bug fixes
File size: 27.9 KB

Line
1	MODULE steplf
2	!!======================================================================
3	!! * MODULE step *
4	!! Time-stepping : manager of the ocean, tracer and ice time stepping
5	!!======================================================================
6	!! History : OPA ! 1991-03 (G. Madec) Original code
7	!! - ! 1991-11 (G. Madec)
8	!! - ! 1992-06 (M. Imbard) add a first output record
9	!! - ! 1996-04 (G. Madec) introduction of dynspg
10	!! - ! 1996-04 (M.A. Foujols) introduction of passive tracer
11	!! 8.0 ! 1997-06 (G. Madec) new architecture of call
12	!! 8.2 ! 1997-06 (G. Madec, M. Imbard, G. Roullet) free surface
13	!! - ! 1999-02 (G. Madec, N. Grima) hpg implicit
14	!! - ! 2000-07 (J-M Molines, M. Imbard) Open Bondary Conditions
15	!! NEMO 1.0 ! 2002-06 (G. Madec) free form, suppress macro-tasking
16	!! - ! 2004-08 (C. Talandier) New trends organization
17	!! - ! 2005-01 (C. Ethe) Add the KPP closure scheme
18	!! - ! 2005-11 (G. Madec) Reorganisation of tra and dyn calls
19	!! - ! 2006-01 (L. Debreu, C. Mazauric) Agrif implementation
20	!! - ! 2006-07 (S. Masson) restart using iom
21	!! 3.2 ! 2009-02 (G. Madec, R. Benshila) reintroduicing z*-coordinate
22	!! - ! 2009-06 (S. Masson, G. Madec) TKE restart compatible with key_cpl
23	!! 3.3 ! 2010-05 (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
24	!! - ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
25	!! 3.4 ! 2011-04 (G. Madec, C. Ethe) Merge of dtatem and dtasal
26	!! 3.6 ! 2012-07 (J. Simeon, G. Madec. C. Ethe) Online coarsening of outputs
27	!! 3.6 ! 2014-04 (F. Roquet, G. Madec) New equations of state
28	!! 3.6 ! 2014-10 (E. Clementi, P. Oddo) Add Qiao vertical mixing in case of waves
29	!! 3.7 ! 2014-10 (G. Madec) LDF simplication
30	!! - ! 2014-12 (G. Madec) remove KPP scheme
31	!! - ! 2015-11 (J. Chanut) free surface simplification (remove filtered free surface)
32	!! 4.0 ! 2017-05 (G. Madec) introduction of the vertical physics manager (zdfphy)
33	!! 4.1 ! 2019-08 (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
34	!!----------------------------------------------------------------------
35
36	!!----------------------------------------------------------------------
37	!! stplf : OPA system time-stepping
38	!!----------------------------------------------------------------------
39	USE step_oce ! time stepping definition modules
40	!
41	USE iom ! xIOs server
42	USE domqe
43	USE traatfqco ! time filtering (tra_atf_qco routine)
44	USE dynatfqco ! time filtering (dyn_atf_qco routine)
45	USE dynspg_ts ! surface pressure gradient: split-explicit scheme (define un_adv)
46	USE bdydyn ! ocean open boundary conditions (define bdy_dyn)
47
48	IMPLICIT NONE
49	PRIVATE
50
51	PUBLIC stplf ! called by nemogcm.F90
52
53	!!----------------------------------------------------------------------
54	!! time level indices
55	!!----------------------------------------------------------------------
56	INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !! used by nemo_init
57	# include "domzgr_substitute.h90"
58	!!----------------------------------------------------------------------
59	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
60	!! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
61	!! Software governed by the CeCILL license (see ./LICENSE)
62	!!----------------------------------------------------------------------
63	CONTAINS
64
65	#if defined key_agrif
66	RECURSIVE SUBROUTINE stplf( )
67	INTEGER :: kstp ! ocean time-step index
68	#else
69	SUBROUTINE stplf( kstp )
70	INTEGER, INTENT(in) :: kstp ! ocean time-step index
71	#endif
72	!!----------------------------------------------------------------------
73	!! * ROUTINE stplf *
74	!!
75	!! ** Purpose : - Time stepping of OPA (momentum and active tracer eqs.)
76	!! - Time stepping of SI3 (dynamic and thermodynamic eqs.)
77	!! - Time stepping of TRC (passive tracer eqs.)
78	!!
79	!! ** Method : -1- Update forcings and data
80	!! -2- Update ocean physics
81	!! -3- Compute the t and s trends
82	!! -4- Update t and s
83	!! -5- Compute the momentum trends
84	!! -6- Update the horizontal velocity
85	!! -7- Compute the diagnostics variables (rd,N2, hdiv,w)
86	!! -8- Outputs and diagnostics
87	!!----------------------------------------------------------------------
88	INTEGER :: ji, jj, jk ! dummy loop indice
89	INTEGER :: indic ! error indicator if < 0
90	!!gm kcall can be removed, I guess
91	INTEGER :: kcall ! optional integer argument (dom_vvl_sf_nxt)
92	!! ---------------------------------------------------------------------
93	#if defined key_agrif
94	kstp = nit000 + Agrif_Nb_Step()
95	Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
96	IF( lk_agrif_debug ) THEN
97	IF( Agrif_Root() .and. lwp) WRITE(,) '---'
98	IF(lwp) WRITE(,) 'Grid Number', Agrif_Fixed(),' time step ', kstp, 'int tstep', Agrif_NbStepint()
99	ENDIF
100	IF( kstp == nit000 + 1 ) lk_agrif_fstep = .FALSE.
101	# if defined key_iomput
102	IF( Agrif_Nbstepint() == 0 ) CALL iom_swap( cxios_context )
103	# endif
104	#endif
105	!
106	IF( ln_timing ) CALL timing_start('stplf')
107	!
108	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
109	! model timestep
110	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
111	!
112	IF( l_1st_euler ) THEN
113	! start or restart with Euler 1st time-step
114	rDt = rn_Dt
115	r1_Dt = 1._wp / rDt
116	ENDIF
117	!
118	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
119
120	! update I/O and calendar
121	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
122	indic = 0 ! reset to no error condition
123
124	IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
125	CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including passible AGRIF zoom)
126	IF( lk_diamlr ) CALL dia_mlr_iom_init ! with additional setup for multiple-linear-regression analysis
127	CALL iom_init_closedef
128	IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid
129	ENDIF
130	IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
131	CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
132	IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp
133
134	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
135	! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
136	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
137	IF( ln_tide ) CALL tide_update( kstp ) ! update tide potential
138	IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
139	IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries
140	IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
141	CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice)
142
143	! !!st !!!!!!!!!!!!!!!!!!!!!!!
144	!
145	! emp = 0._wp
146	! emp_b = 0._wp
147	! qns = 0._wp
148	! qsr = 0._wp
149	! qns_b = 0._wp
150	!
151	! !!st
152
153	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
154	! Update stochastic parameters and random T/S fluctuations
155	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
156	IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters
157	IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn) ) ! Random T/S fluctuations
158
159	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
160	! Ocean physics update
161	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
162	! THERMODYNAMICS
163	CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points
164	CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points
165	CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
166	CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency
167
168	! VERTICAL PHYSICS
169	CALL zdf_phy( kstp, Nbb, Nnn, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
170
171	! LATERAL PHYSICS
172	!
173	IF( l_ldfslp ) THEN ! slope of lateral mixing
174	CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) ) ! before in situ density
175
176	IF( ln_zps .AND. .NOT. ln_isfcav) &
177	& CALL zps_hde ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient
178	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
179
180	IF( ln_zps .AND. ln_isfcav) &
181	& CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
182	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
183	IF( ln_traldf_triad ) THEN
184	CALL ldf_slp_triad( kstp, Nbb, Nnn ) ! before slope for triad operator
185	ELSE
186	CALL ldf_slp ( kstp, rhd, rn2b, Nbb, Nnn ) ! before slope for standard operator
187	ENDIF
188	ENDIF
189	! ! eddy diffusivity coeff.
190	IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp, Nbb, Nnn ) ! and/or eiv coeff.
191	IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff.
192
193	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
194	! Ocean dynamics : hdiv, ssh, e3, u, v, w
195	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
196
197	CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor)
198	IF( .NOT.ln_linssh ) CALL dom_qe_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )
199	IF( .NOT.ln_linssh ) CALL dom_h_nxt ( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors
200	!IF( .NOT.ln_linssh ) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors
201	CALL wzv ( kstp, Nbb, Nnn, Naa, ww ) ! Nnn cross-level velocity
202	IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning
203	CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) ) ! now in situ density for hpg computation
204
205
206	uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
207	vv(:,:,:,Nrhs) = 0._wp
208
209	IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) &
210	& CALL dyn_asm_inc ( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! apply dynamics assimilation increment
211	IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs ) ! bdy damping trends
212	#if defined key_agrif
213	IF(.NOT. Agrif_Root()) &
214	& CALL Agrif_Sponge_dyn ! momentum sponge
215	#endif
216	CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
217	CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS
218	CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing
219	IF( ln_zdfosm ) CALL dyn_osm( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes ==> RHS
220	CALL dyn_hpg( kstp, Nnn , uu, vv, Nrhs ) ! horizontal gradient of Hydrostatic pressure
221	CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa ) ! surface pressure gradient
222
223	! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
224	IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated
225	CALL div_hor ( kstp, Nbb, Nnn ) ! Horizontal divergence (2nd call in time-split case)
226	IF(.NOT.ln_linssh) CALL dom_qe_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
227	!IF(.NOT.ln_linssh) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 ) ! after vertical scale factors (update depth average component)
228	IF(.NOT.ln_linssh) CALL dom_h_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 ) ! after vertical scale factors (update depth average component)
229	ENDIF
230	CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion
231	IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated
232	CALL wzv ( kstp, Nbb, Nnn, Naa, ww ) ! Nnn cross-level velocity
233	IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning
234	ENDIF
235
236
237	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
238	! cool skin
239	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240	IF ( ln_diurnal ) CALL diurnal_layers( kstp )
241
242	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
243	! diagnostics and outputs
244	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
245	IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats
246	IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics
247	CALL dia_hth ( kstp, Nnn ) ! Thermocline depth (20 degres isotherm depth)
248	IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports
249	CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag
250	CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics
251	CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs
252	IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output
253	IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics
254	IF( lk_diamlr ) CALL dia_mlr ! Update time used in multiple-linear-regression analysis
255
256	#if defined key_top
257	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
258	! Passive Tracer Model
259	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
260	CALL trc_stp ( kstp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping
261	#endif
262
263	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
264	! Active tracers
265	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
266	ts(:,:,:,:,Nrhs) = 0._wp ! set tracer trends to zero
267
268	IF( lk_asminc .AND. ln_asmiau .AND. &
269	& ln_trainc ) CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs ) ! apply tracer assimilation increment
270	CALL tra_sbc ( kstp, Nnn, ts, Nrhs ) ! surface boundary condition
271	IF( ln_traqsr ) CALL tra_qsr ( kstp, Nnn, ts, Nrhs ) ! penetrative solar radiation qsr
272	IF( ln_isf ) CALL tra_isf ( kstp, Nnn, ts, Nrhs ) ! ice shelf heat flux
273	IF( ln_trabbc ) CALL tra_bbc ( kstp, Nnn, ts, Nrhs ) ! bottom heat flux
274	IF( ln_trabbl ) CALL tra_bbl ( kstp, Nbb, Nnn, ts, Nrhs ) ! advective (and/or diffusive) bottom boundary layer scheme
275	IF( ln_tradmp ) CALL tra_dmp ( kstp, Nbb, Nnn, ts, Nrhs ) ! internal damping trends
276	IF( ln_bdy ) CALL bdy_tra_dmp( kstp, Nbb, ts, Nrhs ) ! bdy damping trends
277	#if defined key_agrif
278	IF(.NOT. Agrif_Root()) &
279	& CALL Agrif_Sponge_tra ! tracers sponge
280	#endif
281	CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS
282	IF( ln_zdfosm ) CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS
283	IF( lrst_oce .AND. ln_zdfosm ) &
284	& CALL osm_rst ( kstp, Nnn, 'WRITE' ) ! write OSMOSIS outputs + ww (so must do here) to restarts
285	CALL tra_ldf ( kstp, Nbb, Nnn, ts, Nrhs ) ! lateral mixing
286
287	CALL tra_zdf ( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing and after tracer fields
288	IF( ln_zdfnpc ) CALL tra_npc ( kstp, Nnn, Nrhs, ts, Naa ) ! update after fields by non-penetrative convection
289
290	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
291	! Set boundary conditions, time filter and swap time levels
292	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
293	!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
294	!! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
295	!! If so:
296	!! (i) no need to call agrif update at initialization time
297	!! (ii) no need to update "before" fields
298	!!
299	!! Apart from creating new tra_swp/dyn_swp routines, this however:
300	!! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
301	!! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
302	!! e.g. a shift of the feedback interface inside child domain.
303	!! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
304	!! place.
305	!!
306	CALL zdyn_ts ( Nnn, Naa, uu, vv ) ! barotrope ajustment
307	CALL finalize_sbc ( kstp, Nbb, Naa, uu, vv, ts ) ! boundary condifions
308	CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height
309	CALL dom_qe_r3c ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh
310	CALL tra_atf_qco ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays
311	CALL dyn_atf_qco ( kstp, Nbb, Nnn, Naa, uu, vv ) ! time filtering of "now" velocities and scale factors
312	r3t(:,:,Nnn) = r3t_f(:,:)
313	r3u(:,:,Nnn) = r3u_f(:,:)
314	r3v(:,:,Nnn) = r3v_f(:,:)
315
316	!
317	! Swap time levels
318	Nrhs = Nbb
319	Nbb = Nnn
320	Nnn = Naa
321	Naa = Nrhs
322	!
323	!IF(.NOT.ln_linssh) CALL dom_qe_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors
324	IF(.NOT.ln_linssh) CALL dom_h_update ( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors
325	!
326	IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics
327
328	!!gm : This does not only concern the dynamics ==>>> add a new title
329	!!gm2: why ouput restart before AGRIF update?
330	!!
331	!!jc: That would be better, but see comment above
332	!!
333	IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file
334	IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters
335
336	#if defined key_agrif
337	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
338	! AGRIF
339	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
340	Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
341	CALL Agrif_Integrate_ChildGrids( stplf ) ! allows to finish all the Child Grids before updating
342
343	IF( Agrif_NbStepint() == 0 ) THEN
344	CALL Agrif_update_all( ) ! Update all components
345	ENDIF
346	#endif
347	IF( ln_diaobs ) CALL dia_obs ( kstp, Nnn ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
348
349	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
350	! Control
351	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
352	CALL stp_ctl ( kstp, Nbb, Nnn, indic )
353
354	IF( kstp == nit000 ) THEN ! 1st time step only
355	CALL iom_close( numror ) ! close input ocean restart file
356	IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
357	IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
358	ENDIF
359
360	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
361	! Coupled mode
362	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
363	!!gm why lk_oasis and not lk_cpl ????
364	IF( lk_oasis ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
365	!
366	#if defined key_iomput
367	IF( kstp == nitend .OR. indic < 0 ) THEN
368	CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
369	IF(lrxios) CALL iom_context_finalize( crxios_context )
370	IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
371	ENDIF
372	#endif
373	!
374	IF( l_1st_euler ) THEN ! recover Leap-frog timestep
375	rDt = 2._wp * rn_Dt
376	r1_Dt = 1._wp / rDt
377	l_1st_euler = .FALSE.
378	ENDIF
379	!
380	IF( ln_timing ) CALL timing_stop('stplf')
381	!
382	END SUBROUTINE stplf
383
384
385	SUBROUTINE zdyn_ts (Kmm, Kaa, puu, pvv)
386	!!----------------------------------------------------------------------
387	!! * ROUTINE zdyn_ts *
388	!!
389	!! ** Purpose : Finalize after horizontal velocity.
390	!!
391	!! ** Method : * Ensure after velocities transport matches time splitting
392	!! estimate (ln_dynspg_ts=T)
393	!!
394	!! ** Action : puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa) now and after horizontal velocity
395	!!----------------------------------------------------------------------
396	INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices
397	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities
398	!
399	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve
400	!
401	INTEGER :: jk ! dummy loop indices
402	!!----------------------------------------------------------------------
403
404	IF ( ln_dynspg_ts ) THEN
405	ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) )
406	! Ensure below that barotropic velocities match time splitting estimate
407	! Compute actual transport and replace it with ts estimate at "after" time step
408	zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
409	zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
410	DO jk = 2, jpkm1
411	zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
412	zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
413	END DO
414	DO jk = 1, jpkm1
415	puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
416	pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
417	END DO
418	!
419	IF( .NOT.ln_bt_fw ) THEN
420	! Remove advective velocity from "now velocities"
421	! prior to asselin filtering
422	! In the forward case, this is done below after asselin filtering
423	! so that asselin contribution is removed at the same time
424	DO jk = 1, jpkm1
425	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )umask(:,:,jk)
426	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )vmask(:,:,jk)
427	END DO
428	ENDIF
429	!
430	DEALLOCATE( zue, zve )
431	!
432	ENDIF
433	!
434	END SUBROUTINE zdyn_ts
435
436
437	SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts)
438	!!----------------------------------------------------------------------
439	!! * ROUTINE finalize_sbc *
440	!!
441	!! ** Purpose : Apply the boundary condition on the after velocity
442	!!
443	!! ** Method : * Apply lateral boundary conditions on after velocity
444	!! at the local domain boundaries through lbc_lnk call,
445	!! at the one-way open boundaries (ln_bdy=T),
446	!! at the AGRIF zoom boundaries (lk_agrif=T)
447	!!
448	!! ** Action : puu(Kaa),pvv(Kaa) after horizontal velocity and tracers
449	!!----------------------------------------------------------------------
450	INTEGER , INTENT(in ) :: kt ! ocean time-step index
451	INTEGER , INTENT(in ) :: Kbb, Kaa ! before and after time level indices
452	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered
453	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers
454	!
455	! Update after tracer and velocity on domain lateral boundaries
456	!
457	#if defined key_agrif
458	CALL Agrif_tra ! AGRIF zoom boundaries
459	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
460	#endif
461	! ! local domain boundaries (T-point, unchanged sign)
462	CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. &
463	& , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) !* local domain boundaries
464	!
465	! !* BDY open boundaries
466	IF( ln_bdy ) THEN
467	CALL bdy_tra( kt, Kbb, pts, Kaa )
468	IF( ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
469	IF( ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
470	ENDIF
471	!
472	END SUBROUTINE finalize_sbc
473
474	!
475	!!======================================================================
476	END MODULE steplf

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

Download in other formats: