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 @ 12644

Last change on this file since 12644 was 12644, checked in by techene, 4 years ago
stepLF: add e3 substitute and remove pe3, domqe: new routines without e3 computation
File size: 27.4 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	! update I/O and calendar
110	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
111	indic = 0 ! reset to no error condition
112
113	IF( kstp == nit000 ) THEN ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
114	CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including passible AGRIF zoom)
115	IF( lk_diamlr ) CALL dia_mlr_iom_init ! with additional setup for multiple-linear-regression analysis
116	CALL iom_init_closedef
117	IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" ) ! for coarse grid
118	ENDIF
119	IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
120	CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
121	IF( ln_crs ) CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" ) ! tell IOM we are at time step kstp
122
123	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
124	! Update external forcing (tides, open boundaries, ice shelf interaction and surface boundary condition (including sea-ice)
125	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
126	IF( ln_tide ) CALL tide_update( kstp ) ! update tide potential
127	IF( ln_apr_dyn ) CALL sbc_apr ( kstp ) ! atmospheric pressure (NB: call before bdy_dta which needs ssh_ib)
128	IF( ln_bdy ) CALL bdy_dta ( kstp, Nnn ) ! update dynamic & tracer data at open boundaries
129	IF( ln_isf ) CALL isf_stp ( kstp, Nnn )
130	CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition (including sea-ice)
131
132	! !!st !!!!!!!!!!!!!!!!!!!!!!!
133	!
134	! emp = 0._wp
135	! emp_b = 0._wp
136	! qns = 0._wp
137	! qsr = 0._wp
138	! qns_b = 0._wp
139	!
140	! !!st
141
142	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
143	! Update stochastic parameters and random T/S fluctuations
144	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
145	IF( ln_sto_eos ) CALL sto_par( kstp ) ! Stochastic parameters
146	IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn) ) ! Random T/S fluctuations
147
148	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
149	! Ocean physics update
150	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
151	! THERMODYNAMICS
152	CALL eos_rab( ts(:,:,:,:,Nbb), rab_b, Nnn ) ! before local thermal/haline expension ratio at T-points
153	CALL eos_rab( ts(:,:,:,:,Nnn), rab_n, Nnn ) ! now local thermal/haline expension ratio at T-points
154	CALL bn2 ( ts(:,:,:,:,Nbb), rab_b, rn2b, Nnn ) ! before Brunt-Vaisala frequency
155	CALL bn2 ( ts(:,:,:,:,Nnn), rab_n, rn2, Nnn ) ! now Brunt-Vaisala frequency
156
157	! VERTICAL PHYSICS
158	CALL zdf_phy( kstp, Nbb, Nnn, Nrhs ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
159
160	! LATERAL PHYSICS
161	!
162	IF( l_ldfslp ) THEN ! slope of lateral mixing
163	CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) ) ! before in situ density
164
165	IF( ln_zps .AND. .NOT. ln_isfcav) &
166	& CALL zps_hde ( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, & ! Partial steps: before horizontal gradient
167	& rhd, gru , grv ) ! of t, s, rd at the last ocean level
168
169	IF( ln_zps .AND. ln_isfcav) &
170	& CALL zps_hde_isf( kstp, Nnn, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv, gtui, gtvi, & ! Partial steps for top cell (ISF)
171	& rhd, gru , grv , grui, grvi ) ! of t, s, rd at the first ocean level
172	IF( ln_traldf_triad ) THEN
173	CALL ldf_slp_triad( kstp, Nbb, Nnn ) ! before slope for triad operator
174	ELSE
175	CALL ldf_slp ( kstp, rhd, rn2b, Nbb, Nnn ) ! before slope for standard operator
176	ENDIF
177	ENDIF
178	! ! eddy diffusivity coeff.
179	IF( l_ldftra_time .OR. l_ldfeiv_time ) CALL ldf_tra( kstp, Nbb, Nnn ) ! and/or eiv coeff.
180	IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff.
181
182	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
183	! Ocean dynamics : hdiv, ssh, e3, u, v, w
184	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
185
186	CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor)
187	IF( .NOT.ln_linssh ) CALL dom_qe_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) )
188	IF( .NOT.ln_linssh ) CALL dom_h_nxt ( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors
189	!IF( .NOT.ln_linssh ) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa ) ! after vertical scale factors
190	CALL wzv ( kstp, Nbb, Nnn, ww, Naa ) ! now cross-level velocity
191	IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning
192	CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) ) ! now in situ density for hpg computation
193
194
195	uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
196	vv(:,:,:,Nrhs) = 0._wp
197
198	IF( lk_asminc .AND. ln_asmiau .AND. ln_dyninc ) &
199	& CALL dyn_asm_inc ( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! apply dynamics assimilation increment
200	IF( ln_bdy ) CALL bdy_dyn3d_dmp ( kstp, Nbb, uu, vv, Nrhs ) ! bdy damping trends
201	#if defined key_agrif
202	IF(.NOT. Agrif_Root()) &
203	& CALL Agrif_Sponge_dyn ! momentum sponge
204	#endif
205	CALL dyn_adv( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
206	CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS
207	CALL dyn_ldf( kstp, Nbb, Nnn , uu, vv, Nrhs ) ! lateral mixing
208	IF( ln_zdfosm ) CALL dyn_osm( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes ==> RHS
209	CALL dyn_hpg( kstp, Nnn , uu, vv, Nrhs ) ! horizontal gradient of Hydrostatic pressure
210	CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa ) ! surface pressure gradient
211
212	! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well
213	IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated
214	CALL div_hor ( kstp, Nbb, Nnn ) ! Horizontal divergence (2nd call in time-split case)
215	IF(.NOT.ln_linssh) CALL dom_qe_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
216	!IF(.NOT.ln_linssh) CALL dom_qe_sf_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 ) ! after vertical scale factors (update depth average component)
217	IF(.NOT.ln_linssh) CALL dom_h_nxt ( kstp, Nbb, Nnn, Naa, kcall=2 ) ! after vertical scale factors (update depth average component)
218	ENDIF
219	CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion
220	IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated
221	CALL wzv ( kstp, Nbb, Nnn, ww, Naa ) ! now cross-level velocity
222	IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning
223	ENDIF
224
225
226	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
227	! cool skin
228	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
229	IF ( ln_diurnal ) CALL diurnal_layers( kstp )
230
231	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
232	! diagnostics and outputs
233	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
234	IF( ln_floats ) CALL flo_stp ( kstp, Nbb, Nnn ) ! drifting Floats
235	IF( ln_diacfl ) CALL dia_cfl ( kstp, Nnn ) ! Courant number diagnostics
236	CALL dia_hth ( kstp, Nnn ) ! Thermocline depth (20 degres isotherm depth)
237	IF( ln_diadct ) CALL dia_dct ( kstp, Nnn ) ! Transports
238	CALL dia_ar5 ( kstp, Nnn ) ! ar5 diag
239	CALL dia_ptr ( kstp, Nnn ) ! Poleward adv/ldf TRansports diagnostics
240	CALL dia_wri ( kstp, Nnn ) ! ocean model: outputs
241	IF( ln_crs ) CALL crs_fld ( kstp, Nnn ) ! ocean model: online field coarsening & output
242	IF( lk_diadetide ) CALL dia_detide( kstp ) ! Weights computation for daily detiding of model diagnostics
243	IF( lk_diamlr ) CALL dia_mlr ! Update time used in multiple-linear-regression analysis
244
245	#if defined key_top
246	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
247	! Passive Tracer Model
248	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
249	CALL trc_stp ( kstp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping
250	#endif
251
252	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
253	! Active tracers
254	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255	ts(:,:,:,:,Nrhs) = 0._wp ! set tracer trends to zero
256
257	IF( lk_asminc .AND. ln_asmiau .AND. &
258	& ln_trainc ) CALL tra_asm_inc( kstp, Nbb, Nnn, ts, Nrhs ) ! apply tracer assimilation increment
259	CALL tra_sbc ( kstp, Nnn, ts, Nrhs ) ! surface boundary condition
260	IF( ln_traqsr ) CALL tra_qsr ( kstp, Nnn, ts, Nrhs ) ! penetrative solar radiation qsr
261	IF( ln_isf ) CALL tra_isf ( kstp, Nnn, ts, Nrhs ) ! ice shelf heat flux
262	IF( ln_trabbc ) CALL tra_bbc ( kstp, Nnn, ts, Nrhs ) ! bottom heat flux
263	IF( ln_trabbl ) CALL tra_bbl ( kstp, Nbb, Nnn, ts, Nrhs ) ! advective (and/or diffusive) bottom boundary layer scheme
264	IF( ln_tradmp ) CALL tra_dmp ( kstp, Nbb, Nnn, ts, Nrhs ) ! internal damping trends
265	IF( ln_bdy ) CALL bdy_tra_dmp( kstp, Nbb, ts, Nrhs ) ! bdy damping trends
266	#if defined key_agrif
267	IF(.NOT. Agrif_Root()) &
268	& CALL Agrif_Sponge_tra ! tracers sponge
269	#endif
270	CALL tra_adv ( kstp, Nbb, Nnn, ts, Nrhs ) ! hor. + vert. advection ==> RHS
271	IF( ln_zdfosm ) CALL tra_osm ( kstp, Nnn, ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS
272	IF( lrst_oce .AND. ln_zdfosm ) &
273	& CALL osm_rst ( kstp, Nnn, 'WRITE' ) ! write OSMOSIS outputs + ww (so must do here) to restarts
274	CALL tra_ldf ( kstp, Nbb, Nnn, ts, Nrhs ) ! lateral mixing
275
276	CALL tra_zdf ( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing and after tracer fields
277	IF( ln_zdfnpc ) CALL tra_npc ( kstp, Nnn, Nrhs, ts, Naa ) ! update after fields by non-penetrative convection
278
279	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
280	! Set boundary conditions, time filter and swap time levels
281	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
282	!!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap
283	!! (and time filtering) after Agrif update. Then restart would be done after and would contain updated fields.
284	!! If so:
285	!! (i) no need to call agrif update at initialization time
286	!! (ii) no need to update "before" fields
287	!!
288	!! Apart from creating new tra_swp/dyn_swp routines, this however:
289	!! (i) makes boundary conditions at initialization time computed from updated fields which is not the case between
290	!! two restarts => restartability issue. One can circumvent this, maybe, by assuming "interface separation",
291	!! e.g. a shift of the feedback interface inside child domain.
292	!! (ii) requires that all restart outputs of updated variables by agrif (e.g. passive tracers/tke/barotropic arrays) are done at the same
293	!! place.
294	!!
295	!!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine
296	CALL zdyn_ts ( Nnn, Naa, uu, vv ) ! barotrope ajustment
297	CALL finalize_sbc ( kstp, Nbb, Naa, uu, vv, ts ) ! boundary condifions
298	CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height
299	CALL dom_qe_r3c ( ssh(:,:,Nnn), r3t_f, r3u_f, r3v_f ) ! "now" ssh/h_0 ratio from filtrered ssh
300	CALL tra_atf_qco ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays
301	CALL dyn_atf_qco ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors
302	r3t(:,:,Nnn) = r3t_f(:,:)
303	r3u(:,:,Nnn) = r3u_f(:,:)
304	r3v(:,:,Nnn) = r3v_f(:,:)
305
306	!
307	! Swap time levels
308	Nrhs = Nbb
309	Nbb = Nnn
310	Nnn = Naa
311	Naa = Nrhs
312	!
313	!IF(.NOT.ln_linssh) CALL dom_qe_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors
314	IF(.NOT.ln_linssh) CALL dom_h_update ( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors
315	!
316	IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics
317
318	!!gm : This does not only concern the dynamics ==>>> add a new title
319	!!gm2: why ouput restart before AGRIF update?
320	!!
321	!!jc: That would be better, but see comment above
322	!!
323	IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nnn ) ! write output ocean restart file
324	IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters
325
326	#if defined key_agrif
327	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
328	! AGRIF
329	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
330	Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
331	CALL Agrif_Integrate_ChildGrids( stplf ) ! allows to finish all the Child Grids before updating
332
333	IF( Agrif_NbStepint() == 0 ) THEN
334	CALL Agrif_update_all( ) ! Update all components
335	ENDIF
336	#endif
337	IF( ln_diaobs ) CALL dia_obs ( kstp, Nnn ) ! obs-minus-model (assimilation) diagnostics (call after dynamics update)
338
339	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
340	! Control
341	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
342	CALL stp_ctl ( kstp, Nbb, Nnn, indic )
343
344	IF( kstp == nit000 ) THEN ! 1st time step only
345	CALL iom_close( numror ) ! close input ocean restart file
346	IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
347	IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
348	ENDIF
349
350	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
351	! Coupled mode
352	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
353	!!gm why lk_oasis and not lk_cpl ????
354	IF( lk_oasis ) CALL sbc_cpl_snd( kstp, Nbb, Nnn ) ! coupled mode : field exchanges
355	!
356	#if defined key_iomput
357	IF( kstp == nitend .OR. indic < 0 ) THEN
358	CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF
359	IF(lrxios) CALL iom_context_finalize( crxios_context )
360	IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) !
361	ENDIF
362	#endif
363	!
364	IF( ln_timing ) CALL timing_stop('stplf')
365	!
366	END SUBROUTINE stplf
367
368
369	SUBROUTINE zdyn_ts (Kmm, Kaa, puu, pvv)
370	!!----------------------------------------------------------------------
371	!! * ROUTINE zdyn_ts *
372	!!
373	!! ** Purpose : Finalize after horizontal velocity.
374	!!
375	!! ** Method : * Ensure after velocities transport matches time splitting
376	!! estimate (ln_dynspg_ts=T)
377	!!
378	!! ** Action : puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa) now and after horizontal velocity
379	!!----------------------------------------------------------------------
380	INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices
381	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities
382	!
383	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve
384	!
385	INTEGER :: jk ! dummy loop indices
386	!!----------------------------------------------------------------------
387
388	IF ( ln_dynspg_ts ) THEN
389	ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) )
390	! Ensure below that barotropic velocities match time splitting estimate
391	! Compute actual transport and replace it with ts estimate at "after" time step
392	zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1)
393	zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1)
394	DO jk = 2, jpkm1
395	zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk)
396	zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk)
397	END DO
398	DO jk = 1, jpkm1
399	puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)
400	pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)
401	END DO
402	!
403	IF( .NOT.ln_bt_fw ) THEN
404	! Remove advective velocity from "now velocities"
405	! prior to asselin filtering
406	! In the forward case, this is done below after asselin filtering
407	! so that asselin contribution is removed at the same time
408	DO jk = 1, jpkm1
409	puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )umask(:,:,jk)
410	pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )vmask(:,:,jk)
411	END DO
412	ENDIF
413	!
414	DEALLOCATE( zue, zve )
415	!
416	ENDIF
417	!
418	END SUBROUTINE zdyn_ts
419
420
421	SUBROUTINE finalize_sbc (kt, Kbb, Kaa, puu, pvv, pts)
422	!!----------------------------------------------------------------------
423	!! * ROUTINE finalize_sbc *
424	!!
425	!! ** Purpose : Apply the boundary condition on the after velocity
426	!!
427	!! ** Method : * Apply lateral boundary conditions on after velocity
428	!! at the local domain boundaries through lbc_lnk call,
429	!! at the one-way open boundaries (ln_bdy=T),
430	!! at the AGRIF zoom boundaries (lk_agrif=T)
431	!!
432	!! ** Action : puu(Kaa),pvv(Kaa) after horizontal velocity and tracers
433	!!----------------------------------------------------------------------
434	INTEGER , INTENT(in ) :: kt ! ocean time-step index
435	INTEGER , INTENT(in ) :: Kbb, Kaa ! before and after time level indices
436	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered
437	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers
438	!
439	! Update after tracer and velocity on domain lateral boundaries
440	!
441	#if defined key_agrif
442	CALL Agrif_tra ! AGRIF zoom boundaries
443	CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries
444	#endif
445	! ! local domain boundaries (T-point, unchanged sign)
446	CALL lbc_lnk_multi( 'finalize_sbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. &
447	& , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) !* local domain boundaries
448	!
449	! !* BDY open boundaries
450	IF( ln_bdy ) THEN
451	CALL bdy_tra( kt, Kbb, pts, Kaa )
452	IF( ln_dynspg_exp ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa )
453	IF( ln_dynspg_ts ) CALL bdy_dyn( kt, Kbb, puu, pvv, Kaa, dyn3d_only=.true. )
454	ENDIF
455	!
456	END SUBROUTINE finalize_sbc
457
458	!
459	!!======================================================================
460	END MODULE steplf

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

Download in other formats: